Interneuronal mechanisms of hippocampal theta oscillations in a full-scale model of the rodent CA1 circuit
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.001Introduction
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, 2004, 2005; Ferraguti 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 2—source data 1
- https://doi.org/10.7554/eLife.18566.006
-
Figure 2—source data 2
- https://doi.org/10.7554/eLife.18566.007
-
Figure 2—source data 3
- https://doi.org/10.7554/eLife.18566.008
-
Figure 2—source data 4
- https://doi.org/10.7554/eLife.18566.009
-
Figure 2—source data 5
- https://doi.org/10.7554/eLife.18566.010
-
Figure 2—source data 6
- https://doi.org/10.7554/eLife.18566.011
-
Figure 2—source data 7
- https://doi.org/10.7554/eLife.18566.012
-
Figure 2—source data 8
- https://doi.org/10.7554/eLife.18566.013
-
Figure 2—source data 9
- https://doi.org/10.7554/eLife.18566.014
-
Figure 2—source data 10
- https://doi.org/10.7554/eLife.18566.015
-
Figure 2—source data 11
- https://doi.org/10.7554/eLife.18566.016
-
Figure 2—source data 12
- https://doi.org/10.7554/eLife.18566.017
-
Figure 2—source data 13
- https://doi.org/10.7554/eLife.18566.018
-
Figure 2—source data 14
- https://doi.org/10.7554/eLife.18566.019
-
Figure 2—source data 15
- https://doi.org/10.7554/eLife.18566.020
-
Figure 2—source data 16
- https://doi.org/10.7554/eLife.18566.021
-
Figure 2—source data 17
- https://doi.org/10.7554/eLife.18566.022
-
Figure 2—source data 18
- https://doi.org/10.7554/eLife.18566.023
-
Figure 2—source data 19
- https://doi.org/10.7554/eLife.18566.024
-
Figure 2—source data 20
- https://doi.org/10.7554/eLife.18566.025
-
Figure 2—source data 21
- https://doi.org/10.7554/eLife.18566.026
-
Figure 2—source data 22
- https://doi.org/10.7554/eLife.18566.027
-
Figure 2—source data 23
- https://doi.org/10.7554/eLife.18566.028
-
Figure 2—source data 24
- https://doi.org/10.7554/eLife.18566.029
-
Figure 2—source data 25
- https://doi.org/10.7554/eLife.18566.030
-
Figure 2—source data 26
- https://doi.org/10.7554/eLife.18566.031
-
Figure 2—source data 27
- https://doi.org/10.7554/eLife.18566.032
An important set of constraints was the electrophysiology and other properties of individual cells and synapses (Figure 2; Figure 2—source data 3–27; 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 1–2). 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.
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 1–3 and Figure 4—source data 1–2). 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.
-
Figure 3—source data 1
- https://doi.org/10.7554/eLife.18566.038
-
Figure 3—source data 2
- https://doi.org/10.7554/eLife.18566.039
-
Figure 3—source data 3
- https://doi.org/10.7554/eLife.18566.040
-
Figure 4—source data 1
- https://doi.org/10.7554/eLife.18566.042
-
Figure 4—source data 2
- 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 1–11 ). 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 0/360 of the LFP analog theta rhythm (Figure 5A).
-
Figure 5—source data 1
- https://doi.org/10.7554/eLife.18566.046
-
Figure 5—source data 2
- https://doi.org/10.7554/eLife.18566.047
-
Figure 5—source data 3
- https://doi.org/10.7554/eLife.18566.048
-
Figure 5—source data 4
- https://doi.org/10.7554/eLife.18566.049
-
Figure 5—source data 5
- https://doi.org/10.7554/eLife.18566.050
-
Figure 5—source data 6
- https://doi.org/10.7554/eLife.18566.051
-
Figure 5—source data 7
- https://doi.org/10.7554/eLife.18566.052
-
Figure 5—source data 8
- https://doi.org/10.7554/eLife.18566.053
-
Figure 5—source data 9
- https://doi.org/10.7554/eLife.18566.054
-
Figure 5—source data 10
- https://doi.org/10.7554/eLife.18566.055
-
Figure 5—source data 11
- https://doi.org/10.7554/eLife.18566.056
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 3–5) 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 1–2), 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—source data 1
- https://doi.org/10.7554/eLife.18566.062
-
Figure 6—source data 2
- 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).
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 GABA postsynaptic events (Price et al., 2005), whereas the model ivy cells only triggered GABA IPSPs (in agreement with a lack of evidence for ivy cell-evoked GABA IPSPs). Could the slow kinetics of GABA IPSPs contribute to the pacing of the theta oscillations? Indeed, when we selectively removed the GABA component of all neurogliaform cell outgoing synaptic connections, theta power was strongly reduced (Figure 6H). To test whether the contribution of the GABA receptors was due to their slow kinetics, we artificially sped up the GABA IPSPs so that they had GABA kinetics but conserved their characteristic large charge transfer. This alteration was implemented by scaling up the GABA synaptic conductance at neurogliaform cell output synapses to achieve a similar total charge transfer as the control GABA 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 GABA component, only a scaled up fast GABA IPSP with a charge transfer equivalent to the mixed GABA synapses. Therefore, muting the neurogliaform cells strongly disrupted the theta oscillations not because the theta oscillations required the slow kinetics of GABA 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).
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).
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).
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 GABA 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 GABA synapses were not necessary for theta as long as their large charge was carried by the fast GABA synapses. However, we had to increase the conductance of the GABA synapse almost 300 times to achieve a similar charge transfer as that conveyed by the GABA synapse. Such a large conductance is not biologically realistic, indicating that the key role for GABA synapses may be to allow the large synaptic charge transfer via a temporal distribution. Indeed, the importance of GABA 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, GABA, 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; 2011a; Ferguson 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 protocolThe 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:
Calculate the distances between every presynaptic cell and postsynaptic cell of the respective types;
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;
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 protocolThe 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 protocolAs 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 protocolWe 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 protocolWe 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 protocolWe 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 protocolTo 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 protocolWe 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 protocolThe 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 protocolFor the GABA-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 protocolTo 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 protocolOur 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.
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 , is calculated as the average membrane potential during a current injection of 0 . 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 (), is the input resistance calculated from the least hyperpolarized current injection level.
Sag amplitude
Units of , 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 , also computed from the most hyperpolarized current injection level, as the time constant required for an equation of the form 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 , computed from the least hyperpolarized current injection level, as the time constant required for an equation of the form to best fit the potential trajectory from the onset of the current injection until the trace reaches steady state.
Rheobase
Units of , 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 , the average time interval between spike threshold time points for the least depolarized current injection level where the cell spiked regularly.
Threshold
Units of , 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 exceeds some cutoff value (Cooper et al., 2003; Metz et al., 2005); in our case the cutoff was 28 . 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 , meaning the derivative of potential with time exceeds two standard deviations of the average (Atherton and Bevan, 2005).
Spike amplitude
Units of , 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 , 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).
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 2–4, 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
Pyramidal cell: principal cell (311500 Cells)
Model and experimental connectivity
Experimental connection constraints
Model synapse parameters
Physiological characterization of model connections
Axo-axonic cell: fast-spiking axonic inhibitor (1470 Cells)
Model and experimental electrophysiology
Model and experimental connectivity
Experimental connection constraints
Note:No experimental constraints available for incoming synapses to Axo-axonic cells.
Model synapse parameters
Physiological characterization of model connections
Bistratified cell: fast-spiking dendritic inhibitor (2210 Cells)
Model and Experimental Electrophysiology
Model and experimental connectivity
Experimental connection constraints
Model synapse parameters
Physiological characterization of model connections
CCK+ basket cell: regular-spiking somatic inhibitor (3600 cells)
Model and experimental electrophysiology
Model and experimental connectivity
Experimental connection constraints
Note:No experimental constraints available for incoming synapses to CCK+ Basket cells.
Model synapse parameters
Physiological characterization of model connections
Ivy cell: late-spiking cell (8810 cells)
Model and experimental electrophysiology
Model and experimental connectivity
Experimental connection constraints
Model Synapse Parameters
Physiological characterization of model connections
Neurogliaform cell: late-spiking feed forward cell (3580 cells)
Model and experimental electrophysiology
Model and experimental connectivity
Experimental connection constraints
Model synapse parameters
Physiological characterization of model connections
O-LM cell: feed back cell (1640 cells)
Model and experimental electrophysiology
Model and experimental connectivity
Experimental connection constraints
Note:No experimental constraints available for incoming synapses to O-LM cells.
Model synapse parameters
Physiological characterization of model connections
PV+ basket cell: fast-spiking somatic inhibitor (5530 cells)
Model and experimental electrophysiology
Model and experimental connectivity
Experimental connection constraints
Model synapse parameters
Physiological characterization of model connections
Schaffer collateral-associated cell: regular-spiking dendritic inhibitor (400 cells)
Model and experimental electrophysiology
Model and experimental connectivity
Experimental connection constraints
Model synapse parameters
Physiological characterization of model connections
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.
Ion channel descriptions
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.
Ca 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 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.
Ca Channel
Jaffe et al. (1994) also developed this channel, using the same preparation as for the Ca 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.
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.
HCN 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.
HCN Channel
Cutsuridis et al. (2010) developed this channel model based experimental data from adult Sprague-Dawley rats at 23 or 33 (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.
Kchannel
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.
Kchannel
To better fit the behavior of experimental neurogliaform cells, we shifted by −10 mV the voltage dependence of the K channels within the neurogliaform and ivy cells, which usually have a higher threshold than other cells.
Kchannel
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).
K 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.
K 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).
Kchannel
Migliore et al. (1995) also developed this channel model, which differs slightly from the K channel model in timing and voltage.
Kchannel
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.
Kchannel
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 Ca gated. KvGroup was voltage-dependent activation but had no inactivating component.
Kchannel
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).
K 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 Ca concentration for the calcium-gated potassium channels (Appendix 1—figure 28).
Kchannel
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.
Nachannel
Yuen and Durand (1991) originally developed this channel model, which was then modified by Aradi and Holmes (1999).
Nachannel
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.
Nachannel
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.
Nachannel
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.
Nachannel
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.
Nachannel
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:
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.
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 = temperature of model, 34 in this network
T = temperature of model converted to Kelvins, 307.15 in this network (or 307.16 in some channels that convert the temperature differently).
= max. conductance, set when the cell is defined
ms = time step
mV = membrane potential local to the ion channel
kC = (Faraday’s constant)
Temperature may be listed with units of Celsius () or units of Kelvin () 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 (or 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 (), are instead displayed with in the equations of this section.
The universal gas constant , commonly listed in units of Joules per K*mol, is listed here in the alternative units of Coulomb * Volt per K*mol () 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:
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 in the section of the membrane where the ion channel is located:
The activation of channel conductance is represented by (v mVoltage-dependent activation) and (calcium-dependent activation). The driving force through the channel, , 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 :
For :
For :
CavN
The following equations are solved each time step:
The values of and are also solved each time step, by calculating the following equations:
HCN
Since the HCN channel conducts a mixture of sodium and potassium, the reversal potential is set to lie between their reversal potentials, usually around −30 mV. The following equations are solved each time step:
The value of is also solved each time step, by calculating the following equations:
HCNolm
Since the HCN channel conducts a mixture of sodium and potassium, the reversal potential is set to lie between their reversal potentials, usually around −30 mV. The following equations are solved each time step of simulation:
The value of is also solved each time, by calculating the following equations:
HCNp
A different temperature dependence was calculated for this channel than the default calculation specified at the beginning of the section:
As with the other HCN channels, the reversal potential was set to lie between the sodium and potassium reversal potentials, around −30 mV. Then, the following equations are solved each time step:
The value of is also solved each time step, by calculating the following equations:
KCaS
The following equations are solved each time step of simulation:
The value of is also solved each time, by calculating the following equations:
Kdrfast
The following equations are solved each time step of simulation:
The value of is also solved each time step, by calculating the following equations:
Kdrfastngf
The following equations are solved each time step of simulation:
The value of is also solved each time, by calculating the following equations:
Kdrp
A different temperature dependence was implemented for this channel:
Then, the following equations are solved each time step:
The values of and are also solved each time step, by calculating the following equations, where the ideal gas constant is and the conversion of Volts to milliVolts is given as :
Kdrslow
The following equations are solved each time step of simulation:
The value of is also solved each time, by calculating the following equations:
KvA
A different temperature dependence was implemented for this channel:
Then, the following equations are solved each time step of simulation, for every section within the cell, using the membrane potential in the section of the membrane where the ion channel is located.
The values of and are also solved each time step, by calculating the following equations, where the ideal gas constant is and the conversion of Volts to milliVolts is given as :
KvAdistp
This channel implemented a different temperature dependence than usual:
Then, the following equations are solved each time step of simulation, for every section within the cell, using the membrane potential in the section of the membrane where the ion channel is located.
The values of and are also solved each time, by calculating the following equations, where the ideal gas constant is and the conversion of Volts to milliVolts is given as :
KvAngf
This channel used a different form of temperature dependence:
The following equations are solved each time step of simulation, for every section within the cell, using the membrane potential in the section of the membrane where the ion channel is located.
The values of and are also solved each time, by calculating the following equations:
KvAolm
The following equations are solved each time step of simulation, for every section within the cell, using the membrane potential in the section of the membrane where the ion channel is located:
The values of and are also solved each time, by calculating the following equations:
KvAproxp
This channel implemented a different temperature dependence:
Then, the following equations are solved each time step:
The values of and are also solved each time step, by calculating the following equations, where the ideal gas constant is and the conversion of Volts to milliVolts is given as :
KvCaB
The following equations are solved each time step, where the ideal gas constant is and is Faraday’s constant, . In the code for this channel, is reported in units of rather than , but it is always multiplied by the voltage in units of instead of , so the conversion is implicit ().
The value of is also solved each time, by calculating the following equations that depend on voltage and internal calcium concentration:
KvGroup
The following equations are solved each time step of simulation:
The value of is also solved each time, by calculating the following equations:
KvM
This channel used a different temperature dependence than usual:
Then, the following equations are solved each time step of simulation, for every section within the cell, using the membrane potential in the section of the membrane where the ion channel is located.
The value of is also solved each time, by calculating the following equations:
Leak
The following equations are solved each time step of simulation, for every section within the cell, using the membrane potential in the section of the membrane where the ion channel is located. Since the leak channel conducts mostly potassium, the reversal potential is set close to the potassium reversal potential .
Nav
The following equations are solved each time step:
The values of and are also solved each time, by calculating the following equations:
Navaxonp
Within the ion channel, these equations solve to constants:
Then, the following equations are solved each time step:
The values of and are also solved each time step, by calculating the following equations:
Navbis
The following equations are solved each time step:
The values of and are also solved each time, by calculating the following equations:
Navcck
The following equations are solved each time step:
The values of , (fast inactivation) and (slow inactivation) are also solved each time, by calculating the following equations:
Navngf
The following equations are solved each time step:
The values of and are also solved each time, by calculating the following equations:
Navp
This channel model has a different dependence on temperature:
Then, the following equations are solved each time step:
The values of , and are also solved each time step, by calculating the following equations, where the ideal gas constant is and the conversion of Volts to milliVolts is given as :
Data availability
-
Simulation results from a network model of the isolated hippocampal CA1 subfield in ratPublicly available at the Open Science Framework.
-
Hippocampal CA1 network in parallel NEURON with spontaneous theta, gamma: full scale & network clampPublicly available at SenseLab (accession no: 187604).
-
Firing pattern of O-LM cells in mouse hippocampal CA1Publicly available at the Open Science Framework.
-
Intracellular, in vitro somatic membrane potential recordings from whole cell patch clamped rodent hippocampal CA1 neuronsPublicly available at the Open Science Framework.