Interneuronal mechanisms of hippocampal theta oscillations in a full-scale model of the rodent CA1 circuit

  1. Marianne J Bezaire  Is a corresponding author
  2. Ivan Raikov
  3. Kelly Burk
  4. Dhrumil Vyas
  5. Ivan Soltesz
  1. University of California, Irvine, United States
  2. Stanford University, United States

Abstract

The hippocampal theta rhythm plays important roles in information processing; however, the mechanisms of its generation are not well understood. We developed a data-driven, supercomputer-based, full-scale (1:1) model of the rodent CA1 area and studied its interneurons during theta oscillations. Theta rhythm with phase-locked gamma oscillations and phase-preferential discharges of distinct interneuronal types spontaneously emerged from the isolated CA1 circuit without rhythmic inputs. Perturbation experiments identified parvalbumin-expressing interneurons and neurogliaform cells, as well as interneuronal diversity itself, as important factors in theta generation. These simulations reveal new insights into the spatiotemporal organization of the CA1 circuit during theta oscillations.

https://doi.org/10.7554/eLife.18566.001

Introduction

The hippocampal CA1 area supports diverse cognitive tasks including learning, memory, and spatial processing (Squire, 1992; Remondes and Schuman, 2004; Manns et al., 2007; Moser et al., 2008). These cognitive tasks are thought to require coordination of neuronal activity provided by physiological network oscillations, including the theta rhythm (Buzsáki, 2002; Buzsáki and Moser, 2013). In rodents, hippocampal theta is a 5–10 Hz oscillation in the local field potential (LFP) and neuronal firing probabilities (Soltesz and Deschênes, 1993; Lee et al., 1994; Ylinen et al., 1995; Klausberger and Somogyi, 2008; Varga et al., 2012, 2014), occurring during locomotion and in REM sleep (Buzsáki, 2002). Though several major afferents provide theta-frequency rhythmic input to the CA1 in vivo (Soltesz and Deschênes, 1993; Buzsáki, 2002; Fuhrmann et al., 2015), recent reports indicate the presence of spontaneous theta-frequency LFP oscillations even in the isolated whole CA1 preparation in vitro (Goutagny et al., 2009; Amilhon et al., 2015). Therefore, the latter studies suggest an intrinsic ability of the CA1 circuit to generate some form of theta waves even without rhythmic external inputs. However, the intra-CA1 mechanisms that may contribute to the generation of the theta rhythm are not well understood (Colgin, 2013, 2016).

Here we investigated the ability of the CA1 to generate intrinsic theta oscillations using a uniquely biological data-driven, full-scale computer model of the isolated CA1 network. Recent advances in supercomputing power and high-quality synaptic connectivity data present the intriguing opportunity to develop full-scale models where every biological synapse and neuron is explicitly represented. In principle, such full-scale models of mammalian circuits comprising hundreds of thousands of neurons of distinct types advantageously avoid the connectivity scaling tradeoff that besets reduced-scale models: smaller models of large networks with realistic single cell electrophysiological properties (e.g., input resistance and resting membrane potential) remain silent unless synaptic strengths or numbers are arbitrarily increased beyond the biologically relevant levels to compensate for fewer inputs to their model cells (e.g., [Dyhrfjeld-Johnsen et al., 2007; Sterratt et al., 2011]). Biological relevance may also increase as other network components are modeled in greater detail. However, full-scale models require considerable computational resources. Further, such detailed models have a large parameter space which risks being sub-optimally constrained by neurobiological properties that are only partially quantified (Sejnowski et al., 1988). Because the CA1 area is one of the most extensively studied brain regions, there are abundant anatomical and electrophysiological data about its organization, making it a logical choice for the development of a full-scale model. The CA1 area is also worth modeling at full-scale because of the diverse cognitive tasks it supports. These tasks likely require the simultaneous processing of thousands of incoming and outgoing signals, and full-scale network models, at least in principle, have the potential to match this in vivo processing capacity.

In this paper, we describe the development of a full-scale CA1 computational network model of unprecedented biological detail and its application to gain insights into the roles and temporal organization of CA1 interneurons during theta rhythm. The simulated full-scale CA1 circuit was able to spontaneously generate theta waves as well as phase-locked gamma oscillations. Furthermore, distinct interneuron types discharged at particular phases of theta, demonstrating that phase-preferential firing (Klausberger et al., 2003, 20042005Ferraguti et al., 2005; Jinno et al., 2007; Fuentealba et al., 2008; Klausberger and Somogyi, 2008; Varga et al., 2012; Lapray et al., 2012; Katona et al., 2014; Varga et al., 2014) originates in part within the CA1 network. Perturbation experiments revealed that parvalbumin-expressing (PV+) interneurons, neurogliaform cells, connections between CA1 pyramidal cells, and interneuronal diversity were important for theta generation. These results provide new mechanistic insights into the emergence of the theta rhythm from within the CA1 circuitry and the role of interneurons in theta oscillations.

Results

Development of a data-driven, full-scale model of the isolated CA1

Details of the full-scale model are described in the Methods, and the most important features are illustrated in Figures 1 and 2 and summarized here. Briefly, CA1 model cells were evenly distributed within their respective layers in a 3-dimensional prism with realistic dimensions for the rodent hippocampal CA1 region (Figure 1A and B). The model network contained 338,740 cells (similar to the biological CA1 in rats, including 311,500 pyramidal cells and 27,240 interneurons) (Figure 1D–E and Figure 1—figure supplement 1). In addition, the network also incorporated 454,700 artificial stimulating cells (spiking units with random, Poisson-distributed interspike intervals) to simulate afferents to CA1; the cell type-specific distribution, dendritic position, amplitude and kinetics of the excitatory input synapses were all experimentally constrained by afferent CA3 and entorhinal cortical data. Cell type-specific connectivity data, including cell numbers (Figure 1D) and convergence and divergence values (Figure 1E; Figure 1—figure supplement 1 and Table 1) were taken without alteration from our previously published, in-depth, quantitative assessment of the CA1 circuit (Bezaire and Soltesz, 2013). Anatomical constraints of the connectivity were implemented in the model by accounting for the distribution of the axonal boutons as a function of longitudinal and transverse distance from the presynaptic cell soma (Figure 1—figure supplement 2). The afferent divergence and convergence onto the cells were also anatomically patterned, maintaining the topographical arrangement seen experimentally (Hongo et al., 2015), for a total of 5.19 billion synaptic connections in the model network. In addition, the remaining parameters that could not be constrained by experimental data were documented, with the assumptions used to arrive at them explicitly listed in Table 2 of Bezaire and Soltesz (2013) and additional parameter calculations described in the Appendix of the present paper, section 'Inhibitory connectivity'. To highlight the many constraints applied in the current work and address the unconstrained model parameters, we characterized all model components (constrained and unconstrained) in experimental terms, comparing with experimental data where possible (Figure 2; Appendix). For a four second simulation, the full-scale model required 3–4 terabytes (TB) of RAM and four hours of execution time on a supercomputer using ~3000 processors (or up to 12 hr for simulations calculating a high-accuracy local field potential (LFP) analog). Additional details and data about model performance are available in Table 2 and Bezaire et al. (2016a).

Figure 1 with 2 supplements see all
CA1 network connectivity.

(A) The model network is arranged in a layered prism with the lengths of each dimension similar to the actual dimensions of the CA1 region and its layers. (B) The model cell somata within a small chunk of stratum pyramidale (as depicted in A) are plotted to show the regular distribution of model cells throughout the layer in which they are found. (C) Each pyramidal cell in the network has detailed morphology with realistic incoming synapse placement along the dendrites and soma. (D,E) Diagrams illustrate connectivity between types of cells. (D) The network includes one principal cell type (pyramidal cells) and eight interneuron types. Cell types that may connect are linked by a line colored according to the presynaptic cell type. Most cell types can connect to most other cell types. Total number of cells of each type are displayed, as are the number of local output synapses (boutons) from all cells of each type. (E) The number, position, and cell types of each connection are biologically constrained, as are the numbers and positions of the cells. See Figure 1—figure supplement 1) for details about the convergence onto each cell type. Also see Table 1 and Figure 1—figure supplement 2 for information about the cell-type combinations of the 5 billion connections and the axonal distributions followed by each cell type, as well as detailed connectivity results at http://doi.org/10.6080/K05H7D60.

https://doi.org/10.7554/eLife.18566.002
Electrophysiology of the model network components.

(A) Ion channel densities vary as a function of location (top) in the morphologically detailed pyramidal cell model (bottom; adapted from Poolos et al., 2002). Scale bar: 100 μm and 0.01 μF/cm2. (B–C) The sodium channel found in the pyramidal cell soma is characterized in terms of (B) the activation/inactivation curves and (C) the current-voltage relation at peak (transient) current and steady state. (DG) Current sweeps are shown for four model cell types: (D) PV+ basket cell, (E) CCK+ basket cell, (F) O-LM cell, and (G) neurogliaform cell. Scale bar: 100 ms and 20 mV. (H–J) Electrophysiological properties for each cell type, including (H) input resistance, (I) membrane time constant, and (J) action potential threshold. (K–L) Pyramidal cell synaptic connections are characterized as post-synaptic currents with the postsynaptic cell voltage clamped at −50 mV; (K) synapses made onto the pyramidal cell from all other cell types and (L) synapses made by the pyramidal cell onto all network cell types. Cells represented by same colors as in Figure 1. Source Data available for electrophysiological characterizations shown here. Additional details available in the Methods, Table 3, and the Appendix.

https://doi.org/10.7554/eLife.18566.005
Figure 2—source data 1

Model sodium channel activation.

The ion channel characterized in this figure was an Nav channel, inserted into a single compartment cell of diameter and length 16.8 microns (a soma) with a density such that the maximum, macroscopic conductance was 0.001 μS/cm2. The reversal potential of the channel was + 55 mV and the settings during the characterization protocol were: temperature = 34 degrees Celsius, axial resistance = 210 ohm*cm, [Ca2+]internal=5.0000e-06 mM, specific membrane capacitance = 1 μF/cm2. For activation steps, the cell was held at −120 mV and then stepped to potential levels ranging from −90 mV to + 90 mV. For inactivation steps, the cell was held at various potential levels ranging from −90 mV to + 40 mV for 500 ms and then stepped to + 20 mV. Each current injection step is recorded in three columns, where t: time (ms), i: injection (nA), g: conductance (S/cm2). The column labels are followed by the voltage (hold or step, according to the file), with activation steps being recorded in the Na_Channel_Step.dat and inactivation steps being recorded in the Na_Channel_Hold.dat file.

https://doi.org/10.7554/eLife.18566.006
Figure 2—source data 2

Model sodium channel inactivation.

The ion channel characterized in this figure was an Nav channel, inserted into a single compartment cell of diameter and length 16.8 microns (a soma) with a density such that the maximum, macroscopic conductance was .001 μS/cm2. The reversal potential of the channel was + 55 mV and the settings during the characterization protocol were: temperature = 34 degrees Celsius, axial resistance = 210 ohm*cm, [Ca2+]internal=5.0000e-06 mM, specific membrane capacitance = 1 μF/cm2. For activation steps, the cell was held at −120 mV and then stepped to potential levels ranging from −90 mV to +90 mV. For inactivation steps, the cell was held at various potential levels ranging from −90 mV to +40 mV for 500 ms and then stepped to +20 mV. Each current injection step is recorded in three columns, where t: time (ms), i: injection (nA), g: conductance (S/cm2). The column labels are followed by the voltage (hold or step, according to the file), with activation steps being recorded in the Na_Channel_Step.dat and inactivation steps being recorded in the Na_Channel_Hold.dat file.

https://doi.org/10.7554/eLife.18566.007
Figure 2—source data 3

Model axo-axonic cell current injection sweep.

Simulated current injection sweep in AxoClamp ATF (tab-delimited) file format.

https://doi.org/10.7554/eLife.18566.008
Figure 2—source data 4

Model bistratified cell current injection sweep.

Simulated current injection sweep in AxoClamp ATF (tab-delimited) file format.

https://doi.org/10.7554/eLife.18566.009
Figure 2—source data 5

Model CCK+ basket cell current injection sweep.

Simulated current injection sweep in AxoClamp ATF (tab-delimited) file format.

https://doi.org/10.7554/eLife.18566.010
Figure 2—source data 6

Model ivy cell current injection sweep.

Simulated current injection sweep in AxoClamp ATF (tab-delimited) file format.

https://doi.org/10.7554/eLife.18566.011
Figure 2—source data 7

Model neurogliaform cell current injection sweep.

Simulated current injection sweep in AxoClamp ATF (tab-delimited) file format.

https://doi.org/10.7554/eLife.18566.012
Figure 2—source data 8

Model O-LM cell current injection sweep.

Simulated current injection sweep in AxoClamp ATF (tab-delimited) file format.

https://doi.org/10.7554/eLife.18566.013
Figure 2—source data 9

Model PV+ basket cell current injection sweep.

Simulated current injection sweep in AxoClamp ATF (tab-delimited) file format.

https://doi.org/10.7554/eLife.18566.014
Figure 2—source data 10

Model pyramidal cell current injection sweep.

Simulated current injection sweep in AxoClamp ATF (tab-delimited) file format.

https://doi.org/10.7554/eLife.18566.015
Figure 2—source data 11

Model Schaffer Collateral-Associated cell current injection sweep.

Simulated current injection sweep in AxoClamp ATF (tab-delimited) file format.

https://doi.org/10.7554/eLife.18566.016
Figure 2—source data 12

Model paired recording of an Axo-axonic cell to Pyramidal cell connection.

Simulated paired recordings where the postsynaptic cell was voltage-clamped at −50 mV and the reversal potential of the synapse was kept at its physiological potential, as defined in the network model code. Sodium channels in the postsynaptic cell were blocked to prevent a suprathreshold response. A spike was triggered in the presynaptic cell and the current response was measured in the postsynaptic cell at the soma. This recording was repeated 10 times, with a randomly chosen connection location (from anatomically likely locations) each time. Each of the 10 trials are included in this file.

https://doi.org/10.7554/eLife.18566.017
Figure 2—source data 13

Model paired recording of a Bistratified cell to Pyramidal cell connection.

Simulated paired recordings where the postsynaptic cell was voltage-clamped at −50 mV and the reversal potential of the synapse was kept at its physiological potential, as defined in the network model code. Sodium channels in the postsynaptic cell were blocked to prevent a suprathreshold response. A spike was triggered in the presynaptic cell and the current response was measured in the postsynaptic cell at the soma. This recording was repeated 10 times, with a randomly chosen connection location (from anatomically likely locations) each time. Each of the 10 trials are included in this file.

https://doi.org/10.7554/eLife.18566.018
Figure 2—source data 14

Model paired recording of a CA3 cell to Pyramidal cell connection.

Simulated paired recordings where the postsynaptic cell was voltage-clamped at −50 mV and the reversal potential of the synapse was kept at its physiological potential, as defined in the network model code. Sodium channels in the postsynaptic cell were blocked to prevent a suprathreshold response. A spike was triggered in the presynaptic cell and the current response was measured in the postsynaptic cell at the soma. This recording was repeated 10 times, with a randomly chosen connection location (from anatomically likely locations) each time. Each of the 10 trials are included in this file.

https://doi.org/10.7554/eLife.18566.019
Figure 2—source data 15

Model paired recording of a CCK+ basket cell to Pyramidal cell connection.

Simulated paired recordings where the postsynaptic cell was voltage-clamped at −50 mV and the reversal potential of the synapse was kept at its physiological potential, as defined in the network model code. Sodium channels in the postsynaptic cell were blocked to prevent a suprathreshold response. A spike was triggered in the presynaptic cell and the current response was measured in the postsynaptic cell at the soma. This recording was repeated 10 times, with a randomly chosen connection location (from anatomically likely locations) each time. Each of the 10 trials are included in this file.

https://doi.org/10.7554/eLife.18566.020
Figure 2—source data 16

Model paired recording of an ECIII cell to Pyramidal cell connection.

Simulated paired recordings where the postsynaptic cell was voltage-clamped at −50 mV and the reversal potential of the synapse was kept at its physiological potential, as defined in the network model code. Sodium channels in the postsynaptic cell were blocked to prevent a suprathreshold response. A spike was triggered in the presynaptic cell and the current response was measured in the postsynaptic cell at the soma. This recording was repeated 10 times, with a randomly chosen connection location (from anatomically likely locations) each time. Each of the 10 trials are included in this file.

https://doi.org/10.7554/eLife.18566.021
Figure 2—source data 17

Model paired recording of an Ivy cell to Pyramidal cell connection.

Simulated paired recordings where the postsynaptic cell was voltage-clamped at −50 mV and the reversal potential of the synapse was kept at its physiological potential, as defined in the network model code. Sodium channels in the postsynaptic cell were blocked to prevent a suprathreshold response. A spike was triggered in the presynaptic cell and the current response was measured in the postsynaptic cell at the soma. This recording was repeated 10 times, with a randomly chosen connection location (from anatomically likely locations) each time. Each of the 10 trials are included in this file.

https://doi.org/10.7554/eLife.18566.022
Figure 2—source data 18

Model paired recording of a Pyramidal cell to Pyramidal cell connection.

Simulated paired recordings where the postsynaptic cell was voltage-clamped at −50 mV and the reversal potential of the synapse was kept at its physiological potential, as defined in the network model code. Sodium channels in the postsynaptic cell were blocked to prevent a suprathreshold response. A spike was triggered in the presynaptic cell and the current response was measured in the postsynaptic cell at the soma. This recording was repeated 10 times, with a randomly chosen connection location (from anatomically likely locations) each time. Each of the 10 trials are included in this file.

https://doi.org/10.7554/eLife.18566.023
Figure 2—source data 19

Model paired recording of a Neurogliaform cell to Pyramidal cell connection.

Simulated paired recordings where the postsynaptic cell was voltage-clamped at −50 mV and the reversal potential of the synapse was kept at its physiological potential, as defined in the network model code. Sodium channels in the postsynaptic cell were blocked to prevent a suprathreshold response. A spike was triggered in the presynaptic cell and the current response was measured in the postsynaptic cell at the soma. This recording was repeated 10 times, with a randomly chosen connection location (from anatomically likely locations) each time. Each of the 10 trials are included in this file.

https://doi.org/10.7554/eLife.18566.024
Figure 2—source data 20

Model paired recording of an O-LM cell to Pyramidal cell connection.

Simulated paired recordings where the postsynaptic cell was voltage-clamped at −50 mV and the reversal potential of the synapse was kept at its physiological potential, as defined in the network model code. Sodium channels in the postsynaptic cell were blocked to prevent a suprathreshold response. A spike was triggered in the presynaptic cell and the current response was measured in the postsynaptic cell at the soma. This recording was repeated 10 times, with a randomly chosen connection location (from anatomically likely locations) each time. Each of the 10 trials are included in this file.

https://doi.org/10.7554/eLife.18566.025
Figure 2—source data 21

Model paired recording of a PV+ basket cell to Pyramidal cell connection.

Simulated paired recordings where the postsynaptic cell was voltage-clamped at −50 mV and the reversal potential of the synapse was kept at its physiological potential, as defined in the network model code. Sodium channels in the postsynaptic cell were blocked to prevent a suprathreshold response. A spike was triggered in the presynaptic cell and the current response was measured in the postsynaptic cell at the soma. This recording was repeated 10 times, with a randomly chosen connection location (from anatomically likely locations) each time. Each of the 10 trials are included in this file.

https://doi.org/10.7554/eLife.18566.026
Figure 2—source data 22

Model paired recording of a Pyramidal cell to Axo-axonic cell connection.

Simulated paired recordings where the postsynaptic cell was voltage-clamped at −50 mV and the reversal potential of the synapse was kept at its physiological potential, as defined in the network model code. Sodium channels in the postsynaptic cell were blocked to prevent a suprathreshold response. A spike was triggered in the presynaptic cell and the current response was measured in the postsynaptic cell at the soma. This recording was repeated 10 times, with a randomly chosen connection location (from anatomically likely locations) each time. Each of the 10 trials are included in this file.

https://doi.org/10.7554/eLife.18566.027
Figure 2—source data 23

Model paired recording of a Pyramidal cell to Bistratified cell connection.

Simulated paired recordings where the postsynaptic cell was voltage-clamped at −50 mV and the reversal potential of the synapse was kept at its physiological potential, as defined in the network model code. Sodium channels in the postsynaptic cell were blocked to prevent a suprathreshold response. A spike was triggered in the presynaptic cell and the current response was measured in the postsynaptic cell at the soma. This recording was repeated 10 times, with a randomly chosen connection location (from anatomically likely locations) each time. Each of the 10 trials are included in this file.

https://doi.org/10.7554/eLife.18566.028
Figure 2—source data 24

Model paired recording of a Pyramidal cell to Ivy cell connection.

Simulated paired recordings where the postsynaptic cell was voltage-clamped at −50 mV and the reversal potential of the synapse was kept at its physiological potential, as defined in the network model code. Sodium channels in the postsynaptic cell were blocked to prevent a suprathreshold response. A spike was triggered in the presynaptic cell and the current response was measured in the postsynaptic cell at the soma. This recording was repeated 10 times, with a randomly chosen connection location (from anatomically likely locations) each time. Each of the 10 trials are included in this file.

https://doi.org/10.7554/eLife.18566.029
Figure 2—source data 25

Model paired recording of a Pyramidal cell to O-LM cell connection.

Simulated paired recordings where the postsynaptic cell was voltage-clamped at −50 mV and the reversal potential of the synapse was kept at its physiological potential, as defined in the network model code. Sodium channels in the postsynaptic cell were blocked to prevent a suprathreshold response. A spike was triggered in the presynaptic cell and the current response was measured in the postsynaptic cell at the soma. This recording was repeated 10 times, with a randomly chosen connection location (from anatomically likely locations) each time. Each of the 10 trials are included in this file.

https://doi.org/10.7554/eLife.18566.030
Figure 2—source data 26

Model paired recording of a Pyramidal cell to PV+ basket cell connection.

Simulated paired recordings where the postsynaptic cell was voltage-clamped at −50 mV and the reversal potential of the synapse was kept at its physiological potential, as defined in the network model code. Sodium channels in the postsynaptic cell were blocked to prevent a suprathreshold response. A spike was triggered in the presynaptic cell and the current response was measured in the postsynaptic cell at the soma. This recording was repeated 10 times, with a randomly chosen connection location (from anatomically likely locations) each time. Each of the 10 trials are included in this file.

https://doi.org/10.7554/eLife.18566.031
Figure 2—source data 27

Model paired recording of a Pyramidal cell to Schaffer Collateral-Associated cell connection.

Simulated paired recordings where the postsynaptic cell was voltage-clamped at −50 mV and the reversal potential of the synapse was kept at its physiological potential, as defined in the network model code. Sodium channels in the postsynaptic cell were blocked to prevent a suprathreshold response. A spike was triggered in the presynaptic cell and the current response was measured in the postsynaptic cell at the soma. This recording was repeated 10 times, with a randomly chosen connection location (from anatomically likely locations) each time. Each of the 10 trials are included in this file.

https://doi.org/10.7554/eLife.18566.032
Table 1

Number of synapses between each cell type. Connections between cells generally comprise 1–10 synapses each. Presynaptic cells are listed down the first column (corresponding to each row) and postsynaptic cells are listed along the first row (corresponding to each column).

https://doi.org/10.7554/eLife.18566.033
Pre/PostAxoBisCCK+BIvyNGFO-LMPyrPV+BSC-A
Axo0.00e + 000.00e + 000.00e + 000.00e + 000.00e + 000.00e + 001.12e + 070.00e + 000.00e + 00
Bis2.35e + 053.54e + 055.76e + 052.64e + 050.00e + 006.40e + 053.12e + 078.85e + 056.80e + 04
CCK+B1.41e + 052.12e + 059.79e + 055.64e + 050.00e + 002.62e + 053.24e + 075.31e + 058.32e + 04
Ivy3.53e + 055.30e + 053.42e + 062.11e + 061.00e + 062.23e + 061.28e + 081.33e + 064.08e + 05
NGF0.00e + 000.00e + 000.00e + 000.00e + 006.09e + 050.00e + 004.36e + 070.00e + 000.00e + 00
O-LM1.18e + 051.77e + 051.44e + 060.00e + 004.65e + 059.84e + 042.49e + 074.42e + 051.60e + 05
Pyr7.19e + 052.43e + 060.00e + 002.38e + 050.00e + 001.17e + 076.14e + 077.03e + 061.26e + 05
PV+B5.73e + 048.62e + 041.37e + 057.05e + 040.00e + 000.00e + 005.83e + 072.16e + 059.60e + 03
SC-A8.82e + 031.33e + 041.30e + 051.06e + 050.00e + 001.97e + 043.74e + 063.32e + 041.44e + 04
CA31.23e + 072.56e + 071.44e + 073.39e + 070.00e + 000.00e + 003.73e + 096.69e + 071.55e + 06
ECIII1.43e + 061.91e + 064.02e + 060.00e + 003.75e + 060.00e + 008.09e + 080.00e + 004.58e + 05
Table 2

Simulation time, exchange time, and load balance for simulations executed on various supercomputers and numbers of processors.

https://doi.org/10.7554/eLife.18566.034
Supercomputer# ProcessorsSim time (s)Exchange time (s)Load balance
Comet16802610.281.050.999
Comet17042566.760.650.999
Comet17282601.220.860.999
Comet via NSG17282060.880.830.999
Stampede via NSG20482471.641.711.000
Stampede20482578.320.291.000
Stampede25282189.561.780.999
Stampede30081844.220.910.999
Stampede34881641.910.860.999

An important set of constraints was the electrophysiology and other properties of individual cells and synapses (Figure 2; Figure 2—source data 327; Tables 3 and 4) that were based on experimental data (Lee et al., 2016; Quattrocolo and Maccaferri, 2016). Briefly, our pyramidal cell model (Poolos et al., 2002) contained 200 compartments in a realistic morphology and six fully characterized ion channel types with kinetics and densities based on anatomical location within the cell (Figure 2A–C; Figure 2—source data 12). We included eight model interneuron types (Klausberger and Somogyi, 2008; Soltesz, 2006; Armstrong and Soltesz, 2012): PV+ basket cells (these fast-spiking cells synapse on the somata and proximal dendrites of CA1 pyramidal cells), cholecystokinin+ (CCK+) basket cells (these regular-spiking cells also innervate the somata and proximal dendrites, but have properties and functions distinct from the PV+ basket cells), bistratified cells (these PV+ and somatostatin+ (SOM+) fast-spiking cells innervate the basal and apical dendritic trees), axo-axonic cells (these PV+ fast-spiking cells synapse only on the axon initial segments of pyramidal cells and are also known as chandelier cells), Schaffer Collateral-Associated (SC-A) cells (these CCK+, regular-spiking cells innervate dendrites in the stratum radiatum), oriens-lacunosum-moleculare (O-LM) cells (these SOM+ cells project to the distal dendrites in the stratum lacunosum-moleculare though their somata are located in the stratum oriens), neurogliaform cells (these cells have relatively small dendrites and a dense axonal cloud, and they innervate distal dendrites in the stratum lacunosum-moleculare), and ivy cells (these cells are similar to neurogliaform cells, but innervate proximal dendrites) (Figure 2D–E). Some interneurons in the model, as in the biological network, also innervated other interneurons (Table 1). For greater detail of model connectivity, including convergence per single cell, synaptic amplitude, and other factors, see the Appendix. These cell types collectively comprise the majority (~70%) of known CA1 interneurons (Bezaire and Soltesz, 2013). The remaining 30% of the interneurons were not included in the model due to paucity of quantitative data (Bezaire and Soltesz, 2013). We differentiated the interneurons by their electrophysiological profiles, connectivity patterns, synaptic properties, and anatomical abundance (Gulyás et al., 1991; Hájos and Mody, 1997; Maccaferri et al., 2000; Megías et al., 2001; Lee et al., 2010; Krook-Magnuson et al., 2011; Bezaire and Soltesz, 2013; Lee et al., 2014). The synaptic connections were implemented using double exponential mechanisms to better fit experimental data on rise and decay time constants. We used experimental data to constrain the synaptic kinetics, amplitudes, and locations on the postsynaptic cell (Figure 1E, 2K and L). We implemented the model in parallel NEURON (Carnevale and Hines, 2005) and executed the simulations on several supercomputers. All model results, characterizations, and experimental comparisons are publicly available.

Table 3

Electrophysiological characteristics of each model cell type. For more information about model electrophysiology, see the Appendix.

https://doi.org/10.7554/eLife.18566.035
ConditionPyrPV+BCCK+BSC-AAxoBisO-LMIvyNGF
Resting Membrane Potential (mV)−63.0−65.0−70.6−70.5−65.0−67.0−71.5−60.0−60.0
Input Resistance (MΩ)62.252.0211.0272.452.098.7343.8100.0100.0
Membrane Tau (ms)4.86.922.624.47.014.722.421.121.1
Rheobase (pA)250.0300.060.040.0200.0350.050.0160.0170.0
Threshold (mV)52.0−36.6−40.6−43.1−41.6−28.1100.2−27.6−27.7
Delay to 1st Spike (ms)12.474.6166.6127.743.528.48.9173.3119.0
Half-Width (ms)80.70.91.91.60.60.5112.90.60.6
Table 4

Current injection levels used to characterize interneuron current sweeps in Figure 2D–G.

https://doi.org/10.7554/eLife.18566.036
Cell typeHyper. (pA)Step size (pA)Depol. (pA)
PV+ B.−30050+500
CCK+ B.−10020+80
O-LM−13030+80
NGF−13020+190

Emergence of spontaneous theta and gamma oscillations in the full-scale model in the absence of rhythmic external inputs

First, we examined whether the well-constrained, biologically detailed, full-scale CA1 model could oscillate spontaneously within the physiological range. Based on reports of spontaneous theta-frequency LFP oscillations in the isolated CA1 preparation (Goutagny et al., 2009), we expected a sufficiently constrained CA1 model to generate spontaneous theta rhythm when given tonic, arrhythmic excitation. We varied the magnitude of arrhythmic, tonic excitation to the network (by systematically changing the mean spiking frequency of the artificial stimulating cells, see above) and identified excitation levels where the network developed a stable, spontaneous theta rhythm (5–10 Hz; Figure 3 and 4; Figure 3—source data 13 and Figure 4—source data 12). The pyramidal cell spikes (Figures 3C and D) exhibited peak power around the theta frequency of 7.8 Hz (Figure 4 and Table 7). Importantly, every measure of network activity showed theta oscillations, including the somatic intracellular membrane potential from individual cells (Figure 3D), the spike times of individual cells and all cells collectively (Figure 3C), and aggregate measures such as the spike density function (Szucs, 1998) per cell type and the LFP analog (Figure 3A and 4; see also Figure 4—figure supplement 1). In all of these measures of network activity, theta was apparent within one theta period of the simulation start. The theta oscillation was stable, maintaining a steady power level throughout the duration of the oscillation (Figure 4A). To our knowledge, this is the first strictly data-driven, full-scale computational network model of the CA1 that exhibits spontaneous theta rhythm without rhythmic synaptic inputs.

Detailed network activity.

(A–D) One second of network activity is shown. (A–B) The LFP analog, filtered at (A) the theta range of 5–10 Hz and (B) the low gamma range of 25–40 Hz, shows consistent theta and gamma signals. Scale bar represents 100 ms and 0.2 mV (theta) or 0.27 mV (gamma) for filtered LFP traces. (C) Raster of all spikes from cells within 100 μm of the reference electrode point. (D) Representative intracellular somatic membrane potential traces from cells near the reference electrode point. Scale bar represents 100 ms and 50 mV for the intracellular traces.

https://doi.org/10.7554/eLife.18566.037
Figure 3—source data 1

Filtered analog local field potential of model network.

The theta-filtered and gamma-filtered local field potential (LFP) analog traces.

https://doi.org/10.7554/eLife.18566.038
Figure 3—source data 2

Spike Raster.

Spike times for the length of the entire simulation, from the specific cells displayed in raster shown in Figure 3 (spike times of every single cell in the network are available in the CRCNS repository entry for Bezaire et al. [2016b]).

https://doi.org/10.7554/eLife.18566.039
Figure 3—source data 3

Somatic membrane potential recordings.

Full duration, intracellular somatic membrane potential recordings from the specific cells shown in Figure 3.

https://doi.org/10.7554/eLife.18566.040
Figure 4 with 1 supplement see all
Spectral analysis of model activity.

(A) A spectrogram of the local pyramidal-layer LFP analog (including contributions from all pyramidal cells within 100 μm of the reference electrode and 10% of pyramidal cells outside that radius) shows the stability and strength of the theta oscillation over time. The oscillation also featured strong harmonics at multiples of the theta frequency of 7.8 Hz. (B,D) Welch’s periodogram of the spike density function for each cell type, normalized by cell type and by displayed frequency range, shows the dominant network frequencies of (B) theta (7.8 Hz) and (D) gamma (71 Hz). Power is normalized to the peak power displayed in the power spectrum for each cell type. (C) Cross-frequency coupling between theta and gamma components of the LFP analog shows that the gamma oscillation is theta modulated. The gamma envelope is a function of the theta phase with the largest amplitude gamma oscillations occurring at the trough of the theta oscillation. Following convention, the theta trough was designated 0°/360°; see e.g., Varga et al. (2012). A graphical explanation of the relation between a spike train and its spike density function is shown in Figure 4—figure supplement 1.

https://doi.org/10.7554/eLife.18566.041
Figure 4—source data 1

Raw analog local field potential of model network.

The raw local field potential (LFP) analog calculated from the network activity, as detailed in the Materials and methods section.

https://doi.org/10.7554/eLife.18566.042
Figure 4—source data 2

Spike Density Functions of each cell type in control network.

The power of the Spike Density Functions was calculated from a one-sided periodogram using Welch’s method where segments have a 50% overlap with a Hamming Window.

https://doi.org/10.7554/eLife.18566.043

In addition to theta rhythm, the model network displayed gamma oscillations (25–80 Hz; Figures 3B and 4D), as expected based on in vivo data (Soltesz and Deschênes, 1993; Tort et al., 2009; Colgin and Moser, 2010) and in vitro slice data showing 65–75 Hz gamma oscillations arising in response to theta rhythmic network stimulation (Butler et al., 2016). The amplitude envelope of the gamma oscillation was phase-locked to the theta rhythm (Figures 3A,B and 4C), as it is in the biological CA1, representing cross-frequency coupling (Soltesz and Deschênes, 1993; Bragin et al., 1995; Buzsáki et al., 2003; Jensen and Colgin, 2007; Belluscio et al., 2012). The highest amplitude of the gamma oscillations in the model was observed at the theta trough (0°/360°) in the pyramidal layer LFP analog (Figure 4C). Because the current study focused primarily on theta oscillations and experimental data from the isolated CA1 are available only for the theta rhythm (Goutagny et al., 2009; Amilhon et al., 2015), the gamma oscillations were not examined further in the present study.

These results demonstrate that, in spite of gaps in our knowledge, our model was sufficiently well-constrained by experimental data that it generated theta and gamma oscillations on its own, without extrinsic rhythmic inputs or deliberate tuning of intrinsic parameters.

Although we generally refrained from deliberately compensating for missing parameters in this paper, it is of course possible to do so. For example, as mentioned above, no sufficiently detailed information was available for certain interneuron types. Therefore, these lesser-known interneurons were not included in the model, which meant that inhibition received by the pyramidal cells was probably weaker than in the biological situation. Indeed, the pyramidal cells in our model described above (Figures 3 and 4) tended to fire more than they typically do during theta oscillations in vivo (e.g., [Soltesz and Deschênes, 1993; Robbe et al., 2006]). Is the higher firing frequency of the pyramidal cells related to the weaker inhibition? To answer the latter question, in a subset of the simulations we artificially scaled up inhibition in the model to match the inhibitory synapse numbers on CA1 pyramidal cells that were expected from electron microscopic reconstructions of pyramidal cell dendrites and somata (Megías et al., 2001; Bezaire and Soltesz, 2013). The rationale for scaling up inhibition in this way was that, as described in detail in Bezaire and Soltesz (2013), the estimates of local inhibitory inputs to pyramidal cells were different when based on experimental observations of presynaptic anatomy (local boutons available for synapsing from distinct types of intracellularly filled and reconstructed interneurons) as opposed to postsynaptic anatomy (inhibitory post-synaptic densities on pyramidal cell dendrites). In simulations with the model containing this rationally scaled up inhibition, only 1% of the pyramidal cells were active, and they fired at a low rate of 1.8 Hz (data not shown), closely resembling the in vivo condition (Soltesz and Deschênes, 1993; Robbe et al., 2006). Therefore, the model was capable of reproducing the experimentally observed relatively low-firing frequencies for the principal cells during theta oscillations in vivo. However, because the source of the additional inhibition onto CA1 principal cells has not yet been experimentally identified, we used the connectivity estimates as constrained by experimental observations of axonal boutons and lengths in the full scale model (without the scaled-up inhibition) described above (Figures 3 and 4) in the subsequent computational experiments.

Mechanism of theta generation and phase-preferential firing of interneurons in the full-scale model of the isolated CA1

Next, we examined the onset of the theta rhythm and the firing patterns of the various cell types in the model circuit during theta oscillations (Figure 5, Table 5, and Figure 5—source data 111 ). As mentioned above, distinct interneuronal types, defined based on their selective axonal innervation patterns of the postsynaptic domains of pyramidal cells, exhibit characteristic, cell-type-specific preferred phases of firing during theta oscillations in vivo (Klausberger et al., 2003; 2004, 2005; Ferraguti et al., 2005; Jinno et al., 2007; Fuentealba et al., 2008; Varga et al., 2012; Lapray et al., 2012; Katona et al., 2014; Varga et al., 2014). Importantly, this fundamental property emerged spontaneously from the full-scale model, without purposeful tuning of parameters except the mean spiking frequency and synaptic strength of the artificial stimulating cells to set the incoming excitation levels from afferents (see Materials and methods for details). As expected, the numerically dominant pyramidal cells, whose intracellular membrane potential oscillations to a large extent generate and underlie the extracellular LFP signal during theta oscillations (Buzsáki et al., 2012), preferentially discharged around the trough 0o/360o of the LFP analog theta rhythm (Figure 5A).

Figure 5 with 2 supplements see all
Model and experimental cell theta phases.

All model results are based on the spiking of the cells within 100 μm of the reference electrode. (A–B) Firing probability by cell type as a function of theta phase for (A) model and (B) experimental cells under anesthesia (histograms adapted with permission from Figure 2, Figure 5B left, and Figure 6F respectively from Klausberger and Somogyi, 2008; Fuentealba et al., 2008; Fuentealba et al., 2010). The model histograms are normalized; see Figure 5—figure supplement 1 for firing rates. (C) Theta phase preference and theta modulation level were correlated; better modulated cell types spiked closer to the LFP analog trough near the phase preference of pyramidal cells. (D) Theta phase preference plotted on an idealized LFP wave for model data (base of arrow signifies the model phase preference and head of the arrow shows the distance to anesthetized, experimental phase preference).

https://doi.org/10.7554/eLife.18566.045
Figure 5—source data 1

Spike times of axo-axonic cells.

All spike times from all axo-axonic cells, as well as the calculated theta phases (relative to the theta-filtered LFP analog) of each spike.

https://doi.org/10.7554/eLife.18566.046
Figure 5—source data 2

Spike times of bistratified cells.

All spike times from all bistratified cells, as well as the calculated theta phases (relative to the theta-filtered LFP analog) of each spike.

https://doi.org/10.7554/eLife.18566.047
Figure 5—source data 3

Spike times of proximal afferent cells.

All spike times from all proximal afferent cells, as well as the calculated theta phases (relative to the theta-filtered LFP analog) of each spike.

https://doi.org/10.7554/eLife.18566.048
Figure 5—source data 4

Spike times of CCK+ basket cells.

All spike times from all CCK+ basket cells, as well as the calculated theta phases (relative to the theta-filtered LFP analog) of each spike.

https://doi.org/10.7554/eLife.18566.049
Figure 5—source data 5

Spike times of distal afferent cells.

All spike times from all distal afferent cells, as well as the calculated theta phases (relative to the theta-filtered LFP analog) of each spike.

https://doi.org/10.7554/eLife.18566.050
Figure 5—source data 6

Spike times of ivy cells.

All spike times from all ivy cells, as well as the calculated theta phases (relative to the theta-filtered LFP analog) of each spike.

https://doi.org/10.7554/eLife.18566.051
Figure 5—source data 7

Spike times of neurogliaform cells.

All spike times from all neurogliaform cells, as well as the calculated theta phases (relative to the theta-filtered LFP analog) of each spike.

https://doi.org/10.7554/eLife.18566.052
Figure 5—source data 8

Spike times of O-LM cells.

All spike times from all O-LM cells, as well as the calculated theta phases (relative to the theta-filtered LFP analog) of each spike.

https://doi.org/10.7554/eLife.18566.053
Figure 5—source data 9

Spike times of PV+ basket cells.

All spike times from all PV+ basket cells, as well as the calculated theta phases (relative to the theta-filtered LFP analog) of each spike.

https://doi.org/10.7554/eLife.18566.054
Figure 5—source data 10

Spike times of pyramidal cells.

All spike times from all pyramidal cells, as well as the calculated theta phases (relative to the theta-filtered LFP analog) of each spike.

https://doi.org/10.7554/eLife.18566.055
Figure 5—source data 11

Spike times of Schaffer Collateral-associated cells.

All spike times from all Schaffer Collateral-associated cells, as well as the calculated theta phases (relative to the theta-filtered LFP analog) of each spike.

https://doi.org/10.7554/eLife.18566.056
Table 5

Preferred theta firing phases for each model cell type.

https://doi.org/10.7554/eLife.18566.059
Cell typeFiring rate (Hz)ModulationPhase (0o=trough)
Levelp
Axo.8.90.074.58e − 130163.4
Bis.18.00.760.00e + 00340.0
CCK+ B.54.40.100.00e + 00202.8
Ivy43.30.330.00e + 00142.1
NGF.55.10.071.46e − 32176.3
O-LM17.40.760.00e + 00334.7
Pyr.6.00.740.00e + 00339.7
PV+ B.0.90.460.00e + 00356.8
S.C.-A.5.20.031.13e − 07197.9
Table 6

Firing rates and theta phase preferences for various cell types in various conditions. Theta phase is relative to the LFP recorded in the pyramidal layer, where 0o and 360o are at the trough of the oscillation. non: non-theta/non-SWR state. SWR: sharp wave/ripple. u+k and x: urethane + supplemental doses of ketamine and xylazine.

https://doi.org/10.7554/eLife.18566.060
Cell typeFiring rate (Hz)Theta phase (o)State of animalAnimalRef.
ThetaNonSWR
ADI8.600.060.25156anesth: u+k and xrat(Klausberger et al., 2005)
Axo-axonic17.103.502.95185anesth: u+k and xrat(Klausberger et al., 2003)
Axo-axonic2727251awake, head restraintmouse(Varga et al., 2014)
Bistratified5.900.9042.801anesth: u+k and xrat(Klausberger et al., 2004)
Bistratified34360awake, head restraintmouse(Varga et al., 2014)
Bistratified30.4227.6535.822awakerat(Katona et al., 2014)
CCK+ Basket9.401.602.70174anesth: u+k and xrat(Klausberger et al., 2005)
Ivy0.701.700.8031anesth: u+k and xrat(Fuentealba et al., 2008)
Ivy2.802.105.2046awake, freerat(Lapray et al., 2012)
Ivy2.403.006.70awake, freerat(Fuentealba et al., 2008)
NGF6.002.652.30196anesth: u+k and xrat(Fuentealba et al., 2010)
O-LM4.902.300.2319anesth: u+k and xrat(Klausberger et al., 2003)
O-LM29.8010.4025.40346awake, head restraintmouse(Varga et al., 2012)
O-LM17.3011.8818.95342awakerat(Katona et al., 2014)
PPA5.751.951.50100anesth: u+k and xrat(Klausberger et al., 2005)
PV+ Basket7.302.7432.68271anesth: u+k and xrat(Klausberger et al., 2003)
PV+ Basket234anesth: u+k and xrat(Klausberger et al., 2005)
PV+ Basket21.006.50122.00289awake, freerat(Lapray et al., 2012)
PV+ Basket25.008.2075.00307awake, head restraintmouse(Varga et al., 2012)
PV+ Basket2877310awake, head restraintmouse(Varga et al., 2014)
Pyramidal20anesth: u+k and xrat(Klausberger et al., 2003)
Trilaminar0.200.1069.00troughanesth: u+k and xrat(Ferraguti et al., 2005)
Double Proj.0.900.5526.9377anesth: u+k and xrat(Jinno et al., 2007)
Oriens Retro.0.530.3753.3728anesth: u+k and xrat(Jinno et al., 2007)
Radiatum Retro.5.151.900.70298anesth: u+k and xrat(Jinno et al., 2007)

Interneurons in the model preferentially fired at specific phases of theta oscillations, depending on the cell type. Their phase preferences fell into two broad categories (Figure 5A). The cells belonging to the first group, including the PV+ basket cells, bistratified cells and O-LM cells, were most likely to fire at the theta trough compared to other theta phases. Since these cells received substantial excitatory inputs from local CA1 pyramidal cells both in the biological state and in the model (Bezaire and Soltesz, 2013), their firing in the isolated CA1 model was probably driven by the pyramidal cell discharges around the theta trough. In contrast, the second group of cells, including the ivy and neurogliaform cells, the CCK+ basket cells and the axo-axonic cells, fired least around the theta trough, leading to an inverted firing probability distribution relative to the first group of interneurons (Figure 5A). Their differing phase preferences were most likely due to a combination of weak or non-existent excitatory inputs from local CA1 pyramidal cells and inhibition from the interneurons that prominently discharged around the theta trough. In general agreement with the first group of cells being strongly and rhythmically driven by the local pyramidal cells, there was a correlation between the phase preference and the strength of modulation (Figure 5C; see Materials and methods), with the cells discharging around the trough all showing strong modulation of firing.

These results were in line with recent data from the isolated CA1 preparation in vitro (Ferguson et al., 2015) which showed that cells belonging to the broadly defined SOM+ and PV+ classes (identified using genetic drivers) displayed phase preferences similar to the O-LM, PV+ basket and bistratified cells in our model (note that Ferguson and colleagues used LFP theta recorded in the stratum radiatum as reference, which is approximately 180 degrees out of phase with the pyramidal cell layer theta used in this paper). In addition, the interneuronal phase preferences in the model were also remarkably similar to in vivo data from anesthetized animals (Figure 5B; because no data are available on the phase preferential firing of morphologically identified interneurons from the isolated CA1 preparation, comparison is made here with results from anesthetized animals, from which the most complete data sets are available; see also Discussion). Specifically, the majority (71%; 5/7) of the interneuron types for which there were experimental data, including the CCK+ basket, axo-axonic, bistratified, O-LM and neurogliaform cells, showed similar preferential maxima in their firing probabilities in the model (Figure 5A) and in vivo (Figure 5B). The largest differences between the model and the in vivo phase-preferential firing occurred for the PV+ basket cells and the ivy cells, suggesting that during theta oscillations in vivo these cells may be strongly driven by CA3 afferents active during the late falling phase of the theta cycle (Colgin and Moser, 2010); note that PV+ cells receive a high number of excitatory inputs on their dendrites compared to other interneuron classes (Gulyás et al., 1999). A comparison of the model and the anesthetized in vivo data is illustrated in Figure 5D, where the arrows indicate the shift required for the model phase preferences (Figure 5A) to equal the in vivo (Figure 5B) phase preferences; note that the required shifts (arrows) are small for all interneuron types except PV+ basket and ivy cells. A clear majority of the interneuronal types in the model showed phase preferences similar to the in vivo condition where rhythmically discharging afferent inputs are present, indicating that theta-preferential discharges are to a large extent determined by the wiring properties of the CA1 circuit itself.

Perturbation experiments indicate a key role for interneuronal diversity in the emergence of spontaneous theta

Importantly, the ability to generate theta oscillations, phase-locked gamma oscillations, and theta-related phase-preferential firing of distinct interneuronal subtypes was not a universal property of the model. As shown in Figure 6A, our strongly constrained model only exhibited spontaneous theta oscillations at certain levels of afferent excitation. The results described above (Figures 35) were obtained with an afferent excitation level of 0.65 Hz (labeled as ‘Control’ in Figure 6A), meaning that each excitatory afferent cell excited the model network with a Poisson-distributed spike train having a Poisson mean interspike interval (ISI) corresponding to a firing rate of 0.65 Hz. When the excitation level decreased below 0.65 Hz, the theta rhythm fell apart, and when the excitation level increased beyond 0.80 Hz, theta power also started to drop significantly as the oscillation frequency rose out of theta range (Figure 6 and Figure 6—figure supplement 1; Figure 6—source data 12), evolving into a beta oscillation (Engel and Fries, 2010). These data indicate that while synaptic-cellular organization of the CA1 circuit enables the intrinsic, within-CA1 generation of theta waves, the circuit is predisposed to exhibit theta oscillations only under particular excitatory input conditions. The observation that, under certain conditions the model network can oscillate at frequencies between 12 and 20 Hz, is in agreement with recent experimental findings that rhythmic driving of septal PV+ cells can reliably entrain the hippocampus in a 1:1 ratio up to frequencies of 20 Hz (Dannenberg et al., 2015).

Figure 6 with 2 supplements see all
Altered network configurations.

Oscillation power (in mV22/Hz) of the spike density function (SDF) for pyramidal cells within 100 μm of the reference electrode, at the peak frequency within theta range (5–10 Hz) in altered network configurations. For corresponding peak frequencies, see Figure 6—figure supplement 1. (A) Theta is present at some excitation levels. (B) Muting each cell type’s output caused a range of effects. (C) The stability and frequency of spontaneous theta in the network was sensitive to the presence and number of recurrent connections between CA1 pyramidal cells. (D) Partially muting the broad classes of PV+ or SOM+ cells by 50% showed that PV+ muting disrupted the network more than SOM+ muting. (E) Theta falls apart when all interneurons are given the same electrophysiological profile, whether it be of a PV+ basket, CCK+ basket, neurogliaform, or O-LM cell. (F) Gradually setting all interneuron properties to those of PV+ basket cells did not restore theta. From left to right: control network; PV+ basket cell electrophysiology; also weights of incoming synapses; also numbers of incoming synapses; then all interneurons being PV+ basket cells (with the addition of the output synapse numbers, weights, and kinetics); then variable RMP (normal distribution with standard deviation of 8 mV). (G) A wide range in excitation was unable to produce theta in the PV+ B. network. (H) Removing the GABAB component from the neurogliaform synapses onto other neurogliaform cells and pyramidal cells showed a significant drop in theta power. Massively increasing the weight of the GABAA component to produce a similar amount of charge transfer restored theta power (compare the IPSCs corresponding to each condition in Figure 6—figure supplement 2). Standard deviations (n = 3) shown; significance (p=1.8e-05).

https://doi.org/10.7554/eLife.18566.061
Figure 6—source data 1

Simulation name mapping.

Map the names of the simulations (used in the header of SDF_All_Conditions.txt) to the bar labels in the graphs of Figure 6.

https://doi.org/10.7554/eLife.18566.062
Figure 6—source data 2

SDF of each network condition.

The full length pyramidal cell Spike Density Function computed at a resolution of 1000 Hz from the spikes of all pyramidal cells within the local range of the electrode point in the model network, for each network condition studied in Figure 6.

https://doi.org/10.7554/eLife.18566.063

Does the parameter sensitivity of the theta rhythm also apply to recurrent excitation from pyramidal cells and inhibition from CA1 interneurons? In order to answer the latter question, we tested whether the theta rhythm was differentially sensitive to the contribution of each inhibitory cell type (Figure 6B). We characterized the contribution of each local CA1 cell type to the theta rhythm by muting the output of the cell type so that its activity had no effect on the network. First, we studied the role of the recurrent collaterals of pyramidal cells, which contact mostly interneurons and, less frequently, other pyramidal cells (Bezaire and Soltesz, 2013). When we muted all the outputs from pyramidal cells, theta rhythm disappeared (bar labeled ‘Pyr’ in Figure 6B), indicating that the recurrent collaterals of pyramidal cells play a key role in theta oscillations.

Interestingly, muting the relatively rare CA1 pyramidal cell to pyramidal cell excitatory connections alone (each pyramidal cell contacts 197 other pyramidal cells in the CA1; Bezaire and Soltesz, 2013) was sufficient to collapse the theta rhythm (bar labeled ‘None’ in Figure 6C); key roles for inter-pyramidal cell excitatory synapses within CA1 have been suggested for sharp wave ripple oscillations as well (Maier et al., 2011). Furthermore, the parameter-sensitivity of the theta rhythm was also apparent when examining the role of pyramidal cell to pyramidal cell connections, because theta power dramatically decreased when these connections were either increased (doubled) or decreased (halved) from the biologically observed 197 (Figure 6C). Next, we investigated the effects of muting the output from each interneuron type. Silencing the output from any of the fast-spiking, PV family interneurons (PV+ basket, axo-axonic, or bistratified cells), CCK+ basket cells, or neurogliaform cells also strongly reduced theta power in the network (Figure 6B). In contrast, muting other interneuronal types (S.C.-A cells, O-LM cells, or ivy cells) had no effect on this form of theta oscillations generated by the intra-CA1 network (Figure 6B). In additional disinhibition studies simulating optogenetic experimental configurations, partial muting of all PV+ outputs (PV+ basket, bistratified, and axo-axonic cells together) had a larger effect than partial muting of all SOM+ outputs (O-LM and bistratified cells); see Figure 6D. Reassuringly, these results were in overall agreement with experimental data from the isolated CA1 preparation indicating that optogenetic silencing of PV+ cells, but not SOM+ cells such as the O-LM cells, caused a marked reduction in theta oscillations (Amilhon et al., 2015). The differential effects of silencing PV+ versus SOM+ cells could also be obtained in a rationally simplified model called the Network Clamp, where a single pyramidal cell was virtually extracted from the full-scale CA1 network with all of its afferent synapses intact (for further details, see Bezaire et al., 2016a).

Since the diverse sources of inhibition from the various interneuronal types are believed to enable networks to achieve more complex behaviors, including oscillations (Soltesz, 2006; Rotstein et al., 2005; Kepecs and Fishell, 2014), we next tested if reducing the diversity of interneurons in the model would affect its ability to produce spontaneous theta oscillations. Surprisingly, giving all interneurons a single electrophysiological profile appeared to create conditions that were not conducive for the appearance of spontaneous theta oscillations regardless of which interneuronal profile was used (Figure 6E; note that the cells still differed in the strengths, distribution, and identities of their incoming and outgoing connections after this manipulation). To probe this finding further, we focused on PV+ basket cells, which have been implicated in theta generation in vivo (Soltesz and Deschênes, 1993; Buzsáki, 2002; Stark et al., 2013; Hu et al., 2014) and exhibited strong theta power in their spiking in the control network model (Figure 4B). We gradually altered (‘morphed’) the properties of all other model interneuron types until they became PV+ basket cells, by first converging their electrophysiological profiles, then additionally their synaptic kinetics and incoming synapse weights, then also their incoming synapse numbers, and finally their outgoing synaptic weights and numbers (Figure 6F; Table 7). Theta was not apparent in any intermediate steps nor in the final network where all interneurons had become PV+ basket cells (‘All PV+B’ in Figure 6F). Furthermore, introduction of cell to cell variability in the resting membrane potential of interneurons in the ‘All PV+B’ configuration at the biologically observed values for PV+ basket cells also failed to restore theta (‘Var PV+B’ in Figure 6F shows results with standard deviation of (SD) = 8 mV in the resting membrane potential; SD = 5 mV and SD = 2 mV also yielded no theta; biological SD value: approximately 5 mV in Tricoire et al. (2011) and 2 mV in Mercer et al. (2012)). Therefore, although PV basket cells appear to be important for theta-generation both in the biological and the model CA1 network, endowing all interneurons with PV basket cell-like properties does not lead to a network configuration conducive to theta oscillations (Hendrickson et al., 2015).

Table 7

Peak, theta and gamma frequencies and powers of the pyramidal cell spike density function using Welch’s Periodogram. As in Figure 6—figure supplement 1, networks where no pyramidal cells spiked - resulting in zero power within the spectral analysis of the pyramidal cell spike density function - have their peak frequencies listed as ‘n/a’ for ‘not available’.

https://doi.org/10.7554/eLife.18566.066
ThetaGammaOverall
ConditionFrequencyPowerFrequencyPowerFrequencyPower
Tonic excitation level (Hz)
0.20n/a0.0e + 00n/a0.0e + 00n/a0.0e + 00
0.405.95.6e + 0425.44.1e + 0413.76.5e + 04
0.509.88.1e + 0425.41.0e + 0519.55.6e + 05
0.65 (Ctrl.)7.85.0e + 0525.42.0e + 057.85.0e + 05
0.809.87.8e + 0529.32.6e + 059.87.8e + 05
1.009.86.8e + 0529.31.4e + 059.86.8e + 05
1.209.85.1e + 0533.21.8e + 0511.78.2e + 05
1.409.81.9e + 0525.43.4e + 0511.78.6e + 05
Single Interneuron E’phys. Profile
Ctrl7.85.0e + 0525.42.0e + 057.85.0e + 05
O-LMn/a0.0e + 00n/a0.0e + 00n/a0.0e + 00
CCK+B9.85.7e + 0362.56.9e + 0562.56.9e + 05
PV+Bn/a0.0e + 00n/a0.0e + 00n/a0.0e + 00
NGF5.92.6e + 0439.12.4e + 0639.12.4e + 06
Inh. Cells Converge to PV+ B. Cells
Ctrl7.85.0e + 0525.42.0e + 057.85.0e + 05
E’phys.n/a0.0e + 00n/a0.0e + 00n/a0.0e + 00
+input wgt7.86.8e + 0244.91.6e + 0621.53.4e + 06
+input #9.86.1e + 0331.31.1e + 0615.62.0e + 06
All PV+Bn/a0.0e + 00n/a0.0e + 00n/a0.0e + 00
Var. PV+Bn/a0.0e + 00n/a0.0e + 00n/a0.0e + 00
Outputs Muted
Ctrl7.85.0e + 0525.42.0e + 057.85.0e + 05
SOM7.84.7e + 0527.31.4e + 057.84.7e + 05
PV9.83.2e + 0427.38.1e + 0513.71.5e + 06
Pyr to Pyr
2.0x9.81.1e + 0525.47.3e + 0513.71.0e + 06
1.0x (Ctrl.)7.85.0e + 0525.42.0e + 057.85.0e + 05
0.5x7.88.0e + 0429.32.2e + 0529.32.2e + 05
None9.81.1e + 0070.33.7e + 0170.33.7e + 01
Outputs Muted From Each Cell Type
Ctrl7.85.0e + 0525.42.0e + 057.85.0e + 05
Pyr7.81.1e + 0070.33.8e + 0170.33.8e + 01
PV+B9.88.8e + 0329.31.9e + 0629.31.9e + 06
SC-A9.84.9e + 0527.31.8e + 059.84.9e + 05
O-LM7.85.1e + 0525.48.3e + 047.85.1e + 05
NGF9.85.2e + 0327.39.1e + 0513.71.6e + 06
Ivy7.85.3e + 0525.42.0e + 057.85.3e + 05
CCK+B5.95.5e + 0325.43.3e + 033.95.7e + 03
Bis5.91.3e + 0429.31.7e + 0629.31.7e + 06
Axo7.84.0e + 0333.21.2e + 0615.61.9e + 06
Pyr & PV+ B. Network: Tonic Excitation Level (Hz)
0.01n/a0.0e + 00n/a0.0e + 00n/a0.0e + 00
0.05n/a0.0e + 00n/a0.0e + 00n/a0.0e + 00
0.10n/a0.0e + 00n/a0.0e + 00n/a0.0e + 00
0.20n/a0.0e + 00n/a0.0e + 00n/a0.0e + 00
0.405.92.3e + 0225.41.2e + 023.92.4e + 02
0.65n/a0.0e + 00n/a0.0e + 00n/a0.0e + 00
0.80n/a0.0e + 00n/a0.0e + 00n/a0.0e + 00
1.20n/a0.0e + 00n/a0.0e + 00n/a0.0e + 00
Ctrl7.85.0e + 0525.42.0e + 057.85.0e + 05

To rule out the possibility that the lack of theta could be due to an inappropriate excitation level in these reduced diversity configurations, we subjected the ‘All PV+ B’ network to a wide range of incoming excitation levels (Figure 6G). Theta rhythm did not appear at any of these excitation levels. While we could not rule out a hypothetical theta regime somewhere in the parameter space of such low-diversity configurations, any theta solution space would likely be smaller and more elusive than we were able to determine in the control configuration (Figure 6A).

Taken together, these results indicated, for the first time, that interneuronal diversity itself is an important factor in the emergence of spontaneous theta oscillations from the CA1 network.

Neurogliaform cell signaling and theta generation in the isolated CA1 model

In agreement with previous predictions (Capogna, 2011), the perturbation experiments described above suggested that neurogliaform cells were a necessary component for spontaneous theta to arise in the isolated CA1. We wondered why muting the output from neurogliaform cells, but not the closely related ivy cells, affected theta oscillations (Figure 6B), especially since there were fewer neurogliaform cells than ivy cells, and they were less theta modulated (Figure 5A). These two model interneuron groups mainly differed in that the neurogliaform cells evoked mixed GABAA,B postsynaptic events (Price et al., 2005), whereas the model ivy cells only triggered GABAA IPSPs (in agreement with a lack of evidence for ivy cell-evoked GABAB IPSPs). Could the slow kinetics of GABAB IPSPs contribute to the pacing of the theta oscillations? Indeed, when we selectively removed the GABAB component of all neurogliaform cell outgoing synaptic connections, theta power was strongly reduced (Figure 6H). To test whether the contribution of the GABAB receptors was due to their slow kinetics, we artificially sped up the GABAB IPSPs so that they had GABAA kinetics but conserved their characteristic large charge transfer. This alteration was implemented by scaling up the GABAA synaptic conductance at neurogliaform cell output synapses to achieve a similar total charge transfer as the control GABAA,B mixed synapse (Figure 6—figure supplement 2). As shown in Figure 6H (green bar), theta activity was restored when the neurogliaform cell output synapses had no slow GABAB component, only a scaled up fast GABAA IPSP with a charge transfer equivalent to the mixed GABAA,B synapses. Therefore, muting the neurogliaform cells strongly disrupted the theta oscillations not because the theta oscillations required the slow kinetics of GABAB IPSPs specifically, but because the slow kinetics enabled a large total charge transfer.

Discussion

Emergence of theta oscillations from a biological data-driven, full-scale model of the CA1 network

We produced a biologically detailed, full-scale CA1 network model constrained by extensive experimental data (Bezaire and Soltesz, 2013). When excited with arrhythmic inputs at physiologically relevant levels (see below), the model displayed spontaneous theta (and gamma) oscillations with phase preferential firing across the nine model cell types (pyramidal cells and eight interneuron classes). Consistent with experimental results (Goutagny et al., 2009; Amilhon et al., 2015), these oscillations emerged from the network model without explicit encoding, rhythmic inputs or purposeful tuning of intra-CA1 parameters (all anatomical connectivity parameters were exactly as previously published in Bezaire and Soltesz (2013)). Cell type-specific perturbations of the network showed that each interneuronal type contributed uniquely to the spontaneous theta oscillation, and that the presence of diverse inhibitory dynamics was a necessary condition for sustained theta oscillations. In addition to characterizing roles for specific network components, these model results generally suggest that the presence of diverse interneuronal types and the intrinsic circuitry of the CA1 network are sufficient and necessary to enable the isolated CA1 to oscillate at spontaneous theta rhythms while supporting distinct phase preferences of each class of hippocampal neuron. These abilities may serve to maintain the stability and robustness of the theta oscillation mechanism as it operates in vivo in diverse behavioral states. The theta rhythm is thought to be important for organizing disparate memory tasks (Lisman and Idiart, 1995; Hasselmo et al., 2002; Hasselmo, 2005; Lisman and Jensen, 2013; Siegle and Wilson, 2014), and a CA1 network which has evolved a predisposition to oscillate at theta and gamma frequencies may enable more efficient processing of the phasic input it receives in vivo (Akam and Kullmann, 2012; Fries, 2015). In turn, phase preferential firing may aid information processing tasks by providing order and allowing multiple channels of information to be processed in parallel (Jensen and Lisman, 2000; Hasselmo et al., 2002; Womelsdorf et al., 2007; Schomburg et al., 2014; Jeewajee et al., 2014; Maris et al., 2016).

Importantly, theta oscillations appeared only within certain levels of excitatory afferent activity, around 0.65 Hz for the average firing rate of the Poisson-distributed spike trains. When the 454,700 stimulating afferents in the model (representing the CA3 and entorhinal synapses; calculated in Bezaire and Soltesz (2013)) are active at a Poisson mean of 0.65 Hz, they generate approximately 37,900 incoming spikes / theta cycle, given a theta frequency of 7.8 Hz (Equation 1).

(1) 454,700 afferents0.65 spikes/s7.8 theta cycles/s=37,892 spikes/cycle

Is the latter number of spikes in the afferents to the CA1 network within a physiologically plausible range? The biological CA1 network receives most of its input from CA3 and entorhinal cortical layer III (ECIII), and it has been estimated that about 4% of CA3 pyramidal cells fire up to four spikes per theta wave (Gasparini and Magee, 2006). We previously estimated 204,700 pyramidal cells in ipsilateral CA3 (Bezaire and Soltesz, 2013), giving an estimated 32,750 spikes from ipsilateral CA3 per theta cycle (Equation 2).

(2) 204,700 cells0.04 cell fraction4 spikes/cell=32,752 spikes

About 250,000 principal cells from ipsilateral ECIII synapse onto the CA1 region (Andersen et al., 2006), and approximately 2% of these cells are active per theta cycle at a low firing rate (Csicsvari et al., 1999; Mizuseki et al., 2009). Therefore, ECIII cells could provide 5000 input spikes to ipsilateral CA1 (Equation 3).

(3) 250,000 cells0.02 cell fraction1 spike/cell=5,000 spikes

Therefore, about 37,750 spikes per theta cycle arrive from ipsilateral CA3 and entorhinal cortex to the CA1 network in vivo, which is reassuringly close to the our modeling results indicating that robust theta emerged when the CA1 network model received approximately 37,900 afferent spikes per theta cycle. Thus, the model has the capacity to process a biologically realistic number of spike inputs per cycle while maintaining the theta rhythm.

Our results obtained using the 0.65 Hz excitation indicated that the CA1 model network exhibited phenomena that corresponded well with experimental results, for example, on the differential roles of PV+ basket cells and OLM cells. In addition, the simulations unexpectedly revealed that interneuronal diversity itself may also be important in theta generation, since conversion of all interneurons into fast spiking PV+ basket cells did not result in a network that was conducive for the emergence of theta, in spite of the key role of the PV+ basket cells in hippocampal oscillations. The modeling results also provided the interesting insight that GABAB receptors may play important roles in slow oscillations such as the theta rhythm not because their slow kinetics pace the oscillations, but because their slow kinetics enable a massive charge transfer. This insight was illuminated by the fact that slow GABAB synapses were not necessary for theta as long as their large charge was carried by the fast GABAA synapses. However, we had to increase the conductance of the GABAA synapse almost 300 times to achieve a similar charge transfer as that conveyed by the GABAB synapse. Such a large conductance is not biologically realistic, indicating that the key role for GABAB synapses may be to allow the large synaptic charge transfer via a temporal distribution. Indeed, the importance of GABAB receptors has also been indicated by a number of recent experimental studies, for example, in the modulation of theta and gamma oscillations (Kohl and Paulsen, 2010), setting of spike timing of neuron types during theta (Kohl and Paulsen, 2010), and playing a role in cortical oscillations and memory processes (Craig and McBain, 2014).

In addition to identifying key roles for certain inhibitory components (PV+ interneurons, neurogliaform cells, GABAB, and interneuron diversity), our results also highlighted the importance of the recurrent excitatory collaterals from CA1 pyramidal cells in theta generation in the model of the isolated CA1 network. While it may be expected that isolated theta generation would require local pyramidal cells to provide rhythmic, recurrent excitation to interneurons, our simulations additionally showed that the relatively rare pyramidal cell to pyramidal cell local excitatory connections were also required.

Based on our results, we hypothesize that the inhibitory and excitatory connections within CA1 that were identified to be critical in our perturbation (‘muting’) simulations (Figure 6B) interact to generate the theta waves in the model as follows. Pyramidal cells preferentially discharge at the trough of the LFP analog, strongly recruiting especially the PV+ basket and bistratified cells (green and brown raster plots in Figure 3C), which, in turn, cause a silencing of the pyramidal cells (blue raster plot in Figure 3C) for about the first third of the rising half (i.e., from 0° to about 60°) of the LFP analog theta cycle. As the pyramidal cells begin to emerge from this period of strong inhibition, initially only a few, then progressively more and more pyramidal cells reach firing threshold, culminating in the highest firing probability at the theta trough, completing the cycle. The progressive recruitment of pyramidal cells during the theta cycle appears to be paced according to gamma (see blue raster plot in Figure 3C), and it is likely that the intra-CA1 collaterals of the discharging pyramidal cells play key roles in the step-wise (gamma-paced) recruitment of more and more pyramidal cells as the cycle approaches the following trough. The predicted key roles for physiological pyramidal cell to pyramidal cell connections in theta-gamma generation during running may be tested in future experiments.

Rationale for bases of comparison between modeling results with experimental data

Because our model represented the isolated CA1 network, the modeling results were compared with experimental data from the isolated CA1 preparation when possible. Modeling results for which no corresponding experimental data were available from the isolated CA1 preparation, such as the phase preferential firing of individual interneuron types during theta oscillations, were compared with in vivo data from anesthetized animals (Figure 5B). Experimental results from anesthetized animals offered the most complete data set (e.g., no experimental data were available on CCK basket cells and neurogliaform cells from awake animals, see Figure 5—figure supplement 2). Out of the four interneuronal types for which in vivo data were available from both the awake and anesthetized conditions (Figure 5—figure supplement 2), the phase preference of the axo-axonic cell in the model (163°) was closer to the anesthetized phase (185°) than to the awake phase (251°), whereas the PV+ basket cells in the model displayed phase preferential firing (357°) closer to data reported from awake (289°–310°) than anesthetized animals (234°–271°); the precise reasons underlying these differences are not yet clear. In contrast, bistratified and O-LM cells fired close to the trough in the model, under anesthesia and in awake animals, potentially indicating the primary importance of pyramidal cell inputs in driving these interneurons to fire during theta oscillations under all conditions.

While our model is fundamentally a model of the rat CA1 (e.g., in terms of cell numbers and connectivity; see Table 3 in Bezaire and Soltesz [2013]), some of the electrophysiology data used for constructing the single cell models (Appendix) came from the mouse. In addition, the experimental data on the isolated CA1 preparation were obtained from both rat (Goutagny et al., 2009) and mouse (Amilhon et al., 2015), similar to the experimental results on the phase specific firing in vivo (e.g., awake rat: Lapray et al., 2012; awake mouse: Varga et al., 2014). Because there is no reported evidence for major, systematic differences in key parameters such as the phase specific firing of rat and mouse interneurons in vivo, we did not compare our modeling results with rat and mouse data separately.

A final point concerns the nature of the theta rhythm that emerged in our model. In general, the in vivo theta rhythm has been reported to be either atropine resistant or atropine sensitive, where the former is typically associated with walking and may not be dependent on neuromodulatory inputs, while the latter requires intact, rhythmic cholinergic inputs (Kramis et al., 1975). Given that our model did not explicitly represent neuromodulatory inputs, it is likely that the theta that emerged from our model most closely resembled the atropine resistant form. However, it also plausible that both forms of theta benefit from occurring in a network that is predisposed to oscillate at the theta frequency, as the model network results suggested.

An accessible approach to modeling that balances detail, scale, flexibility and performance

Our results from the strictly data-driven, full-scale CA1 model are consistent with those of earlier models that elegantly demonstrated the basic ingredients capable of producing emergent network oscillations at a range of frequencies in microcircuits and small networks (Rotstein et al., 2005; Siekmeier, 2009; Neymotin et al., 2011b; 2011aFerguson et al., 2013). In addition, our modeling approach also provides a full-scale option to advance the recent studies of network activity propagation and information processing during theta (Cutsuridis et al., 2010; Cutsuridis and Hasselmo, 2012; Taxidis et al., 2013; Saudargiene et al., 2015). Here, we demonstrated that emergent theta and gamma oscillations and theta phase preferential firing are possible even as additional interneuron types are incorporated and the network is scaled up to full size with realistic connectivity including 5 billion synapses between the 300,000-plus cells of our network model.

This work is one step in our broader effort to build a 1:1 model of the entire temporal lobe using a hypothesis-driven model development process, where at each stage of model development the models are used to address specific questions. For example, here we employed our newly developed full-scale CA1 model to gain mechanistic insights into the ability of the intra-CA1 circuitry to generate theta oscillations (Goutagny et al., 2009). The current CA1 network model can be developed into a whole hippocampal or temporal lobe model by replacing the designed CA3 and entorhinal cortical afferents with biophysically detailed CA3, ECIII, and septal networks. While we design our model networks with the motivation to answer a particular question, we keep in mind their potential usage for a broad range of questions. Previously, we built a dentate gyrus model to study epileptic network dynamics (Santhakumar et al., 2005; Morgan and Soltesz, 2008) that was then used by several groups to study disparate topics including epilepsy, network mechanisms of inhibition and excitability, simulation optimization, and modeling software (Migliore et al., 2006; Gleeson et al., 2007; Hines et al., 2008a; 2008b; Hines and Carnevale, 2008; Thomas et al., 2009; Winkels et al., 2009; Cutsuridis et al., 2010; Jedlicka et al., 2010a, 2010b; Thomas et al., 2010; Tejada and Roque, 2014). Our previous model has demonstrated how the resource intensive process of designing a detailed, large-scale model is offset by its potential usage in numerous ways by a multitude of groups. On the other hand, future efforts will be needed to continue to incorporate experimental data obtained by the scientific community on additional, not yet represented parameters into the platform offered by our full-scale CA1 network model, e.g., on cell type-specific gap junctions and short-term plasticity, neuromodulators, diversity of pyramidal cells, glial dynamics, cell to cell variability (e.g., [Schneider et al., 2014]) and others.

We developed a flexible and biologically relevant model that uses computational resources efficiently, positioning the model to be used by the broader community for many future questions. Importantly, the model can be run on the Neuroscience Gateway, an online portal for accessing supercomputers that does not require technical knowledge of supercomputing (https://www.nsgportal.org/). The model is public, well documented, and also well characterized in experimentally relevant terms (See Appendix and online links given in Materials and methods). In addition, all the model configurations and simulation result data sets used in this work are available online (Bezaire et al., 2016b) at (http://doi.org/10.17605/OSF.IO/V4CEH) so the same simulations can easily be repeated with a future, updated model using SimTracker (Bezaire et al., 2016a). Mindful that this model could be used by people with a broad range of modeling experience, we have made freely available our custom software SimTracker (RRID:SCR_014735) that works with the model code to support each step of the modeling process (Bezaire et al., 2016a).

Conclusion and outlook

As highlighted by the BRAIN Initiative, there is an increasing recognition in neurobiology that we must compile our collective experimental observations of the brain into something more cohesive and synergistic than what is being conveyed in individual research articles if we are to fully benefit from the knowledge that we collectively produce (Ramaswamy et al., 2015; Markram et al., 2015). By assimilating our collective knowledge into something as functional as a model, we can further probe the gaps in our experimental studies, setting goals for future experimental work. On the other hand, as powerful new tools are gathering vast quantities of neuroscience data, the extraction and organization of the data itself are becoming a challenge. At least three large programs are undertaking this challenge: the Hippocampome project (for neuroanatomical and electrophysiological data in the hippocampus of mice; [Wheeler et al., 2015]), the Human Brain Project (currently for neuroanatomical and electrophysiological data and models of the rat neocortex, [Ramaswamy et al., 2015]), and NeuroElectro (for electrophysiological data from all species and brain areas; [Tripathy et al., 2014]). These comprehensive databases create the opportunity to build strongly biology-inspired models of entire networks, with all the cells and synapses explicitly represented. Such models are not subject to the connectivity scaling tradeoff wherein smaller networks have unrealistically low levels of input or unrealistically high connectivity between cells. In addition, such models are usable for investigations into an almost infinite number of questions at any level from ion channels, to synapses, to cell types, to microcircuit contributions. This approach represents a new strategy in computational neuroscience, distinct from and complementary to the use of more focused models whose role is to highlight the potential mechanism of a small number of network components.

The scale, flexibility, and accessibility of our strictly data-driven, full-scale CA1 model should aid the modeling of other large scale, detailed, biologically constrained neural networks. The current CA1 network model produces results in agreement with experimental data, but also extends the results to probe the mechanisms of spontaneous theta generation. It provides specific testable predictions that enable focused design of future experiments, as well as providing an accessible resource for the broader community to explore mechanisms of spontaneous theta and gamma generation. Because the model is available at full scale, it is a relevant resource for exploring the transformation of incoming spatial and contextual information to outgoing mnemonic engrams as part of spatial and memory processing, and other pertinent network dynamics.

Materials and methods

All results presented in this work were obtained from simulations of computational models. We implemented our CA1 model in parallel NEURON 7.4, a neural network simulator (Carnevale and Hines, 2005). The model simulations were run with a fixed time step between 0.01 and 0.025 ms, for a simulation duration of 2000 or 4000 ms (except for Figure 6D where one simulation ran for 1600 ms). We executed the simulations on several supercomputers, including Blue Waters at the National Center for Supercomputing Applications at University of Illinois, Stampede and Ranger (retired) at the Texas Advanced Computing Center (TACC), Comet and Trestles at the San Diego Supercomputing Center (SDSC), and the High Performance Computing Cluster at the University of California, Irvine. We used our MATLAB-based SimTracker software tool to design, execute, organize, and analyze the simulations (Bezaire et al., 2016a).

Model development

Request a detailed protocol

The CA1 network model included one type of multicompartmental pyramidal cell with realistic morphology and eight types of interneurons with simplified morphology, including PV+ basket cells, CCK+ basket cells, bistratified cells, axo-axonic cells, O-LM cells, Schaffer Collateral-associated cells, neurogliaform cells, and ivy cells.

Model neurons sometimes behave much differently than expected when subjected to current sweep protocols or synaptic inputs that are outside the range of the original protocols used to construct the model. To ensure the model cells exhibited robust biophysical behavior in a wide range of network conditions, we implemented a standard, thorough characterization strategy for each cell type (Appendix).

The behavior of each cell type was characterized using a current injection sweep that matched experimental conditions reported in the literature. The published experimental data were compared side-by-side with model cell simulation results (Appendix). Model cells were connected via NEURON’s double exponential synapse mechanism (Exp2Syn), with each connection comprising an experimentally observed number of synapses (see Table 1).

The connections between cells were determined with the following algorithm, for each postsynaptic and presynaptic cell type combination:

  1. Calculate the distances between every presynaptic cell and postsynaptic cell of the respective types;

  2. Compute the desired number and distance of connections, as defined by the presynaptic axonal distance distribution and total number of desired connections between these two types; the total number of incoming connections expected by each postsynaptic cell type is divided into radial distance bins and distributed among the bins according to the Gaussian axonal bouton distribution of the presynaptic cell;

  3. Assign each of the possible connections determined in step 2 (connections within the axonal extent of the presynaptic cell) to their respective distance bins, and randomly select a specific number of connections from each bin (the specific number calculated to follow the axonal bouton distribution).

When determining which cells of the model to connect, we distributed all cells evenly within their respective layers in 3D space and enabled random connectivity for cell connections where the postsynaptic cell body fell within the axonal extent of the presynaptic cell (looking in the XY plane only). Each time a connection was established between two cells, the presynaptic cell innervated the experimentally observed number of synapses on the postsynaptic cell. The synapse locations were randomly chosen from all possible places on the cell where the presynaptic cell type had been experimentally observed to innervate. The random number generator used was NEURON’s nrnRan4int.

Biological constraints

Request a detailed protocol

The cell number and connectivity parameters were exactly as we reported previously in our in-depth quantitative assessment of anatomical data about the CA1 (Bezaire and Soltesz, 2013). In the latter paper that formed the data-base for the current full-scale model, we combined immunohistochemical data about laminar distribution and coexpression of markers to estimate the number of each interneuron type in CA1. We then extracted from the experimental literature bouton and input synapse counts for each cell type and multiplied these counts by our estimated number of each cell and determined the available input synapses and boutons in each layer of CA1. The number of connections each cell type was likely to make with every other cell type was based on the results of our quantitative assessment. As the quantitative assessment did not make detailed, interneuron type-specific estimates of connections between interneurons, we performed additional calculations to arrive at the numbers of connections between each type of interneuron in our model. Briefly, we determined the number of inhibitory boutons available for synapsing on interneurons within each layer of CA1. Then, we distributed these connections uniformly across the available incoming inhibitory synapses onto interneurons that we had calculated for that layer. We calculated available incoming synapses by using published experimental observations of inhibitory synapse density on interneuron dendrites by cell class and layer in CA1, which we combined with known anatomical data regarding the dendritic lengths of each interneuron type per layer. We therefore made the following assumption: All available incoming inhibitory synapses onto interneurons in a layer have an equal chance of being innervated by the available inhibitory boutons targeting interneurons in that layer. For further details of the exact calculations, please see the Appendix.

The electrophysiology of each cell was tuned using a combination of manual and optimization techniques. We first fit each cell’s resting membrane potential, capacitance, time constant, and input resistance, followed by hyperpolarized properties such as the sag amplitude and time constant, followed by subthreshold depolarized properties such as a transient peak response, and finally active properties such as spike threshold, rheobase, firing rate, action potential width, height, and afterhyperpolarization. For some cells, we employed the Multiple Run Fitter tool within NEURON to simultaneously fit multiple ion channel conductances. The characterization of each cell type, as well as its comparison to experimental data from the same cell type, is included in the Appendix.

After fitting the cell model properties, we simulated paired recordings to characterize the connections between our model cells. Where experimental data existed for paired recordings, we matched the experimental holding potential and synapse reversal potential, then performed 10 different paired recordings. We characterized the average synapse properties from those 10 runs, including the synaptic amplitude, 10–90% rise time, and decay time constant. Finally, we tuned the synaptic weights and time constants to fit our averages to the experimental data.

To determine the synaptic weights and kinetics for those connections that have not yet been experimentally characterized, we used a novel modeling strategy we call Network Clamp, described in Bezaire et al. (2016a). As experimental paired recording data were not available to directly constrain the synapse properties, we instead constrained the firing rate of the cell in the context of the in vivo network, for which experimental data have been published. We innervated the cell with the connections it was expected to receive in vivo, and then sent artificial spike trains through those connections, ensuring that the properties of the spike trains matched the behavior expected from each cell in vivo during theta (firing rate, level of theta modulation, preferred theta firing phase). Next, we adjusted the weight of the afferent excitatory synapses onto the cell (starting from experimentally observed values for other connections involving that cell type) until the cell achieved a realistic firing rate similar to had been experimentally observed in vivo.

Stimulation

Request a detailed protocol

As none of the model neurons in our model CA1 network are spontaneously active, it was necessary to provide excitatory input to them by stimulating their CA3 and entorhinal cortex synapses. Although the model code is structured to allow the addition of detailed CA3 and cortical inputs, the stimulation patterns used in the present study were not representative of the information content thought to be carried via inputs from those areas, because the focus was on the function of the CA1 network in isolation from rhythmic extra-CA1 influences. In accordance with experimental evidence of spontaneous neurotransmitter release (Kavalali, 2015), we modeled the activation of CA3 and entorhinal synapses as independent Poisson stochastic processes. The model neurons were connected to a subset of these afferents, such that they received a constant level of excitatory synaptic input.

We constrained the synapse numbers and positions of the stimulating afferents using anatomical data. To constrain the afferent synapse weights, we used an iterative process to determine the combinations of synaptic weights that enabled most of the interneurons to fire similar to their observed in vivo firing rates (Figure 5—figure supplement 1 and Table 6). First, we used the output of an initial full-scale simulation to run network clamp simulations on a single interneuron type, altering the incoming afferent synapse weights (but not the incoming spike trains) until the interneuron type fired at a reasonable rate. Then, we applied the synaptic weight to the afferent connections onto that interneuron type in the full-scale model. The resulting simulation then led to a new network dynamic as the constrained activity of that interneuron type caused changes in other interneuron activity. We then performed this exercise for each interneuron type as necessary until we achieved a network where all cell types participated without firing at too high of a level. CCK+ cells had a steep response to the weight of the incoming afferent synapses, remaining silent until the weight was increased significantly and then spiking at a high rate, see Figure 5—figure supplement 1; the particular difficulty in obtaining the in vivo observed firing rate for CCK+ cells in the model may indicate that in vivo they may be strongly regulated by extra-CA1 inhibitory inputs (e.g., from the lateral entorhinal cortex; see Basu et al. (2016) that are not included in the isolated CA1 model).

Analysis of simulation results

Request a detailed protocol

We analyzed the results of each simulation with standard neural data analysis methods provided by our SimTracker software, RRID:SCR_014735, discussed in Bezaire et al. (2016a), including the spike density function (SDF) of all pyramidal neuron spikes (Szucs, 1998), the periodogram of the SDF, and the spectrogram of the LFP analog (Goutagny et al., 2009). We determined the dominant theta and gamma frequencies for the network as the peak in the power spectral density estimate obtained by the spectrogram, and confirmed that those peaks are identical for the SDF and the LFP analog. After finding a dominant theta or gamma frequency, we then analyzed the level of modulation and preferred firing phase for each cell type. Finally, we calculated the firing rate of each cell type.

LFP analog

Request a detailed protocol

We calculated an approximation of the LFP generated by the model neurons based on the method described by Schomburg et al. (2012). For each pyramidal cell within 100 μm of a reference electrode location in stratum pyramidale (coordinates = longitudinal: 200 μm; transverse: 500 μm; height from base of stratum oriens: 120 μm), the contribution to extracellular potential at each point along the dendritic and axonal morphology was recorded using NEURON’s extracellular mechanism and scaled in inverse proportion to the distance from the electrode. In order to reduce the computational load of the simulation, 10% of the pyramidal cells outside the 100 μm radius were randomly selected; their distance-scaled extracellular potentials were scaled up by a factor of 10 and then added to the contributions of the inner cells. We performed reference simulations and LFP analog calculations with the inner radius set to 200 μm and 500 μm and obtained results identical with those in Figures 3 and 4 (where an inner radius of 100 μm was used), except for negligible increases in the theta oscillation power found in the LFP analog spectrogram.

Spike density function

Request a detailed protocol

We calculated the spike density function (SDF) of all pyramidal cell spikes using a Gaussian kernel with a window of 3 ms and a bin size of 1 ms (Szucs, 1998). To see how a cell’s spiking activity is related to its SDF, see Figure 4—figure supplement 1.

Oscillations

Request a detailed protocol

To quantify the frequency and power of the oscillations of the network, we computed a one-sided Welch’s Periodogram of the SDF (Colgin et al., 2009) using a Hamming window with 50% overlap. To characterize the stability of the theta oscillation, we ran the control network for 4 s and then computed the spectrogram of the SDF and of the LFP analog using an analysis script from Goutagny et al. (2009) based on the mtspecgramc function from the Chronux toolbox (http://chronux.org/).

Spike phases and theta modulation

Request a detailed protocol

We calculated the preferred firing theta phases of each cell, using all the spikes of that cell type that occurred after the first 50 ms of the simulation, relative to the filtered LFP analog. The spike times were converted to theta phases, relative to the troughs of the LFP analog theta cycle in which they fired. We then subjected the spike phases to a Rayleigh test to determine the level of theta modulation of the firing of each cell type (Varga et al., 2014).

Firing rates

Request a detailed protocol

The firing rates of the cells were calculated by cropping the first 50 ms of the simulation to remove the initial effects, and then dividing the resulting number of spikes of each cell type by the total number of cells of that type and the duration of the simulation. An alternate average firing rate was calculated by dividing by the number of active cells of that type rather than all of the cells of that type, which gave the average firing rate over all firing cells instead, to better compare with experimentally observed firing rate averages.

Statistical comparison of theta power

Request a detailed protocol

For the GABAB-related simulations, we ran three of each condition and then performed an ANOVA to test for significance in the difference of theta power among the conditions.

Cross correlation of theta and gamma

Request a detailed protocol

To investigate whether a relationship existed between the simultaneous theta and gamma oscillations found in the LFP analog of our control simulation, we filtered the LFP analog signal within the theta range (5–10 Hz) and the gamma range (25–80 Hz). We applied a Hilbert transform to each filtered signal and then compared the phase of the theta-filtered signal with the envelope of the gamma-filtered signal to determine the extent to which theta could modulate the gamma oscillation.

Accessibility

Request a detailed protocol

Our model code is available online at ModelDB, entry #187604 (https://senselab.med.yale.edu/ModelDB/showModel.cshtml?model=187604; code version used to produce results in this work) and Open Source Brain, project nc_ca1 (http://opensourcebrain.org/projects/nc_ca1; most recent code version). Open Source Brain provides tools for users to characterize and inspect model components. The model is also characterized online at http://mariannebezaire.com/models/ca1, along with a graphical explanation of our quantitative assessment used to constrain the model connectivity Bezaire and Soltesz (2013), as well as links to our model code and model results, and detailed instruction manuals for our NEURON code and SimTracker tool, RRID:SCR_014735 (Bezaire et al., 2016a).

For those who wish to view and analyze our simulation results without rerunning the simulation, our simulation results are available on the Open Science Framework (RRID:SCR_003238) at http://doi.org/10.17605/OSF.IO/V4CEH (Bezaire et al., 2016b). Our analyses of these data can be recreated using our publicly available SimTracker tool.

SimTracker is freely available online at http://mariannebezaire.com/simtracker/ and is also listed in SimToolDB, entry #153281 at https://senselab.med.yale.edu/SimToolDB/showTool.cshtml?Tool=153281. The tool is offered both as a stand-alone, compiled version for those without access to MATLAB (for Windows, Mac OS X, and Linux operating systems), and as a collection of MATLAB scripts for those with MATLAB access.

Appendix 1

Experimental cell characterization

These experimental cells were whole cell patch clamped to record the intracellular somatic membrane potential during a range of current injections (‘current sweep’). The firing rates of each cell type are plotted in separate graphs, shown in Appendix 1—figure 1. The electrophysiological properties of each cell type are given in Appendix 1—table 1. In the Appendix, section 'Model cell characterization', the model cells are compared with this experimental data using the same calculations and properties.

Appendix 1—figure 1
Firing Rates of Experimental Cells.

Rebound spiking, which occurs in some O-LM cells at hyperpolarized current injection levels, is not shown in this graph.

https://doi.org/10.7554/eLife.18566.067
Appendix 1—table 1

Intrinsic electrophysiological properties of experimental cells.

https://doi.org/10.7554/eLife.18566.068
Cell typenRMP (mV)Input resistance (MΩ)Sag amplitude (mV)Sag tau (ms)Membrane tau (ms)Rheobase (pA)ISI (ms)Threshold (mV)Spike amplitude (mV)AHP (mV)
From mouse
Pyr17−70.7 ± 1.2139.5 ± 38.87.0 ± 2.234.4 ± 11.021.5 ± 8.6182.4 ± 55.7134.0 ± 44.0−36.7 ± 2.678.2 ± 7.28.6 ± 2.1
Axo3−64.4 ± 4.5122.0 ± 57.51.7 ± 0.645.4 ± 6.911.9 ± 2.2283.3 ± 152.847.8 ± 28.5−31.8 ± 3.444.5 ± 6.716.6 ± 3.5
Bis3−63.6 ± 4.7109.1 ± 30.51.7 ± 0.662.3 ± 13.712.2 ± 0.6333.3 ± 57.724.5 ± 21.8−31.9 ± 4.247.3 ± 6.822.6 ± 0.7
O-LM3−64.8 ± 1.3592.3 ± 97.010.4 ± 3.878.5 ± 22.041.4 ± 11.720.0 ± 0.0101.9 ± 30.1−44.2 ± 2.376.3 ± 6.122.1 ± 4.7
PV+B7−61.4 ± 2.065.2 ± 16.21.8 ± 0.562.9 ± 16.313.3 ± 5.4307.1 ± 109.774.2 ± 36.4−35.3 ± 3.751.1 ± 9.018.0 ± 2.7
From rat
CCK+B1−61.2298.12.772.156.060.0261.0−37.763.715.5
Ivy2−62.3 ± 0.3267.2 ± 107.92.4 ± 2.591.1 ± 120.9171.9 ± 45.680.0 ± 28.374.9 ± 20.6−32.8 ± 0.748.2 ± 5.120.1 ± 2.6
NGF2−66.7 ± 13.4260.0 ± 73.61.8 ± 1.661.7 ± 77.077.2 ± 66.2110.0 ± 70.780.0 ± 28.4−34.0 ± 2.234.7 ± 4.916.2 ± 6.3
SC-A2−57.0 ± 4.3529.9 ± 2.97.9 ± 6.291.1 ± 21.774.2 ± 37.330.0 ± 14.1132.4 ± 29.4−34.3 ± 2.258.7 ± 4.512.6 ± 2.0

The cell references and animals they came from (both rat and mouse (RRID:IMSR_JAX:008069, RRID:IMSR_JAX:007905, RRID:IMSR_JAX:000664), species identified for each experimental cell) are provided here, as well as in two Open Science Framework entries online (O-LM cells: http://dx.doi.org/10.17605/OSF.IO/RA8MW and other cells: http://dx.doi.org/10.17605/OSF.IO/M5EDM) where the raw AxoClamp files of these experiments are also provided (Lee et al., 2016; Quattrocolo and Maccaferri, 2016). The tables of experimental conditions associated with the data entries in Open Science Framework are reproduced here for convenience, in Appendix 1—table 2.

The properties were calculated as follows:

RMP

Resting membrane potential, in units of mV, is calculated as the average membrane potential during a current injection of 0 pA. If 0 was not part of the injection sweep, then the average membrane potential prior to the onset of a different current injection value is used.

Input resistance

Units of MegaOhms (MΩ), is the input resistance calculated from the least hyperpolarized current injection level.

Sag amplitude

Units of mV, computed from the most hyperpolarized current injection level as the difference between the steady state membrane potential towards the end of the current injection and the most hyperpolarized potential achieved towards the beginning of the current injection.

Sag tau

Units of ms, also computed from the most hyperpolarized current injection level, as the time constant required for an equation of the form A*(1-exp(-t/τ))4 to best fit the potential trajectory from the most hyperpolarized point during the current injection until the trace reaches steady state.

Membrane tau

Units of ms, computed from the least hyperpolarized current injection level, as the time constant required for an equation of the form A*(exp(-t/τm)) to best fit the potential trajectory from the onset of the current injection until the trace reaches steady state.

Rheobase

Units of pA, the least depolarized current injection level that resulted in the cell spiking during the current injection (i.e., not as a rebound spike after the injection ends, which can happen for certain cell types after a sufficiently hyperpolarized injection).

ISI

Units of ms, the average time interval between spike threshold time points for the least depolarized current injection level where the cell spiked regularly.

Threshold

Units of mV, the average threshold of the first three spikes for the least depolarized current injection level where the cell spiked regularly. For all experimental and model cells except for the experimental pyramidal cells, the threshold was calculated using CellData’s method #2, where the threshold is the first point where dV/dt exceeds some cutoff value (Cooper et al., 2003; Metz et al., 2005); in our case the cutoff was 28 mV/ms. Because calculating the threshold of the experimental pyramidal cells by this method resulted in a threshold point that was visually too depolarized given the shape of the action potential, CellData’s method #1 was used instead, in which the threshold is the first point where dV/dt>mean(dV/dt)+2*std(dV/dt), meaning the derivative of potential with time exceeds two standard deviations of the average (Atherton and Bevan, 2005).

Spike amplitude

Units of mV, the difference between the membrane potential at the peak of the action potential and the membrane potential at the threshold of the action potential, averaged for the first three spikes of the least depolarized current injection where the cell spiked regularly.

Slow AHP amplitude

Units of mV, also referred to simply as ‘AHP’ in the Appendix, the difference between the membrane potential at the most hyperpolarized potential following the action potential and the membrane potential at the threshold of the action potential, averaged for the first three spikes of the least depolarized current injection where the cell spiked regularly.

Further properties characterized from experimental cells (recorded and published by other labs) are available at NeuroElectro’s website (http://neuroelectro.org), although the data included there are from a wide variety of conditions, animal types, and experimental protocols (and the calculations of properties may have been carried out differently).

Appendix 1—table 2

AxoClamp raw data files. Sch. Coll.-Assoc.: Schaffer Collateral-Associated; Super: superficial. Current sweep injection levels are reported as minimum (most hyperpolarized) : step size : maximum (depolarized) level in units of pA.

https://doi.org/10.7554/eLife.18566.069
Cell typeLabCell nameCurrent inj.Original use and methods reference
SpeciesLevels (pA)
Axo-axonicSolteszCA203LF57mouse−200:50:+500unpublished
Axo-axonicSolteszCA204LF59mouse−200:50:+300unpublished
Axo-axonicSolteszCA204RF59mouse−200:50:+400unpublished
BistratifiedSolteszPV16IMmouse−300:50:+400unpublished
BistratifiedSolteszPV74mouse−300:50:+350unpublished
BistratifiedSolteszPV27IMmouse−300:50:+450unpublished
PV+ BasketSolteszPV34mouse−300:50:+500Lee et al. (2014)
PV+ BasketSolteszPV36mouse−300:50:+800Lee et al. (2014)
PV+ BasketSolteszPV37mouse−300:50:+500Lee et al. (2014)
PV+ BasketSolteszPV38mouse−300:50:+300Lee et al. (2014)
PV+ BasketSolteszPV72mouse−300:50:+400Lee et al. (2014)
PV+ BasketSolteszPV80mouse−300:50:+450Lee et al. (2014)
Deep PyramidalSolteszD1_25abfmouse−400:50:+550Lee et al. (2014)
Deep PyramidalSolteszD1_45abfmouse−400:50:+550Lee et al. (2014)
Deep PyramidalSolteszD2_06abfmouse−400:50:+550Lee et al. (2014)
Deep PyramidalSolteszD2_49abfmouse−400:50:+550Lee et al. (2014)
Deep PyramidalSolteszD3_55abfmouse−400:50:+550Lee et al. (2014)
Deep PyramidalSolteszD4_11abfmouse−400:50:+550Lee et al. (2014)
Deep PyramidalSolteszD5_15abfmouse−400:50:+550Lee et al. (2014)
Deep PyramidalSolteszD6_19abfmouse−400:50:+550Lee et al. (2014)
Deep PyramidalSolteszD7mouse−400:50:+550Lee et al. (2014)
Super. PyramidalSolteszS1_04abfmouse−400:50:+550Lee et al. (2014)
Super. PyramidalSolteszS1_47abfmouse−400:50:+550Lee et al. (2014)
Super. PyramidalSolteszS2_08abfmouse−400:50:+550Lee et al. (2014)
Super. PyramidalSolteszS2_31abfmouse−400:50:+550Lee et al. (2014)
Super. PyramidalSolteszS2_51abfmouse−400:50:+550Lee et al. (2014)
Super. PyramidalSolteszS3_13abfmouse−400:50:+550Lee et al. (2014)
Super. PyramidalSolteszS4mouse−400:50:+550Lee et al. (2014)
Super. PyramidalSolteszS5_21abfmouse−400:50:+550Lee et al. (2014)
IvySoltesz0422–1 (File 5)rat−100:20:+890Krook-Magnuson et al. (2011)
IvySoltesz0428–1 (File 4)rat−100:20:+300Krook-Magnuson et al. (2011)
NeurogliaformSoltesz09o21 (File 4)rat−100:20:+120Krook-Magnuson et al. (2011)
NeurogliaformSoltesz09o27 (File 7)rat−100:20:+490Krook-Magnuson et al. (2011)
CCK+ BasketSolteszsh108_BCrat−100:20:+80Lee et al. (2010)
Sch. Coll.-Assoc.Solteszsh114_SCArat−100:20:+60Lee et al. (2010)
Sch. Coll.-Assoc.Solteszsh153_SCArat−100:20:+60Lee et al. (2010)
O-LMMaccaferri1May2012_P3mouse−100:30:+250Quattrocolo and Maccaferri (2013)
O-LMMaccaferri20Sept2011_P2mouse−100:30:+250Quattrocolo and Maccaferri (2013)
O-LMMaccaferri24October2012_C2mouse−100:30:+250Quattrocolo and Maccaferri (2013)

Model cell characterization

Model cell numbers and structural connectivity are based on Bezaire and Soltesz (2013).

In terms of electrophysiology, each model cell is characterized in experimental terms and compared to the experimental data presented above. A graphical summary of electrophysiological comparison is shown below in Appendix 1—figures 24, and further details of intrinsic physiology and synaptic characterization follows. The detailed information is presented as a single figure spanning two pages per cell, where the subfigure panels may contain figures or tables, to better group and arrange the data. For each cell type, the same information is provided:

A: Model current sweep

The somatic intracellular membrane potential recording for the model cell is shown, in response to all hyperpolarized and the most depolarized current injection

B: Experimental current sweep

The somatic intracellular membrane potential recording for an experimental cell (one of the ones featured in the previous section) is shown, in response to all hyperpolarized and the depolarized current injection level closest to the one shown for the model cell.

C: Model electrophysiological property table

The values of each electrophysiological property are shown for the model cell, measured or calculated in the same way as for the experimental cells in the previous section. All experimental data are taken from the cells used in section 'Experimental cell characterization' of the Appendix. More experimental data, obtained under a wider variety of conditions, are available at http://neuroelectro.org.

D: Firing rates

The firing rate of the model cell as a function of current injection is shown for a range of currents. The firing rates of the experimental cells are shown as well (to identify specific experimental cells, see the firing rate graphs in section 'Experimental cell characterization' of the Appendix).

E: Ion channel table

Each ion channel type present in the model cell type is listed here, along with the maximum conductance density of that channel (which may occur in the soma or in another part of the cell). Further details about the ion channels are available in the Appendix, sections 'Ion channel descriptions' and 'Ion channel equations'.

F: Structural connectivity table

The structural basis of the model connections is provided here in terms of both connections (comprising multiple synapses) and synapses, shown for convergence onto the cell (left side) and divergence emanating from the cell (right side). Connectivity involving pyramidal cells (either pre- or postsynaptically) is based on Bezaire and Soltesz (2013) while connectivity between interneurons is detailed in section 'Inhibitory connectivity' in the Appendix. Blank or missing rows for specific cell types indicate that there were no connections with that cell type.

G,H: Model paired recordings – experimental conditions table

For incoming (G) and outgoing (H) connections involving a given cell type, any experimental constraints used to fit the model connections are cited here, including the conditions of the experiment that the model reproduced (the holding potential of the cell and the reversal potential of the synapse as set by the bath and pipette solutions used in the experiment), along with the model connection characterization (amplitude, 10–90% rise time, and decay time constant) and the percentage difference from the mean experimental properties. The model results were reported as the average from 10 random connections between the two cell types, and wherever possible were compared to an actual experimental mean (leaving out failed responses) rather than an effective experimental mean that factored in the synapse failure rate. For connections that were characterized in current clamp rather than voltage clamp, the data were colored in purple rather than black to indicate their different properties (half-width instead of decay time constant, for example). Because the experimental data did not routinely report or, in a standard way, calculate the junction potential of the experiment, we did not factor in the difference in holding potential (and hence driving force of the synapse) due to junction potential; this resulted in a slight over or underestimation of the driving force and hence the actual synaptic conductance, depending on the junction potential value.

I: Model synapse parameters table

The model parameters used in the synapses (resulting from the experimental tuning above, or firing rate tuning if no experimental connection data were available) are listed here, with the parameters for connections onto that cell type from the other types listed on the left side, and the parameters for connections onto other cell types, from that cell type, listed on the right side.

J: Model physiological connections table

All connections in the model were recharacterized under the same condition, as experimental data were gathered under diverse conditions where different connections could not be directly compared. In this table, all connections were recorded using the physiological reversal potential while holding the postsynaptic cell at −50 mV voltage clamp (and not accounting for any junction potential as would occur if an experiment were to replicate these model simulation conditions). The connections from other cell types to the given cell type are shown on the left side; connections from the given cell type to other cell types are shown on the right side. Blank lines or missing rows indicate there is not a connection between those cell types in the model.

K, L: Model physiological connections

The postsynaptic current (PSC) responses for all the physiological connections as detailed in Table J above are graphed here, with connections from other cell types to the given cell type shown in panel K, and connections from the given cell type to other cells shown in panel L.

The model cell types can be further characterized as desired by using the artificially generated AxoClamp files included in the Source Data for Figure 2. These AxoClamp files were generated from our CellClamp tool within SimTracker, providing the same data format (but with generic header lines for the first ten lines), as the tab-delimited ATF file format.

Model and experimental electrophysiology
Appendix 1—figure 2
Physiological properties of experimental and model cells.

Experimental data are shown with closed markers for the mean and error bars for cell types where n > 1. The model cell properties are plotted as open circles. Calculation of properties is explained in the text. (A) resting membrane potential, (B) threshold, and (C) spike amplitude.

https://doi.org/10.7554/eLife.18566.070
Appendix 1—figure 3
Physiological properties, continued.

(A) sag time constant, (B) sag amplitude, and (C) amplitude of afterhyperpolarization (AHP).

https://doi.org/10.7554/eLife.18566.071
Appendix 1—figure 4
Physiological properties, continued.

(A) rheobase, (B) membrane time constant, (C) interspike interval (ISI), and (D) input resistance.

https://doi.org/10.7554/eLife.18566.072

Pyramidal cell: principal cell (311500 Cells)

Appendix 1—figure 5
Pyramidal (A) model and (B) experimental current sweep.

(C) Firing rates of model and experimental cells.

https://doi.org/10.7554/eLife.18566.073
Appendix 1—table 3

Model Pyramidal cell electrophysiological properties.

https://doi.org/10.7554/eLife.18566.074
PropertyValue
RMP−63.0 mV
Input Resistance76.1 MΩ
Sag Amplitude6.5 mV
Sag Tau9.6 ms
Membrane Tau7.1 ms
Rheobase250.0 pA
ISI80.7 ms
Threshold−39.9 mV
Spike Amplitude80.3 mV
Slow AHP Amplitude14.3 mV
Appendix 1—table 4

Model Pyramidal cell ion channels and conductance at highest density location in cell.

https://doi.org/10.7554/eLife.18566.075
ChannelHighest conductance Gmax (S/cm2)
HCNp4.968e-03
Kdrp3.000e-03
KvAdistp4.682e-02
KvAproxp1.599e-02
Navaxonp6.400e-02
Navp3.200e-02
Model and experimental connectivity
Appendix 1—table 5

Structural connection parameters for Pyramidal cells, based on Bezaire and Soltesz (2013).

https://doi.org/10.7554/eLife.18566.076
Other typeOther cell to pyrPyr to other cell
#Syn.s#Post#Syn.s#Post
Conn.s/Conn.#Loc.Conn.s/Conn.#Loc.
Axo6636axon132apical dendrite
Bis1010100any dendrite337apical dendrite
CCK+B138104any dendrite
Ivy4210420any dendrite030apical dendrite
NGF1410140apical dendrite
O-LM81080apical dendrite13337basal dendrite
Pyr1971197apical dendrite1971197apical dendrite
PV+B1711187soma8322apical dendrite
SC-A030apical dendrite
CA35985211970any dendrite
ECIII129922598any dendrite
Experimental connection constraints
Appendix 1—table 6

Experimental constraints for incoming connections onto Pyramidal cells (clamp: black=voltage; purple=current).

https://doi.org/10.7554/eLife.18566.077
Pre typeExp. ref.

Hold (mV)

Erev (mV)Amp. (pA,mV)Diff. %t10-90 (ms)Diff. %τdecay (ms)Diff. %
AxoMaccaferri et al., 2000−70.07.0323.78+5.10.83+3.111.20+0.0
BisMaccaferri et al., 2000−70.07.0143.21−4.52.22+11.215.40−4.3
CCK+BLee et al., 2010−70.0−26.0118.97+3.10.53−16.76.15−4.9
IvyFuentealba et al., 2008−50.0−88.08.17+2.13.50+25.015.43−3.9
NGFPrice et al., 2008−50.0−89.05.25+7.115.48−3.932.73−34.5
O-LMMaccaferri et al., 2000−70.07.024.35−6.34.68−24.618.88−9.3
PyrDeuchars and Thomson, 1996−67.00.00.60−14.56.00+122.220.55+22.3
PV+BSzabadics et al., 2007−70.0−26.091.94−13.90.50−5.76.70+4.7
SC-ALee et al., 2010−70.0−26.052.42−12.91.63+13.68.55+3.0
Appendix 1—table 7

Experimental constraints for outgoing connections from Pyramidal cells (clamp: black=voltage; purple=current).

https://doi.org/10.7554/eLife.18566.078
Post typeExp. ref.Hold (mV)Erev (mV)Amp. (pA,mV)Diff. %t10-90 (ms)Diff. %τdecay (ms)Diff. %
BisPawelzik et al., 2002−66.00.00.77−19.61.58+31.316.75+41.1
IvyFuentealba et al., 2008−65.8−70.00.06−97.91.38−8.321.35+41.1
PyrDeuchars and Thomson, 1996−67.00.00.60−14.56.00+122.219.05+22.3
PV+BLee et al., 2014−60.00.015.09−67.70.28−72.52.00+22.3
Model synapse parameters
Appendix 1—table 8

Model synaptic parameters for Pyramidal cells in the control network.

https://doi.org/10.7554/eLife.18566.079
TypeOther cell to pyrPyr to other cell
Erev (mV)Gmax (nS)τrise (ms)τdecay (ms)Erev (mV)Gmax (nS)τrise (ms)τdecay (ms)
Axo−60.01.150e-030.288.400.04.000e-050.300.60
Bis−60.05.100e-040.119.700.01.900e-030.110.25
CCK+B−60.05.200e-040.204.20
Ivy−60.04.100e-051.1011.000.04.050e-040.300.60
NGF−60.06.500e-059.0039.00
O-LM−60.03.000e-040.1311.000.02.000e-040.300.60
Pyr0.07.000e-020.101.500.07.000e-020.101.50
PV+B−60.02.000e-040.306.200.07.000e-040.070.20
SC-A0.04.050e-040.300.60
CA30.02.000e-040.503.00
ECIII0.02.000e-040.503.00
Physiological characterization of model connections
Appendix 1—table 9

Model synaptic properties under voltage clamp at −50 mV with physiological reversal potentials

https://doi.org/10.7554/eLife.18566.080
TypeOther cell to pyrPyr to other cell
Hold (mV)Erev (mV)Amp. (pA)t10-90 (ms)τdecay (ms)Hold (mV)Erev (mV)Amp. (pA)t10-90 (ms)τdecay (ms)
Axo−50.0−60.036.450.8511.57−50.00.01.850.782.53
Bis−50.0−60.013.472.1715.20−50.00.064.480.281.42
CCK+B−50.0−60.024.860.526.03
Ivy−50.0−60.01.633.6315.35−50.00.040.700.581.28
NGF−50.0−60.01.1065.580.00
O-LM−50.0−60.00.543.7014.10−50.00.017.470.601.53
Pyr−50.00.022.132.229.65−50.00.022.132.229.65
PV+B−50.0−60.020.560.506.70−50.00.014.750.251.77
SC-A−50.00.017.420.683.05
CA3−50.00.07.151.837.08
ECIII−50.00.01.413.2513.63
Appendix 1—figure 6
Connections onto (A) and (B) from model Pyramidal cells, under voltage clamp at −50 mV with physiological reversal potentials.
https://doi.org/10.7554/eLife.18566.081

Axo-axonic cell: fast-spiking axonic inhibitor (1470 Cells)

Model and experimental electrophysiology
Appendix 1—figure 7
Axo-axonic (A) model and (B) experimental current sweep.

(C) Firing rates of model and experimental cells.

https://doi.org/10.7554/eLife.18566.082
Appendix 1—table 10

Model Axo-axonic cell electrophysiological properties.

https://doi.org/10.7554/eLife.18566.083
PropertyValue
RMP−65.0 mV
Input Resistance52.3 MΩ
Sag Amplitude
Sag Tau
Membrane Tau7.0 ms
Rheobase200.0 pA
ISI57.3 ms
Threshold−42.0 mV
Spike Amplitude94.3 mV
Slow AHP Amplitude33.4 mV
Appendix 1—table 11

Model Axo-axonic cell ion channels and conductance at highest density location in cell.

https://doi.org/10.7554/eLife.18566.084
ChannelGmax (S/cm2)
CavL5.000e-03
CavN8.000e-04
KCaS2.000e-06
Kdrfast1.300e-02
KvA1.500e-04
KvCaB2.000e-07
Nav1.500e-01
leak1.800e-04
Model and experimental connectivity
Appendix 1—table 12

Structural connection parameters for Axo-axonic cells, based on Bezaire and Soltesz (2013).

https://doi.org/10.7554/eLife.18566.085
Other typeOther cell to axoAxo to other cell
#Syn.s#Post#Syn.s#Post
Conn.s/Conn.#Loc.Conn.s/Conn.#Loc.
Bis1610160any dendrite
CCK+B12896any dendrite
Ivy2410240any dendrite
O-LM81080apical dendrite
Pyr1623486apical dendrite127167628axon
PV+B39139soma
SC-A166any dendrite
CA3417028340any dendrite
ECIII4852970any dendrite
Experimental connection constraints

Note:No experimental constraints available for incoming synapses to Axo-axonic cells.

Appendix 1—table 13

Experimental constraints for outgoing connections from Axo-axonic cells (clamp: black=voltage; purple=current).

https://doi.org/10.7554/eLife.18566.086
Post typeExp. ref.Hold (mV)Erev (mV)Amp. (pA,mV)Diff. %t10-90 (ms)Diff. %τdecay (ms)Diff. %
PyrMaccaferri et al., 2000−70.07.0323.78+5.10.83+3.111.20+0.0
Model synapse parameters
Appendix 1—table 14

Model synaptic parameters for Axo-axonic cells in the control network.

https://doi.org/10.7554/eLife.18566.087
TypeOther cell to axoAxo to other cell
Erev (mV)Gmax (nS)τrise (ms)τdecay (ms)Erev (mV)Gmax (nS)τrise (ms)τdecay (ms)
Bis−60.06.000e-040.292.67
CCK+B−60.07.000e-040.434.49
Ivy−60.05.700e-052.903.10
O-LM−60.01.200e-040.7310.00
Pyr0.04.000e-050.300.60−60.01.150e-030.288.40
PV+B−60.01.200e-040.292.67
SC-A−60.06.000e-040.424.99
CA30.01.200e-042.006.30
ECIII0.01.200e-042.006.30
Physiological characterization of model connections
Appendix 1—table 15

Model synaptic properties under voltage clamp at −50 mV with physiological reversal potentials

https://doi.org/10.7554/eLife.18566.088
TypeOther cell to axoAxo to other cell
Hold (mV)Erev (mV)Amp. (pA)t10-90 (ms)τdecay (ms)Hold (mV)Erev (mV)Amp. (pA)t10-90 (ms)τdecay (ms)
Bis−50.0−60.036.770.703.70
CCK+B−50.0−60.047.290.755.27
Ivy−50.0−60.04.342.136.57
O-LM−50.0−60.04.762.5512.03
Pyr−50.00.01.850.782.53−50.0−60.036.450.8511.57
PV+B−50.0−60.01.080.453.13
SC-A−50.0−60.024.001.006.13
CA3−50.00.010.852.308.80
ECIII−50.00.08.743.089.20
Appendix 1—figure 8
Connections onto (A) and (B) from model Axo-axonic cells, under voltage clamp at −50 mV with physiological reversal potentials.
https://doi.org/10.7554/eLife.18566.089

Bistratified cell: fast-spiking dendritic inhibitor (2210 Cells)

Model and Experimental Electrophysiology
Appendix 1—figure 9
Bistratified (A) model and (B) experimental current sweep.

(C) Firing rates of model and experimental cells.

https://doi.org/10.7554/eLife.18566.090
Appendix 1—table 16

Model Bistratified cell electrophysiological properties.

https://doi.org/10.7554/eLife.18566.091
PropertyValue
RMP−67.0 mV
Input Resistance98.8 MΩ
Sag Amplitude0.0 mV
Sag Tau
Membrane Tau14.7 ms
Rheobase350.0 pA
ISI39.0 ms
Threshold−28.1 mV
Spike Amplitude51.2 mV
Slow AHP Amplitude48.8 mV
Appendix 1—table 17

Model Bistratified cell ion channels and conductance at highest density location in cell.

https://doi.org/10.7554/eLife.18566.092
ChannelGmax (S/cm2)
CavL4.000e-03
CavN4.000e-04
KCaS7.000e-07
Kdrfast1.600e-02
KvA5.000e-05
KvCaB7.000e-08
Navbis7.000e-02
leak9.001e-05
Model and experimental connectivity
Appendix 1—table 18

Structural connection parameters for Bistratified cells, based on Bezaire and Soltesz (2013).

https://doi.org/10.7554/eLife.18566.093
Other typeOther cell to bisBis to other cell
#Syn.s#Post#Syn.s#Post
Conn.s/Conn.#Loc.Conn.s/Conn.#Loc.
Axo1110106any dendrite
Bis1610160any dendrite1610160any dendrite
CCK+B12896any dendrite2610260any dendrite
Ivy2410240any dendrite1210119any dendrite
O-LM81080apical dendrite2910289any dendrite
Pyr36631098apical dendrite14101014095any dendrite
PV+B39139soma4010400any dendrite
SC-A166any dendrite31030any dendrite
CA35782211564any dendrite
ECIII4322864any dendrite
Experimental connection constraints
Appendix 1—table 19

Experimental constraints for incoming connections onto Bistratified cells (clamp: black=voltage; purple=current).

https://doi.org/10.7554/eLife.18566.094
Pre typeExp. ref.Hold (mV)Erev (mV)Amp. (pA,mV)Diff. %t10-90 (ms)Diff. %τdecay (ms)Diff. %
PyrPawelzik et al., 2002−66.00.00.77−19.61.58+31.314.68+41.1
PV+BCobb et al., 1997−55.0−70.00.27−27.50.47−52.57.30+30.4
Appendix 1—table 20

Experimental constraints for outgoing connections from Bistratified cells (clamp: black=voltage; purple=current).

https://doi.org/10.7554/eLife.18566.095
Post typeExp. ref.Hold (mV)Erev (mV)Amp. (pA,mV)Diff. %t10-90 (ms)Diff. %τdecay (ms)Diff. %
PyrMaccaferri et al., 2000−70.07.0143.21−4.52.22+11.215.40−4.3
Model synapse parameters
Appendix 1—table 21

Model synaptic parameters for Bistratified cells in the control network.

https://doi.org/10.7554/eLife.18566.096
TypeOther cell to bisBis to other cell
Erev (mV)Gmax (nS)τrise (ms)τdecay (ms)Erev (mV)Gmax (nS)τrise (ms)τdecay (ms)
Axo−60.06.000e-040.292.67
Bis−60.05.100e-040.292.67−60.05.100e-040.292.67
CCK+B−60.07.000e-040.434.49−60.08.000e-040.292.67
Ivy−60.07.700e-052.903.10−60.05.000e-040.292.67
O-LM−60.01.100e-040.6015.00−60.02.000e-051.008.00
Pyr0.01.900e-030.110.25−60.05.100e-040.119.70
PV+B−60.02.900e-030.180.45−60.09.000e-030.292.67
SC-A−60.06.000e-040.424.99−60.08.000e-040.292.67
CA30.01.500e-042.006.30
ECIII0.01.500e-042.006.30
Physiological characterization of model connections
Appendix 1—table 22

Model synaptic properties under voltage clamp at −50 mV with physiological reversal potentials

https://doi.org/10.7554/eLife.18566.097
TypeOther cell to bisBis to other cell
Hold (mV)Erev (mV)Amp. (pA)t10-90 (ms)τdecay (ms)Hold (mV)Erev (mV)Amp. (pA)t10-90 (ms)τdecay (ms)
Axo−50.0−60.036.770.703.70
Bis−50.0−60.034.340.703.72−50.0−60.034.340.703.72
CCK+B−50.0−60.048.130.785.35−50.0−60.048.550.734.15
Ivy−50.0−60.06.392.156.63−50.0−60.043.400.603.17
O-LM−50.0−60.06.312.7017.05−50.0−60.01.861.788.13
Pyr−50.00.064.480.281.42−50.0−60.013.472.1715.20
PV+B−50.0−60.024.450.170.73−50.0−60.0429.340.574.13
SC-A−50.0−60.026.431.026.20−50.0−60.050.350.704.10
CA3−50.00.013.812.388.82
ECIII−50.00.012.043.059.30
Appendix 1—figure 10
Connections onto (A) and (B) from model Bistratified cells, under voltage clamp at −50 mV with physiological reversal potentials.
https://doi.org/10.7554/eLife.18566.098

CCK+ basket cell: regular-spiking somatic inhibitor (3600 cells)

Model and experimental electrophysiology
Appendix 1—figure 11
CCK+ Basket (A) model and (B) experimental current sweep.

(C) Firing rates of model and experimental cells.

https://doi.org/10.7554/eLife.18566.099
Appendix 1—table 23

Model CCK+ Basket cell electrophysiological properties.

https://doi.org/10.7554/eLife.18566.100
PropertyValue
RMP−70.6 mV
Input Resistance222.4 MΩ
Sag Amplitude9.2 mV
Sag Tau45.6 ms
Membrane Tau25.5 ms
Rheobase80.0 pA
ISI180.8 ms
Threshold−38.0 mV
Spike Amplitude65.9 mV
Slow AHP Amplitude32.1 mV
Appendix 1—table 24

Model CCK+ Basket cell ion channels and conductance at highest density location in cell.

https://doi.org/10.7554/eLife.18566.101
ChannelGmax (S/cm2)
CavL2.700e-03
CavN2.000e-05
HCN1.000e-04
KCaS4.000e-06
Kdrfast8.000e-05
KvA4.000e-04
KvCaB4.000e-05
KvGroup2.600e-03
Navcck1.800e-02
leak3.704e-05
Model and experimental connectivity
Appendix 1—table 25

Structural connection parameters for CCK+ Basket cells, based on Bezaire and Soltesz (2013).

https://doi.org/10.7554/eLife.18566.102
Other typeOther cell to CCK+BCCK+B to other cell
#Syn.s#Post#Syn.s#Post
Conn.s/Conn.#Loc.Conn.s/Conn.#Loc.
Axo5839any dendrite
Bis1610160any dendrite7858any dendrite
CCK+B358280any dendrite358280any dendrite
Ivy9610960any dendrite208156any dendrite
O-LM4010400apical dendrite9872any dendrite
Pyr112588998any dendrite
PV+B38138soma188147any dendrite
SC-A6636any dendrite3824any dendrite
CA3200024000any dendrite
ECIII55921118any dendrite
Experimental connection constraints

Note:No experimental constraints available for incoming synapses to CCK+ Basket cells.

Appendix 1—table 26

Experimental constraints for outgoing connections from CCK+ Basket cells (clamp: black=voltage; purple=current).

https://doi.org/10.7554/eLife.18566.103
Post typeExp. ref.Hold (mV)Erev (mV)Amp. (pA,mV)Diff. %t10-90 (ms)Diff. %τdecay (ms)Diff. %
PyrLee et al., 2010−70.0−26.0118.97+3.10.53−16.76.15−4.9
Model synapse parameters
Appendix 1—table 27

Model synaptic parameters for CCK+ Basket cells in the control network.

https://doi.org/10.7554/eLife.18566.104
TypeOther cell to CCK+BCCK+B to other cell
Erev (mV)Gmax (nS)τrise (ms)τdecay (ms)Erev (mV)Gmax (nS)τrise (ms)τdecay (ms)
Axo−60.07.000e-040.434.49
Bis−60.08.000e-040.292.67−60.07.000e-040.434.49
CCK+B−60.04.500e-040.434.49−60.04.500e-040.434.49
Ivy−60.03.700e-052.903.10−60.03.000e-040.434.49
O-LM−60.01.200e-030.7320.20−60.07.000e-041.008.00
Pyr−60.05.200e-040.204.20
PV+B−60.01.200e-030.292.67−60.09.000e-030.434.49
SC-A−60.08.500e-040.424.99−60.07.000e-040.434.49
CA30.06.500e-042.006.30
ECIII0.06.500e-042.006.30
Physiological characterization of model connections
Appendix 1—table 28

Model synaptic properties under voltage clamp at −50 mV with physiological reversal potentials

https://doi.org/10.7554/eLife.18566.105
TypeOther cell to CCK+BCCK+B to other cell
Hold (mV)Erev (mV)Amp. (pA)t10-90 (ms)τdecay (ms)Hold (mV)Erev (mV)Amp. (pA)t10-90 (ms)τdecay (ms)
Axo−50.0−60.047.290.755.27
Bis−50.0−60.048.550.734.15−50.0−60.048.130.785.35
CCK+B−50.0−60.032.190.735.30−50.0−60.032.190.735.30
Ivy−50.0−60.03.002.256.95−50.0−60.022.340.805.05
O-LM−50.0−60.040.323.1028.42−50.0−60.054.981.359.05
Pyr−50.0−60.024.860.526.03
PV+B−50.0−60.011.310.423.08−50.0−60.0523.110.685.70
SC-A−50.0−60.033.811.056.90−50.0−60.049.550.705.38
CA3−50.00.055.242.539.35
ECIII−50.00.043.273.4010.87
Appendix 1—figure 12
Connections onto (A) and (B) from model CCK+ Basket cells, under voltage clamp at −50 mV with physiological reversal potentials.
https://doi.org/10.7554/eLife.18566.106

Ivy cell: late-spiking cell (8810 cells)

Model and experimental electrophysiology
Appendix 1—figure 13
Ivy (A) model and (B) experimental current sweep.

(fig:ivypage:firing) Firing rates of model and experimental cells.

https://doi.org/10.7554/eLife.18566.107
Appendix 1—table 29

Model Ivy cell electrophysiological properties.

https://doi.org/10.7554/eLife.18566.108
PropertyValue
RMP−60.0 mV
Input Resistance100.7 MΩ
Sag Amplitude0.0 mV
Sag Tau
Membrane Tau21.3 ms
Rheobase160.0 pA
ISI305.5 ms
Threshold−27.7 mV
Spike Amplitude54.6 mV
Slow AHP Amplitude20.9 mV
Appendix 1—table 30

Model Ivy cell ion channels and conductance at highest density location in cell.

https://doi.org/10.7554/eLife.18566.109
ChannelGmax (S/cm2)
CavL5.611e-02
CavN5.817e-04
KCaS4.515e-07
Kdrfastngf1.551e-01
KvAngf5.220e-06
KvCaB1.024e-06
Navngf3.786e+00
leak8.471e-05
Model and experimental connectivity
Appendix 1—table 31

Structural connection parameters for Ivy cells, based on Bezaire and Soltesz (2013).

https://doi.org/10.7554/eLife.18566.110
OtherTypeOther cell to ivyIvy to other cell
#Syn.s#Post#Syn.s#Post
Conn.s/Conn.#Loc.Conn.s/Conn.#Loc.
Axo41040any dendrite
Bis31030any dendrite61060any dendrite
CCK+B8864any dendrite3910392any dendrite
Ivy2410240any dendrite2410240any dendrite
NGF1110113any dendrite
O-LM2510253any dendrite
Pyr9327apical dendrite14851014850any dendrite
PV+B818soma1510150any dendrite
SC-A2612any dendrite51046any dendrite
CA3192323846any dendrite
Experimental connection constraints
Appendix 1—table 32

Experimental constraints for incoming connections onto Ivy cells (clamp: black=voltage; purple=current).

https://doi.org/10.7554/eLife.18566.111
Pre typeExp. ref.Hold (mV)Erev (mV)Amp. (pA,mV)Diff. %t10-90 (ms)Diff. %τdecay (ms)Diff. %
PyrFuentealba et al., 2008−65.8−70.00.06−97.91.38−8.3−4.9
Appendix 1—table 33

Experimental constraints for outgoing connections from Ivy cells (clamp: black=voltage; purple=current).

https://doi.org/10.7554/eLife.18566.112
Post typeExp. ref.Hold (mV)Erev (mV)Amp. (pA,mV)Diff. %t10-90 (ms)Diff. %τdecay (ms)Diff. %
PyrFuentealba et al., 2008−50.0−88.08.17+2.13.50+25.015.43−3.9
Model Synapse Parameters
Appendix 1—table 34

Model synaptic parameters for Ivy cells in the control network.

https://doi.org/10.7554/eLife.18566.113
TypeOther cell to ivyIvy to other cell
Erev (mV)Gmax (nS)τrise (ms)τdecay (ms)Erev (mV)Gmax (nS)τrise (ms)τdecay (ms)
Axo−60.05.700e-052.903.10
Bis−60.05.000e-040.292.67−60.07.700e-052.903.10
CCK+B−60.03.000e-040.434.49−60.03.700e-052.903.10
Ivy−60.05.700e-052.903.10−60.05.700e-052.903.10
NGF−60.05.700e-052.903.10
O-LM−60.05.700e-052.903.10
Pyr0.04.050e-040.300.60−60.04.100e-051.1011.00
PV+B−60.01.600e-040.292.67−60.07.000e-042.903.10
SC-A−60.08.500e-040.424.99−60.03.700e-052.903.10
CA30.03.000e-042.006.30
Physiological characterization of model connections
Appendix 1—table 35

Model synaptic properties under voltage clamp at −50 mV with physiological reversal potentials

https://doi.org/10.7554/eLife.18566.114
TypeOther cell to ivyIvy to other cell
Hold (mV)Erev (mV)Amp. (pA)t10-90 (ms)τdecay (ms)Hold (mV)Erev (mV)Amp. (pA)t10-90 (ms)τdecay (ms)
Axo−50.0−60.04.342.136.57
Bis−50.0−60.043.400.603.17−50.0−60.06.392.156.63
CCK+B−50.0−60.022.340.805.05−50.0−60.03.002.256.95
Ivy−50.0−60.05.481.886.42−50.0−60.05.481.886.42
NGF−50.0−60.05.481.886.42
O-LM−50.0−60.05.322.106.33
Pyr−50.00.040.700.581.28−50.0−60.01.633.6315.35
PV+B−50.0−60.01.440.553.13−50.0−60.051.352.056.75
SC-A−50.0−60.046.620.855.58−50.0−60.03.092.226.88
CA3−50.00.029.422.058.60
Appendix 1—figure 14
Connections onto (A) and (B) from model Ivy cells, under voltage clamp at −50 mV with physiological reversal potentials.
https://doi.org/10.7554/eLife.18566.115

Neurogliaform cell: late-spiking feed forward cell (3580 cells)

Model and experimental electrophysiology
Appendix 1—figure 15
Neurogliaform (A) model and (B) experimental current sweep.

(C) Firing rates of model and experimental cells.

https://doi.org/10.7554/eLife.18566.116
Appendix 1—table 36

Model Neurogliaform cell electrophysiological properties.

https://doi.org/10.7554/eLife.18566.117
PropertyValue
RMP−60.0 mV
Input Resistance100.8 MΩ
Sag Amplitude0.0 mV
Sag Tau
Membrane Tau21.3 ms
Rheobase170.0 pA
ISI170.3 ms
Threshold−27.8 mV
Spike Amplitude55.2 mV
Slow AHP Amplitude20.6 mV
Appendix 1—table 37

Model Neurogliaform cell ion channels and conductance at highest density location in cell.

https://doi.org/10.7554/eLife.18566.118
ChannelGmax (S/cm2)
CavL5.611e-02
CavN5.817e-04
KCaS4.515e-07
Kdrfastngf1.551e-01
KvAngf5.220e-06
KvCaB1.024e-06
Navngf3.786e+00
leak8.471e-05
Model and experimental connectivity
Appendix 1—table 38

Structural connection parameters for Neurogliaform cells, based on Bezaire and Soltesz (2013).

https://doi.org/10.7554/eLife.18566.119
Other typeOther cell to NGFNGF to other cell
#Syn.s#Post#Syn.s#Post
Conn.s/Conn.#Loc.Conn.s/Conn.#Loc.
Ivy2810280any dendrite
NGF1710170apical dendrite1710170apical dendrite
O-LM1310130apical dendrite
Pyr12181012181apical dendrite
ECIII52321046any dendrite
Experimental connection constraints
Appendix 1—table 39

Experimental constraints for incoming connections onto Neurogliaform cells (clamp: black=voltage; purple=current).

https://doi.org/10.7554/eLife.18566.120
Pre typeExp. ref.Hold (mV)Erev (mV)Amp. (pA,mV)Diff. %t10-90 (ms)Diff. %τdecay (ms)Diff. %
NGFKarayannis et al., 2010−65.0−11.085.50+0.24.83−1.732.03−46.9
O-LMElfant et al., 2007−50.0−70.018.43−4.01.98−10.211.63+7.6
Appendix 1—table 40

Experimental constraints for outgoing connections from Neurogliaform cells (clamp: black=voltage; purple=current).

https://doi.org/10.7554/eLife.18566.121
Post typeExp. ref.Hold (mV)Erev (mV)Amp. (pA,mV)Diff. %t10-90 (ms)Diff. %τdecay (ms)Diff. %
NGFKarayannis et al., 2010−65.0−11.085.50+0.24.83−1.732.03−46.9
PyrPrice et al., 2008−50.0−89.05.25+7.115.48−3.932.73−34.5
Model synapse parameters
Appendix 1—table 41

Model synaptic parameters for Neurogliaform cells in the control network.

https://doi.org/10.7554/eLife.18566.122
TypeOther cell to NGFNGF to other cell
Erev (mV)Gmax (nS)𝛕rise (ms)𝛕decay (ms)Erev (mV)Gmax (nS)𝛕rise (ms)𝛕decay (ms)
Ivy−60.05.700e-052.903.10
NGF−60.01.600e-043.1042.00−60.01.600e-043.1042.00
O-LM−60.09.800e-051.3010.20
Pyr−60.06.500e-059.0039.00
ECIII0.03.500e-032.006.30
Physiological characterization of model connections
Appendix 1—table 42

Model synaptic properties under voltage clamp at −50 mV with physiological reversal potentials

https://doi.org/10.7554/eLife.18566.123
TypeOther cell to NGFNGF to other cell
Hold (mV)Erev (mV)Amp. (pA)t10-90 (ms)𝛕decay (ms)Hold (mV)Erev (mV)Amp. (pA)t10-90 (ms)𝛕decay (ms)
Ivy−50.0−60.05.481.886.42
NGF−50.0−60.017.525.6714.32−50.0−60.017.525.6714.32
O-LM−50.0−60.09.141.9811.63
Pyr−50.0−60.01.1065.580.00
ECIII−50.00.0324.352.138.80
Appendix 1—figure 16
Connections onto (A) and (B) from model Neurogliaform cells, under voltage clamp at −50 mV with physiological reversal potentials.
https://doi.org/10.7554/eLife.18566.124

O-LM cell: feed back cell (1640 cells)

Model and experimental electrophysiology
Appendix 1—figure 17
O-LM (A) model and (B) experimental current sweep.

(C) Firing rates of model and experimental cells.

https://doi.org/10.7554/eLife.18566.125
Appendix 1—table 43

Model O-LM cell electrophysiological properties.

https://doi.org/10.7554/eLife.18566.126
PropertyValue
RMP−68.0 mV
Input Resistance267.7 MΩ
Sag Amplitude26.5 mV
Sag Tau42.5 ms
Membrane Tau22.7 ms
Rheobase20.0 pA
ISI66.9 ms
Threshold−37.8 mV
Spike Amplitude42.6 mV
Slow AHP Amplitude34.6 mV
Appendix 1—table 44

Model O-LM cell ion channels and conductance at highest density location in cell.

https://doi.org/10.7554/eLife.18566.127
ChannelGmax (S/cm2)
HCNolm5.000e-04
Kdrfast1.174e-01
KvAolm4.950e-03
Nav2.340e-02
leak1.000e-05
Model and experimental connectivity
Appendix 1—table 45

Structural connection parameters for O-LM cells, based on Bezaire and Soltesz (2013).

https://doi.org/10.7554/eLife.18566.128
Other typeOther cell to O-LMO-LM to other cell
#Syn.s#Post#Syn.s#Post
Conn.s/Conn.#Loc.Conn.s/Conn.#Loc.
Axo71071apical dendrite
Bis3910390any dendrite1110107apical dendrite
CCK+B208160any dendrite8810878apical dendrite
Ivy136101360any dendrite
NGF2810283apical dendrite
O-LM61060basal dendrite61060basal dendrite
Pyr237937137basal dendrite15201015195apical dendrite
PV+B2710269apical dendrite
SC-A101097apical dendrite
Experimental connection constraints

Note:No experimental constraints available for incoming synapses to O-LM cells.

Appendix 1—table 46

Experimental constraints for outgoing connections from O-LM cells (clamp: black=voltage; purple=current).

https://doi.org/10.7554/eLife.18566.129
Post typeExp. ref.Hold (mV)Erev (mV)Amp. (pA,mV)Diff. %t10-90 (ms)Diff. %𝛕decay (ms)Diff. %
NGFElfant et al., 2007−50.0−70.018.43−4.01.98−10.211.63+7.6
PyrMaccaferri et al., 2000−70.07.024.35−6.34.68−24.618.88−9.3
SC-AElfant et al., 2007−50.0−70.017.06−12.54.07+114.530.08−3.6
Model synapse parameters
Appendix 1—table 47

Model synaptic parameters for O-LM cells in the control network.

https://doi.org/10.7554/eLife.18566.130
TypeOther cell to O-LMO-LM to other cell
Erev (mV)Gmax (nS)𝛕rise (ms)𝛕decay (ms)Erev (mV)Gmax (nS)𝛕rise (ms)𝛕decay (ms)
Axo−60.01.200e-040.7310.00
Bis−60.02.000e-051.008.00−60.01.100e-040.6015.00
CCK+B−60.07.000e-041.008.00−60.01.200e-030.7320.20
Ivy−60.05.700e-052.903.10
NGF−60.09.800e-051.3010.20
O-LM−60.01.200e-030.257.50−60.01.200e-030.257.50
Pyr0.02.000e-040.300.60−60.03.000e-040.1311.00
PV+B−60.01.100e-030.257.50
SC-A−60.01.500e-040.0729.00
Physiological characterization of model connections
Appendix 1—table 48

Model synaptic properties under voltage clamp at −50 mV with physiological reversal potentials

https://doi.org/10.7554/eLife.18566.131
TypeOther cell to O-LMO-LM to other cell
Hold (mV)Erev (mV)Amp. (pA)t10-90 (ms)𝛕decay (ms)Hold (mV)Erev (mV)Amp. (pA)t10-90 (ms)𝛕decay(ms)
Axo−50.0−60.04.762.5512.03
Bis−50.0−60.01.861.788.13−50.0−60.06.312.7017.05
CCK+B−50.0−60.054.981.359.05−50.0−60.040.323.1028.42
Ivy−50.0−60.05.322.106.33
NGF−50.0−60.09.141.9811.63
O-LM−50.0−60.078.691.059.30−50.0−60.078.691.059.30
Pyr−50.00.017.470.601.53−50.0−60.00.543.7014.10
PV+B−50.0−60.035.531.6510.18
SC-A−50.0−60.07.913.9029.83
Appendix 1—figure 18
Connections onto (A) and (B) from model O-LM cells, under voltage clamp at −50 mV with physiological reversal potentials.
https://doi.org/10.7554/eLife.18566.132

PV+ basket cell: fast-spiking somatic inhibitor (5530 cells)

Model and experimental electrophysiology
Appendix 1—figure 19
PV+ Basket (A) model and (B) experimental current sweep.

(C) Firing rates of model and experimental cells.

https://doi.org/10.7554/eLife.18566.133
Appendix 1—table 49

Model PV+ Basket cell electrophysiological properties.

https://doi.org/10.7554/eLife.18566.134
PropertyValue
RMP−65.0 mV
Input Resistance52.1 MΩ
Sag Amplitude
Sag Tau
Membrane Tau7.0 ms
Rheobase300.0 pA
ISI151.4 ms
Threshold−36.7 mV
Spike Amplitude90.7 mV
Slow AHP Amplitude41.4 mV
Appendix 1—table 50

Model PV+ Basket cell ion channels and conductance at highest density location in cell.

https://doi.org/10.7554/eLife.18566.135
ChannelGmax (S/cm2)
CavL5.000e-03
CavN8.000e-04
KCaS2.000e-06
Kdrfast1.300e-02
KvA1.500e-04
KvCaB2.000e-07
Navaxonp1.500e-01
leak1.800e-04
Model and experimental connectivity
Appendix 1—table 51

Structural connection parameters for PV+ Basket cells, based on Bezaire and Soltesz (2013).

https://doi.org/10.7554/eLife.18566.136
Other typeOther cell to PV+BPV+B to other cell
#Syn.s#Post#Syn.s#Post
Conn.s/Conn.#Loc.Conn.s/Conn.#Loc.
Axo10110soma
Bis1610160any dendrite16115soma
CCK+B12896any dendrite25124soma
Ivy2410240any dendrite13112soma
O-LM81080apical dendrite
Pyr42431272apical dendrite9581110533soma
PV+B39139soma39139soma
SC-A211soma
CA36047212094any dendrite
Experimental connection constraints
Appendix 1—table 52

Experimental constraints for incoming connections onto PV+ Basket cells (clamp: black=voltage; purple=current).

https://doi.org/10.7554/eLife.18566.137
Pre typeExp. ref.Hold (mV)Erev (mV)Amp. (pA,mV)Diff. %t10-90 (ms)Diff. %𝛕decay (ms)Diff. %
PyrLee et al., 2014−60.00.015.09−67.70.28−72.51.83−55.7
PV+BCobb et al., 1997−59.0−70.00.29+14.92.67+105.814.72−45.5
Appendix 1—table 53

Experimental constraints for outgoing connections from PV+ Basket cells (clamp: black=voltage; purple=current).

https://doi.org/10.7554/eLife.18566.138
Post typeExp. ref.Hold (mV)Erev (mV)Amp. (pA,mV)Diff. %t10-90 (ms)Diff. %𝛕decay (ms)Diff. %
BisCobb et al., 1997−55.0−70.00.27−27.50.47−52.511.85+30.4
PyrSzabadics et al., 2007−70.0−26.091.94−13.90.50−5.76.70+4.7
PV+BCobb et al., 1997−59.0−70.00.29+14.92.67+105.813.45−45.5
Model synapse parameters
Appendix 1—table 54

Model synaptic parameters for PV+ Basket cells in the control network.

https://doi.org/10.7554/eLife.18566.139
TypeOther cell to PV+BPV+B to other cell
Erev (mV)Gmax(nS)𝛕rise (ms)𝛕decay (ms)Erev (mV)Gmax (nS)𝛕rise (ms)𝛕decay (ms)
Axo−60.01.200e-040.292.67
Bis−60.09.000e-030.292.67−60.02.900e-030.180.45
CCK+B−60.09.000e-030.434.49−60.01.200e-030.292.67
Ivy−60.07.000e-042.903.10−60.01.600e-040.292.67
O-LM−60.01.100e-030.257.50
Pyr0.07.000e-040.070.20−60.02.000e-040.306.20
PV+B−60.01.600e-030.084.80−60.01.600e-030.084.80
SC-A−60.06.000e-040.292.67
CA30.02.200e-042.006.30
Physiological characterization of model connections
Appendix 1—table 55

Model synaptic properties under voltage clamp at −50 mV with physiological reversal potentials

https://doi.org/10.7554/eLife.18566.140
TypeOther cell to PV+BPV+B to other cell
Hold (mV)Erev (mV)Amp. (pA)t10-90 (ms)𝛕decay (ms)Hold (mV)Erev (mV)Amp. (pA)t10-90 (ms)𝛕decay (ms)
Axo−50.0−60.01.080.453.13
Bis−50.0−60.0429.340.574.13−50.0−60.024.450.170.73
CCK+B−50.0−60.0523.110.685.70−50.0−60.011.310.423.08
Ivy−50.0−60.051.352.056.75−50.0−60.01.440.553.13
O-LM−50.0−60.035.531.6510.18
Pyr−50.00.014.750.251.77−50.0−60.020.560.506.70
PV+B−50.0−60.013.940.235.25−50.0−60.013.940.235.25
SC-A−50.0−60.05.710.423.08
CA3−50.00.019.712.388.78
Appendix 1—figure 20
Connections onto (A) and (B) from model PV+ Basket cells, under voltage clamp at −50 mV with physiological reversal potentials.
https://doi.org/10.7554/eLife.18566.141

Schaffer collateral-associated cell: regular-spiking dendritic inhibitor (400 cells)

Model and experimental electrophysiology
Appendix 1—figure 21
Schaffer Collateral-Associated (A) model and (B) experimental current sweep.

(C) Firing rates of model and experimental cells.

https://doi.org/10.7554/eLife.18566.142
Appendix 1—table 56

Model Schaffer Collateral-Associated cell electrophysiological properties.

https://doi.org/10.7554/eLife.18566.143
PropertyValue
RMP−70.5 mV
Input Resistance300.0 MΩ
Sag Amplitude12.9 mV
Sag Tau41.7 ms
Membrane Tau28.9 ms
Rheobase60.0 pA
ISI115.9 ms
Threshold−36.6 mV
Spike Amplitude80.3 mV
Slow AHP Amplitude35.2 mV
Appendix 1—table 57

Model Schaffer Collateral-Associated cell ion channels and conductance at highest density location in cell.

https://doi.org/10.7554/eLife.18566.144
ChannelGmax (S/cm2)
CavL1.000e-03
CavN2.000e-05
HCN7.000e-05
KCaS1.000e-06
Kdrfast6.000e-05
KvA1.000e-04
KvCaB7.000e-06
KvGroup2.200e-03
Navcck4.000e-02
leak2.857e-05
Model and experimental connectivity
Appendix 1—table 58

Structural connection parameters for Schaffer Collateral-Associated cells, based on Bezaire and Soltesz (2013).

https://doi.org/10.7554/eLife.18566.145
Other typeOther cell to SC-ASC-A to other cell
#Syn.s#Post#Syn.s#Post
Conn.s/Conn.#Loc.Conn.s/Conn.#Loc.
Axo4622any dendrite
Bis1710170any dendrite6633any dendrite
CCK+B278216any dendrite546324any dendrite
Ivy102101020any dendrite446264any dendrite
O-LM4010400apical dendrite
Pyr1053315apical dendrite
PV+B24124soma
CA3194023880any dendrite
ECIII57321146any dendrite
Experimental connection constraints
Appendix 1—table 59

Experimental constraints for incoming connections onto Schaffer Collateral-Associated cells (clamp: black=voltage; purple=current).

https://doi.org/10.7554/eLife.18566.146
Pre typeExp. ref.Hold (mV)Erev (mV)Amp. (pA,mV)Diff. %t10-90 (ms)Diff. %𝛕decay (ms)Diff.%
O-LMElfant et al., 2007−50.0−70.017.06−12.54.07+114.530.08−3.6
SC-APawelzik et al., 2002−58.0−70.01.53−405.62.35−41.322.90−33.2
Appendix 1—table 60

Experimental constraints for outgoing connections from Schaffer Collateral-Associated cells (clamp: black=voltage; purple=current).

https://doi.org/10.7554/eLife.18566.147
Post typeExp. ref.Hold (mV)Erev (mV)Amp. (pA,mV)Diff. %t10-90 (ms)Diff. %𝛕decay (ms)Diff. %
PyrLee et al., 2010−70.0−26.052.42−12.91.63+13.68.55+3.0
SC-APawelzik et al., 2002−58.0−70.01.53−405.62.35−41.327.98−33.2
Model synapse parameters
Appendix 1—table 61

Model synaptic parameters for Schaffer Collateral-Associated cells in the control network.

https://doi.org/10.7554/eLife.18566.148
Other cell to SC-ASC-A to other cell
TypeErev (mV)Gmax (nS)𝛕rise (ms)𝛕decay (ms)Erev (mV)Gmax (nS)𝛕rise (ms)𝛕decay (ms)
Axo−60.06.000e-040.424.99
Bis−60.08.000e-040.292.67−60.06.000e-040.424.99
CCK+B−60.07.000e-040.434.49−60.08.500e-040.424.99
Ivy−60.03.700e-052.903.10−60.08.500e-040.424.99
O-LM−60.01.500e-040.0729.00
Pyr0.04.050e-040.300.60
PV+B−60.06.000e-040.292.67
CA30.03.000e-042.006.30
ECIII0.04.500e-042.006.30
Physiological characterization of model connections
Appendix 1—table 62

Model synaptic properties under voltage clamp at −50 mV with physiological reversal potentials

https://doi.org/10.7554/eLife.18566.149
TypeOther cell to SC-ASC-A to other cell
Hold (mV)Erev (mV)Amp. (pA)t10-90 (ms)𝛕decay (ms)Hold (mV)Erev (mV)Amp. (pA)t10-90 (ms)𝛕decay (ms)
Axo−50.0−60.024.001.006.13
Bis−50.0−60.050.350.704.10−50.0−60.026.431.026.20
CCK+B−50.0−60.049.550.705.38−50.0−60.033.811.056.90
Ivy−50.0−60.03.092.226.88−50.0−60.046.620.855.58
O-LM−50.0−60.07.913.9029.83
Pyr−50.00.017.420.683.05
PV+B−50.0−60.05.710.423.08
CA3−50.00.027.102.359.13
ECIII−50.00.031.823.3810.47
Appendix 1—figure 22
Connections onto (A) and (B) from model Schaffer Collateral-Associated cells, under voltage clamp at −50 mV with physiological reversal potentials.
https://doi.org/10.7554/eLife.18566.150

Inhibitory connectivity

Although the connectivity of the hippocampal CA1 network was assessed in Bezaire and Soltesz (2013), detailed connectivity estimates were only made for pyramidal cells, while the convergence onto inhibitory cells and especially the inhibitory-inhibitory connections were only estimated at a high level due to lack of sufficiently specific experimental data. Here, we performed additional calculations and made use of additional data to arrive at specific estimates for each interneuron type.

First, we gathered previously published morphology data about each interneuron type. As there were no data available for ivy and neurogliaform cells, we performed the necessary experiments in our lab by filling ivy and neurogliaform cells in hippocampal CA1 slices from Wistar rats and then measuring their somatic area and dendrites. While we performed this experiment in slices, these cell types have a relatively compact morphology, allowing us to characterize a significant amount of their dendritic extent. The data for ivy and neurogliaform cells are available in Appendix 1—table 63. The somatic and dendritic lengths used for each interneuron type in this work are given in Appendix 1—table 64.

Next, we gathered data about the somatic and dendritic synaptic densities for each cell class (Appendix 1—table 65) and then multiplied the areas and lengths by their respective synaptic densities to arrive at the estimated excitatory and inhibitory synapses on the somata and dendrites of each cell types (Appendix 1—table 66). Finally, using our previous estimate of the excitatory and inhibitory boutons available for synapsing on interneurons in CA1 Bezaire and Soltesz (2013), we evenly distributed the boutons available for synapsing on each neurite type in each layer across the available postsynaptic densities of that neurite type and transmitter type, so that each interneuron received approximately the same level of coverage of its incoming synapses, while respecting specific observations about interneuron connectivity, such as that CCK+ basket cells have never been observed to receive direct monosynaptic excitement from local CA1 pyramidal cells Lee et al. (2010) and O-LM cells receive almost all of their excitatory connections from local collaterals Blasco-Ibáñez and Freund (1995). The final connectivity between each cell type, including interneurons, is given in the manuscript, Table 1.

Appendix 1—table 63

Measured dendritic lengths and somatic diameters for ivy and neurogliaform cells from the hippocampal CA1 area in Wistar rats, with calculation of somatic surface area included. Cells were characterized in our lab and their function has been reported in Krook-Magnuson et al. (2011). Source Data available in Appendix 1—table 1 - Source Data.

https://doi.org/10.7554/eLife.18566.151
Cell typeCell nameDendritic length (μM))Somatic dia-meter (μm)Calculated synap-tic Area (μm2)
# SectionsSOSPSRSLM
Ivy0217–1 DAB 3_2_10 left slice1129.264.51200.6038.91188.5
Ivy9 n23-7 DAB 12_16_09 left+middle slice2002703.2300.345.21604.6
Ivy9 n23-6 DAB 06_10 left slice175.4133.82115.4036.61052.1
Ivy9 n16-3 DAB 12_29_09 left slice1001015.2052.52164.8
IvyAverage51.1549.5751758.675.0751502.5
Neurogliaform9n 12–5 DAB 1_06_091002097.752534.4929.4
Neurogliaform91021 DAB 3_18_10 second,third,fourth from left slice3001230.7780.128615.8
Neurogliaform9d 8–3 DAB 1_15_10 left and right slice2002328.21382.432.2814.3
NeurogliaformAverage001885.5895.8786.5
Appendix 1—table 64

Estimated or observed somatic area and dendritic length. Experimental observations of the dendritic length of broad interneuron classes were used as the basis for these estimations. The relative lengths for PV+ basket cells and axo-axonic cells were further differentiated based on experimental observations in region CA3 Papp et al. (2013). The observations published in Mátyás et al. (2004) for CCK+ basket cells were also applied to the CCK+ Schaffer Collateral-Associated cells, based on the discussion in Mátyás et al. (2004). The data for ivy and neurogliaform cells were based on measurements from filled cells from slices. Due to the compact nature of their morphology, especially the neurogliaform cells, the dendritic lengths within the slices were assumed to comprise most or all of the dendritic extents of those cells. See section below for raw data. The O-LM cell morphological measurements were taken from Blasco-Ibáñez and Freund (1995).

https://doi.org/10.7554/eLife.18566.152
InterneuronSoma area (100 μm2)Dendritic length (μm2)Reference
TotalSOSPSRSLM
Ivy15021934.451.1549.5751758.675.075See below
Neurogliaform7862781.4001885.5895.8See below
PV+ basket3428435914936971877292(Papp et al., 2013)
Bistratified10064347.751074.57248.282369.24655.66(Gulyás et al., 1999)
Axo-axonic232928255706591259337(Papp et al., 2013)
CCK+ basket9666338.311213.92310.613522.61291.18(Mátyás et al., 2004)
SCA9666338.311213.92310.613522.61291.18(Mátyás et al., 2004)
O-LM3007.784165.684165.68000(Blasco-Ibáñez and Freund, 1995)
Appendix 1—table 65

The synaptic densities (# boutons per 100 μm of dendritic length, or # boutons per 100 μm2 of somatic area) on the soma and dendrites of PV cells, given in Gulyás et al. (1999), were applied to the axo-axonic, PV+ basket, and bistratified cells. The synaptic densities of CCK+ cells Mátyás et al. (2004) were applied to the CCK+ basket and the Schaffer Collateral-Associated cells. For the ivy, neurogliaform, and O-LM cells, there were not sufficient experimental data published to constrain the synaptic density, and so an average of all synaptic densities for all cell classes was computed and applied to these cell types.

https://doi.org/10.7554/eLife.18566.153
Dendritic
SomaticSOSPSRSLM
ReferenceExcInhExcInhExcInhExcInhExcInhRef
Ivy21.816.1172.223.7163.338.5193.825.397.431.7Calc. from average
Neurogliaform21.816.1172.223.7163.338.5193.825.397.431.7Calc. from average
PV+ basket40.718.1342.519.234516.1371.218.3132.228.6(Gulyás et al., 1999)
Bistratified40.718.1342.519.234516.1371.218.3132.228.6(Gulyás et al., 1999)
Axo-axonic40.718.1342.519.234516.1371.218.3132.228.6(Gulyás et al., 1999)
CCK+ basket3.416.184.332.552.787.48237.886.558.8(Mátyás et al., 2004)
SCA3.416.184.332.552.787.48237.886.558.8(Mátyás et al., 2004)
O-LM21.816.1172.223.7163.338.5193.825.397.431.7Calc. from average
Appendix 1—table 66

Estimated numbers of excitatory and inhibitory synapses on each cell type, calculated by multiplying the somatic area or dendritic length by the respective synaptic density. About 20% of synapses onto O-LM cells are GABAergic, while at least 60% are from local excitatory collaterals Kispersky et al. (2012). Therefore, we conserved the total (inhibitory + excitatory) synaptic density of O-LM cells as calculated previously, but set 20% of that total to be inhibitory and the rest to be excitatory synapses.

https://doi.org/10.7554/eLife.18566.154
Dendritic
SomaticSOSPSRSLMTotal
RefExcInhExcInhExcInhExcInhExcInhExcInhRef
Ivy326.9242.388128219340844573243651500Calculated
Neurogliaform171.1126.8000036544778732844527761Calculated
PV+ basket1395.2620.1423023786640944946684018215385925Calculated
Bistratified409.4182368120685640879643486718814200868Calculated
Axo-axonic947.9421.31615908193863383139692109741651Calculated
CCK+ basket32.8155.7102439416427128871332111775951922756Calculated
SCA32.8155.7102439416427128871332111775951922756Calculated
O-LM654.6485.26527163200000065271632Calculated

Ion channel descriptions

Appendix 1—table 67

Ion channels included in the model. GHK: based on Goldman-Hodgkin-Katz equation; Q-O: quasi-ohmic; Hyperpol.-act: Hyperpolarization-activated; Nucleo.-gated: Nucleotide-gated; voltage-act.: voltage activated; voltage-dep.: voltage dependent; Calcium-act.: Calcium-activated; Pyr.: pyramidal; NGF: neurogliaform; dist.: distal; prox.: proximal.

https://doi.org/10.7554/eLife.18566.155
IonModel
ChannelDescriptionTypePyramidalAxo-axonicBistratifiedCCK+ BasketIvyNeurogliaformO-LMPV+ BasketS.C.-Assoc.
Cav,LL-type CalciumGHK
Cav,NN-type CalciumQ-O
HCNHyperpol.-act, Cyclic Nucleo.-gatedQ-O
HCNOLMHyperpol.-act, Cyclic Nucleo.-gated for O-LM cellsQ-O
HCNpHyperpol.-act, Cyclic Nucleo.-gated for Pyr. cellsQ-O
KCa,SSmall (SK) Calcium-activated potassiumQ-O
Kdr,fastFast delayed rectifier potassiumQ-O
Kdr,fast,ngfFast delayed rectifier potassium for NGF-family cellsQ-O
Kdr,pDelayed rectifier potassium for Pyr. cellsQ-O
Kv,AA-type voltage-act. potassiumQ-O
Kv,A,dist,pA-type voltage-act. potassium for dist. Pyr. dendritesQ-O
Kv,A,ngfA-type voltage-act. potassium for NGF-family cellsQ-O
Kv,A,olmA-type voltage-act. potassium for O-LM cellsQ-O
Kv,A,prox,pA-type voltage-act. potassium for prox. Pyr. dendritesQ-O
Kv,Ca,BBig (BK) Calcium-act., voltage-dep. potassiumQ-O
Kv,GroupMultiple slower voltage-dep. potassiumQ-O
leakLeakQ-O
NavVoltage-dep. sodiumQ-O
Nav,bisVoltage-dep. sodium for bistratified cellsQ-O
Nav,cckVoltage-dep. sodium for CCK+ cellsQ-O
Nav,ngfVoltage-dep. sodium for NGF-family cellsQ-O
Nav,pVoltage-dep. sodium for Pyr. cellsQ-O
pasLeakQ-O

Appendix—table 67 lists which model ion channels are found in each model cell type. The channels are further described below, summarized from Bezaire (2015), and their equations are included as well. The activation/inactivation curves and current-voltage relations are reproduced here from Bezaire (2015).

Calcium channels

The calcium channels were adapted from previous Soltesz Lab models (Santhakumar et al., 2005; Dyhrfjeld-Johnsen et al., 2007; Morgan and Soltesz, 2008) and include an L-type and N-type channel (Appendix 1—figure 23); the L-type channel does not inactivate and the N-type channel inactivates.

Cav,L Channel

Jaffe et al. (1994) developed this channel based on activation data from CA1 and CA3 hippocampal neurons in adult guinea pigs, at room temperature. It has been further used in many other models implemented by Migliore. The voltage of half-activation was shifted by -10 mV, accounting for ionic differences in the experimental preparation compared to the model condition. It uses the GHK equation to calculate the driving force through the channel, allowing a mild dependence on calcium concentration.

Cav,N Channel

Jaffe et al. (1994) also developed this channel, using the same preparation as for the Cav,L channel. Aradi and Holmes (1999) then modified the channel code, replacing the GHK calculation with a quasi-ohmic calculation of the driving force. In addition, its behavior was altered somewhat compared to previous implementations such as Morgan and Soltesz (2008) and Santhakumar et al. (2005). Their implementations contained a typo in the channel definition that caused its equations to differ from those presented in Aradi and Holmes (1999), and had the effect of reducing the conductance of the channel below its intended magnitude.

At high levels of activation, the channel conductance decreases slightly (Appendix 1—figure 23), a behavior that resulted from replacing of the GHK calculation with a quasi-ohmic expression. However, it may not have too large of an effect since it only happens at very depolarized potentials, potentials that are likely to be achieved only at the peak of a spike.

Appendix 1—figure 23
Calcium channel currents.
https://doi.org/10.7554/eLife.18566.156

HCN channels

The characteristic behavior of HCN channels, their hyperpolarized voltage-dependent activation is captured in these models, but not their cyclic-nucleotide gating. Because they are hyperpolarization-activated, the protocol used to characterize the inactivating of the other channels was used to characterize the activation of the HCN channels. In Appendix 1—figure 24, the differing behavior of the HCN channels can be seen.

HCN Channel

The HCN channel model was based on experiments carried out in CA1 pyramidal cells of Sprague-Dawley rats at room temperature Chen et al. (2001). The original channel model included fast and slow components and used separate, artificial ion definitions for each. We only retained the slow component as the inclusion of the fast component caused a non-physiological, oscillating sag when included in cells. We also reduced the voltage dependence of the slow component slightly to further decrease the oscillation of the sag.

Appendix 1—figure 24
HCN channel currents.
https://doi.org/10.7554/eLife.18566.157
HCNOLM Channel

Saraga et al. (2003) developed this channel model based on data from young Sprague-Dawley rats at a warm room temperature, and Cutsuridis et al. (2010) included it in a model O-LM cell, which we incorporated into our model.

HCNp Channel

Cutsuridis et al. (2010) developed this channel model based experimental data from adult Sprague-Dawley rats at 23o or 33o (Magee, 1998) and included it in a model pyramidal cell.

Potassium channels

We included several potassium channels: delayed rectifier, A-type potassium channels, calcium-dependent potassium channels, and leak channels. We developed multiple variations, which enabled us to tune their thresholds and voltage dependence to the voltage-dependent behavior of the cells into which we placed them.

Delayed rectifier potassium channels

The delayed rectifier models had a voltage-dependent activation component but not an inactivating component.

In Appendix 1—figure 25, the differing behavior of the delayed rectifiers can be seen.

Kdr,fastchannel

Yuen and Durand (1991) originally implemented this model, based on another model of a fast delayed rectifier in a squid axon. They adjusted its parameters so their model cell action potential wave form matched that produced by an experimental mouse cell. Aradi and Holmes (1999) later modified the model to shift the voltage dependence slightly.

Kdr,fast,ngfchannel

To better fit the behavior of experimental neurogliaform cells, we shifted by −10 mV the voltage dependence of the Kdr,fast channels within the neurogliaform and ivy cells, which usually have a higher threshold than other cells.

Appendix 1—figure 25
Delayed rectifier potassium channel currents.
https://doi.org/10.7554/eLife.18566.158
Kdr,pchannel

This channel, implemented by Migliore based on experimental data from hippocampal cells in rat pups at room temperature (Klee et al., 1995), was used in the pyramidal cell model we included in our model (Poolos et al.,2002).

A-type potassium channels

A-type potassium channels are transient and quickly-inactivating; they activate near the action potential threshold, delaying action potential onset, increasing the action potential threshold, and even modulating early repolarization after an action potential(Storm, 1990).

Kv,A Channel

Migliore et al. (1995) developed this channel model, starting with the equations of Borg-Graham (1991) and modifying them to account for Ficker and Heinemann (1992) and Numann et al. (1987) to get the burst behavior they were investigating.

Kv,A,dist,p Channel

Migliore et al. (1995) also developed this channel for use in the CA1 pyramidal cell model of Poolos et al. (2002), based on rat hippocampal pyramidal cell data from Klee et al. (1995).

Appendix 1—figure 26
A-type potassium channel currents.
https://doi.org/10.7554/eLife.18566.159
Kv,A,prox,pchannel

Migliore et al. (1995) also developed this channel model, which differs slightly from the Kv,A,dist,p channel model in timing and voltage.

Kv,A,ngfchannel

We modified the activation and inactivation equations of the KvA A-type channel developed by Migliore et al. (1995) for use in neurogliaform and ivy cells. We offset the voltage dependence by + 10 mV to better fit the high-threshold neurogliaform cell family.

Kv,A,olmchannel

Saraga et al. (2003) developed this channel based on experimental data from the McBain lab and others, and Cutsuridis et al. (2010) used it in the O-LM cell model that we included in our model.

Other potassium channels

Some of the potassium channels included in the model were not voltage dependent in the same way as the above ones. The leak and KCaS channels are not voltage activated nor inactivating. The KvCaB channel was voltage-dependent, but was also Ca2+ gated. KvGroup was voltage-dependent activation but had no inactivating component.

KCa,Schannel

The model of this calcium-activated potassium channel (known as the small or ‘SK’ channel), was developed by Yuen and Durand (1991) and modified by Aradi and Holmes (1999) based on data from Beck et al. (1997), Latorre et al. (1989), Sah (1996), and Lancaster et al. (1991).

Kv,Ca,B Channel

This big (‘BK’) potassium channel, both voltage dependent and calcium-gated, was implemented by Migliore et al. (1995) based on a model from Moczydlowski and Latorre (1983) and has been modified somewhat.

The behavior of these other channels is shown as a function of voltage (Appendix 1—figure 27) and as a function of internal Ca2+ concentration for the calcium-gated potassium channels (Appendix 1—figure 28).

Appendix 1—figure 27
Other potassium channel currents.

Because they didn’t have a voltage-sensitive inactivation component, only the activation curve, which is equivalent to the IV Peak curve, need be shown here.

https://doi.org/10.7554/eLife.18566.160
Appendix 1—figure 28
Calcium-dependent potassium channel dependence on calcium concentration.

(a) The normalized conductance of the channels are plotted as a function of test voltage step and calcium concentration. (b) and (c) The current-voltage relation is shown at several calcium concentrations for (b) KvCaB channel and (c) KCaS channel. Note that the KCaS channel is only active at the highest calcium concentration and is not dependent on voltage (although the voltage continues to set the driving force) when it is active.

https://doi.org/10.7554/eLife.18566.161
Kv,Groupchannel

We implemented a new potassium channel model, starting with the fast delayed rectifier potassium channel, and adjusting the parameters to match the channel behavior of Lien and Jonas (2003) which was presented as a model of the Kv3.1b channel. However, the methods used in Lien and Jonas (2003) are likely to have included multiple potassium channel types in their channel characterization: they added 300 uM of 4-AP, billed as a low enough dose that it would block only Kv3.1b channels, and subtracted the intracellular potential recording of the cell in the presence of that blocker from the control recording and attributed that entire difference to Kv3.1b. However, it is possible to block other potassium channel types with doses of 4-AP as low as 5-10 μM Campanac et al. (2013), so two or more potassium channel types were likely blocked in Lien and Jonas (2003), meaning that more than one potassium channel type contributes to the dynamics of the model channel fit by Lien and Jonas (2003), so we called it ‘KvGroup’.

Leak channel

The leak channel is a very simple, quasi-ohmic model component that employs a non-specific current with a reversal potential set near but depolarized from the potassium reversal potential. Its conductance is inversely related to the membrane resistance.

Sodium channels

As with the potassium channels, we implemented variations of the fast activating, fast-inactivating sodium channel to enable each cell type to achieve a physiologically realistic threshold.

In Appendix 1—figure 29, the differing behavior of the sodium channels can be seen.

Appendix 1—figure 29
Sodium channel voltage dependence.

The normalized conductance of the sodium channel is plotted as (a) a function of test voltage step to show activation and (b) as a function of holding voltage prior to the test step to show inactivation.

https://doi.org/10.7554/eLife.18566.162
Navchannel

Yuen and Durand (1991) originally developed this channel model, which was then modified by Aradi and Holmes (1999).

Nav,bischannel

We depolarized the threshold of the ch_Nav channel slightly, to better fit the bistratified cell behavior. Additionally, the voltage dependence of its activation has been shifted by −5 mV and inactivation by −3 mV. The coefficients have also been adjusted.

Nav,cckchannel

This channel is modified from ch_Nav to have a higher threshold and a slow inactivation component to help realize the spike adaptation observed in experimental CCK+ cells. The voltage dependence of its activation has been shifted by −1 mV and its coefficients have been adjusted as well.

Nav,ngfchannel

We modified this channel from ch_Nav, giving it a much higher threshold suitable for neurogliaform family cells. We offset its activation and inactivation voltage dependencies have been offset by about −19 mV, and the coefficients have been slightly modified as well.

Nav,cckchannel

We modified this channel from the ch_Nav model to have slightly different kinetics, mainly adjusting the coefficients (with only minor adjustments to the voltage dependence) to achieve more realistic CCK+ cell behavior.

Nav,pchannel

Migliore et al. (1999) developed this channel, using it in their pyramidal cell model that we incorporated into our model Poolos et al. (2002). We only implemented the fast inactivating option of this channel, though a slow option was also available in the original implementation Poolos et al. (2002).

Ion channel equations

For each ion channel, the equations used to compute the conductance are laid out below, with the following conventions:

  • Iion gives the specific current/area through the channels (of that type) within a small area of the membrane; it is multiplied by the area of the neuron section it covers to arrive at a total current estimate.

  • Eion gives the reversal potential of the channel

  • The ions flowing through HCN and leak channels (usually a combination of sodium and potassium) are not explicitly specified in the code, so instead they are denoted with an H or leak.

  • The units of variables and quantities are listed to their right in blue text

  • T Co = temperature of model, 34Co in this network

  • T K = temperature of model converted to Kelvins, 307.15K in this network (or 307.16K in some channels that convert the temperature differently).

  • gmaxnScm2 = max. conductance, set when the cell is defined

  • dt ms = time step

  • tinc ms=dt msq10 

  • v mV = membrane potential local to the ion channel

  • F kC = 96.487kCmol (Faraday’s constant)

  • Temperature may be listed with units of Celsius (Co) or units of Kelvin (K) depending on what was used in the original code equation. For equations in this appendix section listing units of Kelvin, the equation has been simplified so as not to show the conversion of temperature from units of Celsius to Kelvin, even though this conversion takes place explicitly in the equation of the code file itself. Ex: equations in the code files including the term 273.15+T (or 273.16+T if converting relative to the freezing point of water rather than the more correct triple point of water), representing the conversion from Celsius to Kelvin (1K1oC=273.15 degrees), are instead displayed with TK in the equations of this section.

  • The universal gas constant R, commonly listed in units of Joules per K*mol, is listed here in the alternative units of Coulomb * Volt per K*mol (CVK1mol1) to be consistent with the units used for the other terms in the equations.

Unless otherwise specified, for each channel, the temperature dependence can be given by:

(4) q10=3T-34Co10Co

The equations specific to each channel are given below. For each channel, first any equations for constants, which only need to be solved once, are listed. Next, equations that need to be solved each time step are listed. Because of the way the NEURON simulator works, some equations may even be solved multiple times per time step. Also, at each instance of the ion channel (at each point within each section of the cell where the ion channel is inserted), the equations need to be solved again.

CavL

The following equations are solved each time step of simulation, as explained above, using the membrane potential v in the section of the membrane where the ion channel is located:

(5) g Scm2=gmax Scm2m2h
(6) IL mAcm2=g Scm2ghk(v mV,[Ca2+]i mM,[Ca2+]o mM)

The activation of channel conductance is represented by m (v mVoltage-dependent activation) and h (calcium-dependent activation). The driving force through the channel, ghk, is calculated using the Goldman-Hodgkin-Katz (GHK) equation. The values of ghk, m and h are also solved each time step, by calculating the following equations.

For ghk:

(7)f mV=25 mV293.15 KT K2(8)ghk mV=f mV(1([Ca2+]i mM[Ca2+]o mM)exp(v mVf mV))v mVf mVexp(v mVf mV)1

For m:

(9)α ms1=15.69 ms1(1.0v+81.5 mV)exp(1.0v+81.5 mV10.0 mV)1.0(10)β ms1=0.29 ms1exp(v10.86 mV)(11)τm ms=1(α ms1+β ms1)(12)m=α ms1τm ms(13)dmdt ms=mmτm ms

For h:

(14) h=0.001 mM0.001 mM+[Ca2+]i mM

CavN

The following equations are solved each time step:

(15) g Scm2=gmax Scm2c2d
(16) ICa mAcm2=g Scm2(v mVECa mV)

The values of c and d are also solved each time step, by calculating the following equations:

(17)αc ms1=0.19 ms1v mV19.88 mVexp(v mV19.88 mV10 mV)1(18)βc ms1=0.046 ms1exp(v20.73 mV)(19)τc ms=1αc ms1+βc ms1(20)c=αc ms1αc ms1+βc ms1
(21)αd ms1=0.00016 ms1exp(v48.4 mV)(22)βd ms1=1exp(v+39 mV10 mV)+1(23)τd ms=1αd ms1+βd ms1(24)d=αd ms1αd ms1+βd ms1
(25) cexp=1exp(tinc msτc ms),dexp=1exp(tinc msτd ms)
(26) c=c+cexp(cc),d=d+dexp(dd)

HCN

Since the HCN channel conducts a mixture of sodium and potassium, the reversal potential EH is set to lie between their reversal potentials, usually around −30 mV. The following equations are solved each time step:

(27) g Scm2=gmax Scm2h2
(28) IH mAcm2=g Scm2(v mVEH mV)

The value of h is also solved each time step, by calculating the following equations:

(29)τslow ms=(801.5 ms+.75172.7 ms1+exp((v mV+59.3 mV)0.83 mV))1q10(30)h=11+exp(v mV+91 mV10 mV)(31)dhdt ms=hhτslow ms

HCNolm

Since the HCN channel conducts a mixture of sodium and potassium, the reversal potential EH is set to lie between their reversal potentials, usually around −30 mV. The following equations are solved each time step of simulation:

(32) g Scm2=gmax Scm2r
(33) IH mAcm2=g Scm2(v mVEH mV)

The value of r is also solved each time, by calculating the following equations:

(34)r=11+exp(v mV+84.1 mV10.2 mV)(35)τr=100 ms+1 msexp(17.9 mV0.116v mV)+exp(1.84 mV+0.09v mV)(36)rexp=1exp(dt msτr ms)(37)drdt ms=rrτr ms

HCNp

A different temperature dependence was calculated for this channel than the default calculation specified at the beginning of the section:

(38) qt=4.5T-33Co10Co

As with the other HCN channels, the reversal potential EH was set to lie between the sodium and potassium reversal potentials, around −30 mV. Then, the following equations are solved each time step:

(39) g Scm2=gmax Scm2l
(40) IH mAcm2=g Scm2(v mVEH mV)

The value of l is also solved each time step, by calculating the following equations:

(41)α=exp(0.0378 mV12.2(v mV+75 mV))(42)β=exp(0.0378 mV12.20.4(v mV+75 mV))(43)l=11+exp(0.0378 mV14(v mV+90 mV))(44)τl ms=βqt0.011 ms1(1+α)(45)dldt ms=llτl ms

KCaS

The following equations are solved each time step of simulation:

(46) g Scm2=gmax Scm2q2
(47) IK mAcm2=g Scm2(v mVEK mV)

The value of q is also solved each time, by calculating the following equations:

(48)αq ms1=15 mM2ms1([Ca2+]i mM)2(49)βq ms1=0.00025 ms1(50)τq ms=1q10(αq ms1+βq ms1)(51)q=αq ms1τq ms
(52)qexp=1exp(dt msτq ms)(53)q=q+qexp(qq)

Kdrfast

The following equations are solved each time step of simulation:

(54) g Scm2=gmax Scm2n4
(55) IK mAcm2=g Scm2(v mVEK mV)

The value of n is also solved each time step, by calculating the following equations:

(56)αn ms1=0.07 ms1v mV+18 mVexp(v mV+18 mV6 mV)1(57)βn ms1=0.264 ms1exp(v mV+43 mV40 mV)(58)τn ms=1αn ms1+βn ms1(59)n=αn ms1αn ms1+βn ms1
(60) nexp=1exp(tinc msτn ms)
(61) n=n+nexp*(n-n)

Kdrfastngf

The following equations are solved each time step of simulation:

(62) g Scm2=gmax Scm2n4
(63) IK mAcm2=g Scm2(v mVEK mV)

The value of n is also solved each time, by calculating the following equations:

(64)αn ms1=0.07 ms1v mV+8 mVexp(v mV+8 mV6 mV)1(65)βn ms1=0.264 ms1exp(v mV+33 mV40 mV)(66)τn ms=1αn ms1+βn ms1(67)n=αn ms1αn ms1+βn ms1
(68) nexp=1exp(tinc msτn ms)
(69) n=n+nexp*(n-n)

Kdrp

A different temperature dependence was implemented for this channel:

(70) qt=1T-24Co10Co

Then, the following equations are solved each time step:

(71) g Scm2=gmax Scm2n
(72) IK mAcm2=g Scm2(v mVEK mV)

The values of n and l are also solved each time step, by calculating the following equations, where the ideal gas constant R is 8.315 CVK1mol1 and the conversion of Volts to milliVolts is given as 1.0e3 VmV1:

(73)n=11+exp(1.0e3 VmV1(3)(v mV13 mV)9.648e4 Cmol18.315 CVK1mol1T K)(74)τn ms=exp(1.0e3 VmV1(3)0.7(v mV13 mV)9.648e4 Cmol18.315 CVK1mol1T K)qt0.02 ms1(1+exp(1.0e3 VmV1  (3)  (v mV13 mV)9.648e4 Cmol18.315 CVK1mol1  T K))(75)dndt ms=nnτn ms

Kdrslow

The following equations are solved each time step of simulation:

(76) g Scm2=gmax Scm2n4
(77) IK mAcm2=g Scm2(v mVEK mV)

The value of n is also solved each time, by calculating the following equations:

(78)αn ms1=0.028 ms1mV1v mV+30 mVexp(v mV+30 mV6 mV)1(79)βn ms1=0.1056 ms1exp(v mV+55 mV40 mV)(80)τn ms=1αn ms1+βn ms1(81)n=αn ms1αn ms1+βn ms1
(82) nexp=1exp(tinc msτn ms)
(83) n=n+nexp*(n-n)

KvA

A different temperature dependence was implemented for this channel:

(84) q10=3T30 C10 C

Then, the following equations are solved each time step of simulation, for every section within the cell, using the membrane potential v in the section of the membrane where the ion channel is located.

(85)g Scm2=gmax Scm2nl(86)IA mAcm2=g Scm2(v mVEK mV)

The values of n and l are also solved each time step, by calculating the following equations, where the ideal gas constant R is 8.315 CVK1mol1 and the conversion of Volts to milliVolts is given as 1.0e3 VmV1:

(87)n=11+exp(1.0e3 VmV1(3)(v mV+33.6 mV)9.648e4 Cmol18.315 CVK1mol1T K)(88)τn ms=exp(1.0e3 VmV1(3)0.6(v mV+33.6 mV)9.648e4 Cmol18.315 CVK1mol1T K)q100.02 ms1(1+exp(1.0e3 VmV1(3)(v mV+33.6 mV)9.648e4 Cmol18.315 CVK1mol1T K))
(89)l=11+exp(1.0e3 VmV14(v mV+83 mV)9.648e4 Cmol18.315 CVK1mol1T K)(90)τl=exp(1.0e3 VmV14(v mV+83 mV)9.648e4 Cmol18.315 CVK1mol1T K)q100.08 ms1(1+exp(1.0e3 VmV14(v mV+83 mV)9.648e4 Cmol18.315 CVK1mol1T K))
(91)dndt ms=nnτn ms(92)dldt ms=llτl ms

KvAdistp

This channel implemented a different temperature dependence than usual:

(93) qt=3T-24Co10Co

Then, the following equations are solved each time step of simulation, for every section within the cell, using the membrane potential v in the section of the membrane where the ion channel is located.

(94)g Scm2=gmax Scm2nl(95)IA mAcm2=g Scm2(v mVEK mV)

The values of n and l are also solved each time, by calculating the following equations, where the ideal gas constant R is 8.315 CVK1mol1 and the conversion of Volts to milliVolts is given as 1.0e3 VmV1:

(96)dndt ms=nnτn ms(97)dldt ms=llτl ms
(98)αn=exp(1.0e3 VmV1(11+exp(v mV+40 mV5 mV)1.8)(v mV+1 mV)9.648e4 Cmol18.315 CVK1mol1T K)(99)βn=exp(1.0e3 VmV1(11+exp(v mV+40 mV5 mV)1.8)0.39(v mV+1 Cmol18.315 CVK1mol1T K)
(100)n=11+αn(101)τn ms=βnqt0.1 ms1(1+αn)
(102)αl=exp(1.0e3 VmV13(v mV+56 mV)9.648e4 Cmol18.315 CVK1mol1T K)(103)l=11+αl(104)τl ms=0.26 msmV1(v mV+50 mV)

KvAngf

This channel used a different form of temperature dependence:

(105) q10=3T-30Co10Co

The following equations are solved each time step of simulation, for every section within the cell, using the membrane potential v in the section of the membrane where the ion channel is located.

(106) g Scm2=gmax Scm2nl
(107) IA mAcm2=g Scm2(v mVEK mV)

The values of n and l are also solved each time, by calculating the following equations:

(108)n=11+exp(1.0e3 VmV1(3)(v mV+23.6 mV)9.648e4 Cmol18.315 CVK1mol1T K)(109)τn ms=exp(1.0e3 VmV1(3)0.6(v mV+23.6 mV)9.648e4 Cmol18.315 CVK1mol1T K)q100.02 ms1(1+exp(1.0e3 VmV1(3)(v mV+23.6 mV)9.648e4 Cmol18.315 CVK1mol1T K))(110)l=11+exp(1.0e3 VmV14(v mV+83 mV)9.648e4 Cmol18.315 CVK1mol1T K)(111)τl ms=exp(1.0e3 VmV14(v mV+83 mV)9.648e4 Cmol18.315 CVK1mol1T K)q100.08 ms1(1+exp(1.0e3 VmV14(v mV+83 mV)9.648e4 Cmol18.315 CVK1mol1T K))
(112)dndt ms=nnτn ms(113)dldt ms=llτl ms

KvAolm

The following equations are solved each time step of simulation, for every section within the cell, using the membrane potential v in the section of the membrane where the ion channel is located:

(114) g Scm2=gmax Scm2ab
(115) IA mAcm2=g Scm2(v mVEK mV)

The values of a and b are also solved each time, by calculating the following equations:

(116)τa ms=5 ms(117)a=11+exp((v mV+14 mV)16.6 mV)
(118)αb ms1=0.000009 ms1exp(v mV26 mV18.5 mV)(119)βb ms1=0.014 ms1exp(v mV+70 mV11 mV)+0.2(120)τb ms=1αb ms1+βb ms1 (121)b=11+exp(v mv+71 mv7.3 mv)
(122) aexp=1exp(dt msτa ms),bexp=1exp(dt msτb ms)
(123) a=a+aexp(aa),b=b+bexp(bb)

KvAproxp

This channel implemented a different temperature dependence:

(124) q10=3T24 C10 C

Then, the following equations are solved each time step:

(125) g Scm2=gmax Scm2nl
(126) IA mAcm2=g Scm2(v mVEK mV)

The values of n and l are also solved each time step, by calculating the following equations, where the ideal gas constant R is 8.315 CVK1mol1 and the conversion of Volts to milliVolts is given as 1.0e3 VmV1:

(127)dndt ms=nnτn ms(128)dldt ms=llτl ms
(129)αn=exp(1.0e3V mV1 (11+exp(v mV+40 mV5 mV)1.5)(v mV11 mV)9.648e4 Cmol18.315 CVK1mol1T K)(130)βn=exp(1.0e3 VmV1(11+exp(v mV+40 mV5 mV)1.5)0.55(v mV11 mV)9.648e4 Cmol18.315 CVK1mol1T K)
(131)n=11+αn(132)τn ms=βnqt0.05 ms1(1+αn)(133)αl=exp(1.0e3 VmV13(v mV+56 mV)9.648e4 Cmol18.315 CVK1mol1T K)(134)l=11+αl(135)τl ms=0.26 msmV1(v mV+50 mV)

KvCaB

The following equations are solved each time step, where the ideal gas constant R is 8.313424 CVK1mol1 and F is Faraday’s constant, 9.648e4 Cmol1. In the code for this channel, F is reported in units of kC rather than C, but it is always multiplied by the voltage in units of mV instead of V, so the conversion is implicit (mV*kC=V*C).

(136)g Scm2=gmax Scm2n(137)IK mAcm2=g Scm2(v mVEK mV)

The value of n is also solved each time, by calculating the following equations that depend on voltage and internal calcium concentration:

(138)αn ms1=[Ca2+]i mM0.28 ms1[Ca2+]i mM+(0.48e3mMexp(20.84F Cmol1v mVR CVK1mol1T K))(139)βn ms1=0.48 ms11+[Ca2+]i mM0.13e6 mMexp(2F Cmol1v mVR CVK1mol1T K)(140)τn ms=1αn ms1+βn ms1(141)n=αn ms1τn ms
(142)nexp=1exp(tinc msτn ms)(143)n=n+nexp(nn)

KvGroup

The following equations are solved each time step of simulation:

(144) g Scm2=gmax Scm2n
(145) IGroup mAcm2=g Scm2(v mVEK mV)

The value of n is also solved each time, by calculating the following equations:

(146)αn ms1=0.0189324 ms1(v mV4.18371 mV)exp((v mV4.18371 mV)6.42606 mV)1(147)βn ms1=0.015857 ms1  exp(v mV25.4834 mV)(148)τn ms=1αn ms1+βn ms1(149)n=αn ms1αn ms1+βn ms1
(150)nexp=1exp(tinc msτn ms)(151)n=n+nexp(nn)

KvM

This channel used a different temperature dependence than usual:

(152) q10=5T35 C10 C

Then, the following equations are solved each time step of simulation, for every section within the cell, using the membrane potential v in the section of the membrane where the ion channel is located.

(153) g Scm2=gmax Scm2m
(154) IM mAcm2=g Scm2(v mVEK mV)

The value of m is also solved each time, by calculating the following equations:

(155)m=11+exp(v mV+40 mV10 mV)(156)τm ms=120 ms+exp(0.03787.4(v mV+42 mV))0.009 ms1(1+exp(0.03787(v mV+42 mV)))(157)dmdt ms=mmτm ms

Leak

The following equations are solved each time step of simulation, for every section within the cell, using the membrane potential v in the section of the membrane where the ion channel is located. Since the leak channel conducts mostly potassium, the reversal potential Eleak is set close to the potassium reversal potential EK.

(158) g Scm2=gmax Scm2
(159) Ileak mAcm2=g Scm2(v mVEleak mV)

Nav

The following equations are solved each time step:

(160)g Scm2=gmax Scm2m3h(161)INa mAcm2=g Scm2(v mVENa mV)

The values of m and h are also solved each time, by calculating the following equations:

(162)αm ms1=0.3 ms1v mV+43 mVexp(v mV+43 mV5 mV)1(163)βm ms1=0.3 ms1v mV+15 mVexp(v mV+15 mV5 mV)1
(164)τm ms=1αm ms1+βm ms1(165)m=αm ms1αm ms1+βm ms1
(166)αh ms1=0.23 ms1exp(v mV+65 mV20 mV)(167)βh .ms1=3.33 ms11+exp(v mV+12.5 mV10 mV)
(168)τh ms=1αh ms1+βh ms1(169)h=αh ms1αh ms1+βh ms1
(170)mexp=1exp(tinc msτm ms),hexp=1exp(tinc msτh ms)(171)m=m+mexp(mm),h=h+hexp(hh)

Navaxonp

Within the ion channel, these equations solve to constants:

(172) qt=1T24 C10 C

Then, the following equations are solved each time step:

(173)g Scm2=gmax Scm2m3h(174)INa mAcm2=g Scm2(v mVENa mV)

The values of m and h are also solved each time step, by calculating the following equations:

(175)InAct=1(176)dmdt ms=mmτm ms(177)dhdt ms=hhτh ms
(178)αm ms1=Ra ms1(v mV+15 mV)1exp((v mV+15 mV)7.2 mV)(179) betam =Rb ms1(v15 mV)1exp((v15 mV)7.2 mV)
(180)τm ms=1(αm ms1+βm ms1)qt(181)m=αm ms1αm ms1+βm ms1
(182)αh=Rd ms1(v mV+30 mV)1exp((v mV+30 mV)1.5 mV)(183)βh=Rg ms1(v30 mV)1exp((v30 mV)1.5 mV)
(184)τh ms=1(αh ms1+βh ms1)qt(185)h=11+exp(v mV+35 mV4 mV)

Navbis

The following equations are solved each time step:

(187) g Scm2=gmax Scm2m3h
(188) INa mAcm2=g Scm2(v mVENa mV)

The values of m and h are also solved each time, by calculating the following equations:

(189)αm ms1=0.2 ms1v mV+38 mVexp(v mV+38 mV5 mV)1(190)βm ms1=0.5 ms1v mV+10 mVexp(v mV+10 mV5 mV)1(191)τm ms=1αm ms1+βm ms1(192)m=αm ms1αm ms1+βm ms1
(193)αh ms1=0.23 ms1exp(v mV+62 mV20 mV)(194)βh ms1=2 ms11+exp(v mV+9.5 mV10 mV)(195)τh ms=1αh ms1+βh ms1(196)h=αh ms1αh ms1+βh ms1
(197) mexp=1exp(tinc msτm ms),hexp=1exp(tinc msτh ms)
(198) m=m+mexp(mm),h=h+hexp(hh)

Navcck

The following equations are solved each time step:

(199) g Scm2=gmax Scm2m3hs
(200) INa mAcm2=g Scm2(v mVENa mV)

The values of m, h (fast inactivation) and s (slow inactivation) are also solved each time, by calculating the following equations:

(201)αm ms1=0.5 ms1v mV+42 mVexp(v mV+42 mV5 mV)1(202)βm ms1=0.3 ms1v mV+13 mVexp(v mV+13 mV5 mV)1(203)τm ms=1αm ms1+βm ms1(204)m=αm ms1αm ms1+βm ms1
(205)αh ms1=0.6 ms1exp(v mV+65 mV20 mV)(206)βh ms1=1.31 ms11+exp(v mV+12.5 mV10 mV)(207)τh ms=1αh ms1+βh ms1(208)h=αh ms1αh ms1+βh ms1
(209)αs ms1=0.003 ms1exp(v +45 mV6 mV)(210)βs ms1=0.005 ms11+exp(v mV+35 mV20 mV)(211)τs ms=1αs ms1+βs ms1(212)s=αs ms1αs ms1+βs ms1
(213) mexp=1exp(tinc msτm ms),hexp=1exp(tinc msτh ms),sexp=1exp(tinc msτs ms)
(214) m=m+mexp(mm),h=h+hexp(hh),s=s+sexp(ss)

Navngf

The following equations are solved each time step:

(215) g Scm2=gmax Scm2m3h
(216) INa mAcm2=g Scm2(v mVENa mV)

The values of m and h are also solved each time, by calculating the following equations:

(217)αm ms1=0.34133 ms1v mV+24 mVexp(v mV+24 mV5 mV)1(218)βm ms1=0.28483 ms1v mV4 mVexp(v mV4 mV5 mV)1(219)τm ms=1αm ms1+βm ms1(220)m=αm ms1αm ms1+βm ms1
(221)αh ms1=0.29648 ms1exp(v mV+64.4184 mV20 mV)(222)βh ms1=3.0931 ms11+exp(v mV+12.1463 mV10 mV)(223)τh ms=1αh ms1+βh ms1(224)h=αh ms1αh ms1+βh ms1
(225) mexp=1exp(tinc msτm ms),hexp=1exp(tinc msτh ms)
(226) m=m+mexp(mm),h=h+hexp(hh)

Navp

This channel model has a different dependence on temperature:

(227) qt=1T24 C10 C

Then, the following equations are solved each time step:

(228) g Scm2=gmax Scm2m3hs
(229) INa mAcm2=g Scm2(v mVENa mV)

The values of m, h, and s are also solved each time step, by calculating the following equations, where the ideal gas constant R is 8.315CVK-1mol-1 and the conversion of Volts to milliVolts is given as 1.0e3 VmV1:

(230)InAct=1(231)dmdt ms=mmτm ms(232)dhdt ms=hhτh ms(233)dsdt ms=ssτs ms(234)αm ms1=Ra ms1(v mV(15) mV)1exp((v mV(15) mV)7.2 mV)(235)βm ms1=Rb ms1(v(15) mV)1exp((v(15) mV)7.2 mV)(236)τm ms=1(αm ms1+βm ms1)qt(237)m=αmαm+βm(238)αh ms1=Rd ms1(v mV(30) mV)1exp((v mV(30) mV)1.5 mV)(239)βh ms1=Rg ms1(v(30) mV)1exp((v(30) mV)1.5 mV)(240)τh=1(αh ms1+βh ms1)qt(241)h=11+exp(v mV+35 mV4 mV)
(242)s=1(1+exp(v mV+43 mV2 mV))+InAct(11(1+exp(v mV+43 mV2 mV)))(243)τs ms=exp(1.0e3 VmV1120.2(v mV+45 mV)9.648e4 Cmol18.315 CVK1mol1 T K)0.0003 ms1(1+exp(1.0e3 VmV112(v mV+45 mV)  9.648e4 Cmol18.315 CVK1mol1  T K))

Data availability

The following data sets were generated
The following previously published data sets were used

References