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

1The hippocampal CA1 area supports diverse cognitive tasks including learning, memory, and spatial processing. These cognitive tasks are thought to require coordination of neuronal activity provided by physiological network oscillations, including the theta rhythm. In rodents, hippocampal theta is a 5–10 Hz oscillation in the local field potential (LFP) and neuronal firing probabilities, occurring during locomotion and in REM sleep. Though several major afferents provide theta-frequency rhythmic input to the CA1 in vivo, recent reports indicate the presence of spontaneous theta-frequency LFP oscillations even in the isolated whole CA1 preparation in vitro. 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.

2Here 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.,). 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. 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.

3In 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 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.

4Details 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. 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, 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 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.

5An 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. Briefly, our pyramidal cell model 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: 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. The remaining 30% of the interneurons were not included in the model due to paucity of quantitative data. We differentiated the interneurons by their electrophysiological profiles, connectivity patterns, synaptic properties, and anatomical abundance. 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 and executed the simulations on several supercomputers. All model results, characterizations, and experimental comparisons are publicly available.

6First, 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, 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 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.

7In addition to theta rhythm, the model network displayed gamma oscillations (25–80 Hz; Figures 3B and 4D), as expected based on in vivo data and in vitro slice data showing 65–75 Hz gamma oscillations arising in response to theta rhythmic network stimulation. 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. 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, the gamma oscillations were not examined further in the present study.

8These 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.

9Although 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.,). 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. The rationale for scaling up inhibition in this way was that, as described in detail in, 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. 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.

10Next, 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. 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, preferentially discharged around the trough 0/360 of the LFP analog theta rhythm (Figure 5A).

11Interneurons 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, 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.

12These results were in line with recent data from the isolated CA1 preparation in vitro 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; note that PV+ cells receive a high number of excitatory inputs on their dendrites compared to other interneuron classes. 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.

13Importantly, 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. 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.

14Does 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. 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.

15Interestingly, muting the relatively rare CA1 pyramidal cell to pyramidal cell excitatory connections alone (each pyramidal cell contacts 197 other pyramidal cells in the CA1; ) 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. 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. 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 ).

16Since the diverse sources of inhibition from the various interneuronal types are believed to enable networks to achieve more complex behaviors, including oscillations, 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 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 and 2 mV in ). 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.

17To 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).

18Taken 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.

19In agreement with previous predictions, 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, 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.

20We produced a biologically detailed, full-scale CA1 network model constrained by extensive experimental data. 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, 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 ). 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, 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. In turn, phase preferential firing may aid information processing tasks by providing order and allowing multiple channels of information to be processed in parallel.

21Importantly, 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 ) 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).

22Is 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. We previously estimated 204,700 pyramidal cells in ipsilateral CA3, giving an estimated 32,750 spikes from ipsilateral CA3 per theta cycle (Equation 2).

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

24Therefore, 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.

25Our 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, setting of spike timing of neuron types during theta, and playing a role in cortical oscillations and memory processes.

26In 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.

27Based 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.

28Because 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.

29While our model is fundamentally a model of the rat CA1 (e.g., in terms of cell numbers and connectivity; see Table 3 in ), 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 and mouse, similar to the experimental results on the phase specific firing in vivo (e.g., awake rat:; awake mouse: ). 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.

30A 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. 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.

31Our 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. 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. 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.

32This 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. 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 that was then used by several groups to study disparate topics including epilepsy, network mechanisms of inhibition and excitability, simulation optimization, and modeling software. 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.,) and others.

33We 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 at (http://doi.org/10.17605/OSF.IO/V4CEH) so the same simulations can easily be repeated with a future, updated model using SimTracker. 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.

34As 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. 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;), the Human Brain Project (currently for neuroanatomical and electrophysiological data and models of the rat neocortex,), and NeuroElectro (for electrophysiological data from all species and brain areas;). 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.

35The 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.

36All 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. 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.

37The 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.

38Model 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).

39The 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).

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

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

42Compute 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;

43Assign 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).

44When 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.

45The cell number and connectivity parameters were exactly as we reported previously in our in-depth quantitative assessment of anatomical data about the CA1. 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.

46The 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.

47After 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.

48To 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. 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.

49As 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, 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.

50We 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 ) that are not included in the isolated CA1 model).

51We analyzed the results of each simulation with standard neural data analysis methods provided by our SimTracker software, RRID:SCR_014735, discussed in, including the spike density function (SDF) of all pyramidal neuron spikes, the periodogram of the SDF, and the spectrogram of the LFP analog. 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.

52We calculated an approximation of the LFP generated by the model neurons based on the method described by. 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.

53We 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. To see how a cell’s spiking activity is related to its SDF, see Figure 4—figure supplement 1.

54To quantify the frequency and power of the oscillations of the network, we computed a one-sided Welch’s Periodogram of the SDF 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 based on the mtspecgramc function from the Chronux toolbox (http://chronux.org/).

55We 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.

56The 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.

57For 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.

58To 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.

59Our 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, 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.

60For 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. Our analyses of these data can be recreated using our publicly available SimTracker tool.

61SimTracker 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.