Combining Theory, Model, and Experiment to Explain How Intrinsic Theta Rhythms Are Generated in an In Vitro Whole Hippocampus Preparation without Oscillatory Inputs

Abstract Scientists have observed local field potential theta rhythms (3–12 Hz) in the hippocampus for decades, but understanding the mechanisms underlying their generation is complicated by their diversity in pharmacological and frequency profiles. In addition, interactions with other brain structures and oscillatory drives to the hippocampus during distinct brain states has made it difficult to identify hippocampus-specific properties directly involved in theta generation. To overcome this, we develop cellular-based network models using a whole hippocampus in vitro preparation that spontaneously generates theta rhythms. Building on theoretical and computational analyses, we find that spike frequency adaptation and postinhibitory rebound constitute a basis for theta generation in large, minimally connected CA1 pyramidal (PYR) cell network models with fast-firing parvalbumin-positive (PV+) inhibitory cells. Sparse firing of PYR cells and large excitatory currents onto PV+ cells are present as in experiments. The particular theta frequency is more controlled by PYR-to-PV+ cell interactions rather than PV+-to-PYR cell interactions. We identify two scenarios by which theta rhythms can emerge, and they can be differentiated by the ratio of excitatory to inhibitory currents to PV+ cells, but not to PYR cells. Only one of the scenarios is consistent with data from the whole hippocampus preparation, which leads to the prediction that the connection probability from PV+ to PYR cells needs to be larger than from PYR to PV+ cells. Our models can serve as a platform on which to build and develop an understanding of in vivo theta generation.


Introduction
The goals of mathematical modeling in neuroscience are many and varied. For any particular study, modeling goals need to be clear, as they guide our decisions during model development. Developed models can be used to generate new hypotheses and investigate interactions across different scales. In doing this, it is helpful to consider what is meant by an explanation. Aristotle's doctrine of the four causes (Falcon, 2015)-material, formal, efficient and final, as described and interpreted in Brette (2013)-is useful. The efficient cause is what triggers the phenomenon to be explained; the material cause refers to the physical substrate of the phenomenon; the formal cause is the specific pattern ("balance") responsible for the phenomenon; and the final cause is the function of the phenomenon. Typically, as stated by Brette (2013), theoretical approaches tend to focus on formal and final causes, and experimental approaches on material and efficient causes. Although all four causes may be needed to obtain a complete understanding, considering these causes can serve to clarify modeling goals and guide usage of developed models.
Electrical oscillations, as recorded in electroencephalograms and local field potentials (LFPs), are hallmarks of the brain that are linked to normal and pathologic functioning (Buzsaki, 2006). Thus, it is essential to understand the mechanisms underlying their generation. A large part of the challenge in obtaining mechanisms underlying oscillation generation is the multiscale nature of the brain, with its biological complexity and cellular specifics (Cohen and Gulbinaite, 2014). Various building blocks have long been known (Gjorgjieva et al., 2016), and it is clear that, for example, postinhibitory rebound building blocks contribute to the generation of cortical oscillations (McCormick et al., 2015). However, it is unclear if one can identify essential building block combinations and balances that underlie oscillation generation in the mammalian brain. In the Aristotelian sense, we have some insight into efficient causes (building blocks) for oscillation generation, but we do not yet have formal cause explanations.
Almost 80 years ago, Jung and Kornmuller discovered theta (3-12 Hz) rhythms in the rabbit (Colgin, 2013). These dominant rhythms are associated with memory processing and spatial navigation, present when the animal is actively exploring or during REM sleep. They can be parsed into atropine-resistant or atropine-sensitive types, with higher or lower theta frequencies, respectively (Buzsaki, 2002;Colgin, 2013Colgin, , 2016, and in vitro models of theta rhythms have been developed (Gillies et al., 2002;Traub et al., 2004). Further, low or high theta rhythms were found to be elicited in rats with fearful or social stimuli respectively (Tendler and Wagner, 2015). In the human hippocampus, theta rhythms are linked to similar behaviors (Lega et al., 2012), although it may be the case that they are associated with a wider behavioral repertoire relative to rodents, as they are present without sensory input (Qasim and Jacobs, 2016). Theta rhythms are heavily studied, but with multiple forms, pharmacological sensitivities, and interactions between brain structures, it is challenging to have a clear understanding of their generation.
To explain how theta rhythms are generated, we need to have models that can be mapped onto experiments. As discussed by Colgin (2013), it is traditionally thought that the medial septum (MS) is critical for the generation of theta rhythms, since they are disrupted when the MS is lesioned or inactivated. Indeed, to understand theta rhythms, many studies have explored, characterized, and modeled the interactions between MS and hippocampus (e.g., Brazhnik and Fox 1999;Wang 2002;Hajós et al., 2004;Kocsis and Li, 2004;Manseau et al., 2008;Varga et al., 2008;Hangya et al., 2009). However, the hippocampus can exhibit theta rhythms without the MS (Goutagny et al., 2009). Further, distinct inhibitory cell populations, such as parvalbumin-positive (PV ϩ ) cells, fire at unique phases of the theta rhythm and play important roles in their generation (Varga et al., 2014;. Ultimately, to understand the varied functional roles of these dominant rhythms and how they are modulated and controlled, we need to include cellular aspects and be clear about the particular form of theta. From a mathematical modeling perspective, this reduces to deciding what "parameters, parameters, parameters" (Skinner, 2012) and values to use and how to represent the biological system, given that any mathematical model is an approximation of the biology.
In this article, we develop microcircuit models that are mapped to an in vitro whole hippocampus preparation that spontaneously expresses theta rhythms. We take advantage of theoretical insights and the ability to readily do thousands of network simulations with our developed mathematical models. We present an explanation for intrinsic CA1 theta generation that has elements of efficient, material, and formal causes. It involves building blocks of spike frequency adaptation and postinhibitory rebound in large pyramidal cell populations coupled with fast-firing PV ϩ cells, in which there is a larger connection probability from PV ϩ to pyramidal cells relative to the other way.

Materials and Methods
Here we summarize our overall strategy, the experimental context of the whole hippocampus preparation, and our mathematical models and analyses. We also describe previous and motivating modeling work that the results are built on.

Overall strategy
Our goal is to develop experimentally motivated microcircuit models of a hippocampal CA1 network to provide insight into the mechanisms underlying theta rhythm generation. Our approach is shown in the schematic of Fig. 1, where orange and black arrows refer to links in the present or previous work, respectively. Cellular mathematical models of excitatory and inhibitory cells based on whole-cell patch clamp recordings from the whole hippocampus preparation were previously developed (Ferguson et al., , 2015b. These cellular models were used to generate either excitatory (Ferguson et al., 2015a) or inhibitory  network models with sizes and connectivities as appropriate for the experimental context. A mean field theory (MFT) approach was taken advantage of to determine parameter regimes in the excitatory networks (Ferguson et al., 2015a).
In the present work, we combine these excitatory and inhibitory networks and perform a detailed computational analysis of this network. We investigate the dynamic interplay between these two cell populations and their roles in theta generation. Further theoretical analyses are required to fully understand the network dynamics. Overall, our strategy combines experiment, model development, simulation and theory. Our models bring together network size, connectivity, and cellular characteristics in a closed fashion given the experimental context.

Experimental context
In  developed an in vitro whole hippocampus rodent preparation that spontaneously generates theta (3-12 Hz) rhythms in the CA1 region. By blocking transmission across the septo-temporal axis, they identified multiple theta oscillators in the hippocampus (see Supplementary Fig. 11 in ). Further, the presence of these theta oscillations was not dependent on the CA3 region (see Supplementary Fig.10 in Goutagny et al. [2009]) and required GABA A and AMPA receptors (see Supplementary Table 1 in ). Given this, we estimate that the minimum circuitry required for CA1 theta rhythms is contained in ϳ1 mm 3 . Using known cell densities and approximate volumes of axonal innervation (West et al., 1991;Aika et al., 1994;Sik et al., 1995;Jinno and Kosaka, 2006;Hosseini-Sharifabad and Nyengaard, 2007), we approximate that 30,000 excitatory pyramidal (PYR) cells and 500 PV ϩ cells are involved in the spontaneous generation of theta rhythms in the CA1 region of the hippocampus. This size estimate is thoroughly described in previous modeling work (Ferguson et al., , 2015a. Subsequent work by Amilhon et al. (2015) indicates that networks of PV ϩ and PYR cells could encompass the basic (minimal) units required for theta rhythm generation in the whole hippocampus preparation. This is because optogenetically silencing PV ϩ cells eliminates the theta rhythm, whereas silencing somatostatin-positive inhibitory cells does not. Further, it is the case that PV ϩ cells receive large excitatory postsynaptic currents (EPSCs) relative to PYR cells during ongoing theta rhythms. Also, Figure 1. Overall strategy. The three schematic parts (left, right, lower) of theory, simulation, and experiment/mathematical model development are bidirectionally linked by arrows. Theory refers to mean field theory that was used to constrain the parameter sets to examine in simulations, using cellular models derived from experiment. Simulation refers to the computation of thousands of network simulations done. Experiment/mathematical model development refers to the cellular, Izhikevich-type models that were developed in the experimental context of the whole hippocampus preparation. In the middle is a schematic of the whole hippocampus preparation with an added blue square to illustrate the piece of tissue from the CA1 region of the hippocampus that is being modeled. The hippocampus schematic is adapted from Fig. 1 of Huh et al. (2016). Orange arrows are links in the present work and black arrows, in previous work. Dashed orange arrow from theory to experiment is because it is indirect, as theory previously contributed to simulation. simultaneous recordings from PV ϩ or PYR cells with extracellular field recording of the theta rhythm indicates that the majority of PV ϩ cells fire phasically with the rhythm, whereas PYR cells fire sparsely (Huh et al., 2016). We estimate that 20% or less of the PYR cells are firing, based on Huh et al. (2016). Given that PV ϩ cells receive large EPSCs and that PYR cells fire sparsely, ongoing theta rhythms must necessarily be dependent on a large network effect.
During the ongoing theta rhythm, EPSCs in PYR cells are very small and variable (estimated to be Ͻ20 pA), whereas the inhibitory postsynaptic currents (IPSCs) to the PYR cells are larger (estimated to be ϳ200 pA). Conversely, the PV ϩ cells receive very large EPSCs (estimated to be up to 1000 pA) and smaller IPSCs (ϳ200 pA). These estimates are based on whole-cell current recordings from Huh et al. (2016). Given these estimates, EPSC/IPSC ratios for PYR cells are Ͻ1, and EPSC/IPSC ratios for PV ϩ cells are Ͼ1.
Overall, we aim to determine the conditions under which our network models can produce population bursts at theta (3-12 Hz) frequency, given that there is sparse firing of excitatory PYR cells and non-sparse firing of inhibitory PV ϩ cells, and particular excitatory/inhibitory balances (Huh et al., 2016).

Mathematical models Cell model
Previously developed cellular models are based on experimental data from the in vitro whole hippocampus preparation (Ferguson et al., , 2015b. They use the mathematical model structure developed by Izhikevich (2006Izhikevich ( , 2010, in which the subthreshold behavior and the upstroke of the action potential is captured and a reset mechanism to represent the spike's fast downstroke is used. Despite being relatively simple, parameter choices can be made such that they have a well-defined (albeit limited) relationship to the electrophysiological recordings. The structure has a fast variable representing the membrane potential, V (mV), and a variable for the slow "recovery" current, u (pA). We used a slight modification to be able to reproduce the spike width. The model is given by where C m (pF) is the membrane capacitance, v r (mV) is the resting membrane potential, v t (mV) is the instantaneous threshold potential, v peak (mV) is the spike cutoff value, I shift (pA) is a current that shifts the f-I curve laterally to allow the model to easily capture the rheobase current (for the strongly/weakly adapting models, rheobase current is ϳ0/5 pA, respectively), I syn (pA) represents the synaptic input from the presynaptic cell population (further details below), I other (pA) is an (excitatory) current drive to the network that is not directly modeled through I syn (further details below), a (ms -1 ) is the recovery time constant of the adaptation current, b (nS) describes the sensitivity of the adaptation current to subthreshold fluctuations (greater values couple V and u more strongly, resulting in possible subthreshold oscillations and low-threshold spiking dynamics), c (mV) is the voltage reset value, d (pA) is the total amount of outward minus inward currents activated during the spike and affecting the after-spike behavior, and k (nS/mV) represents a scaling factor. Parameter values for the cell models (strongly and weakly adapting PYR and PV ϩ cell models) are given in Table 1.

Synaptic model
Synaptic input is modeled through a chemical synapse represented by where g (nS) is the maximal synaptic conductance of the synapse from a presynaptic neuron to the postsynaptic neuron, E rev (mV) is the reversal potential of the synapse, and V (mV) is the membrane potential of the postsynaptic cell. The gating variable, s, represents the fraction of open synaptic channels and is given by first-order kinetics (Destexhe et al., 1994;Ermentrout and Terman 2010, p.159): The parameters ␣ (in mM -1 ms -1 ) and ␤ (in ms -1 ) in Eq. 3 are related to the inverse of the rise and decay time constants ( R and D in ms).
[T] represents the concentration of transmitter released by a presynaptic spike. Suppose that the time of a spike is t ϭ t 0 and [T] is given by a square pulse of height 1 mM lasting for 1 ms (until t 1 ). Then, we can represent After the pulse of transmitter has gone, s(t) decays as s(t) ϭ s(t 1 )e Ϫ␤(tϪt 1 ) .
Network models Excitatory, inhibitory, and excitatory-inhibitory network models are illustrated in Fig. 2 for PYR cell networks (top), PV ϩ cell networks (middle), and PYR-PV ϩ cell networks (bottom). Networks have either deterministic or noisy "other" input. If deterministic, then I other is a constant, tonic input to individual cells in the network, where I other is chosen from a normal distribution with mean I app (pA) and SD app (pA). If noisy, then I other ϭ Ϫ g e ͑t͒͑V Ϫ E rev ͒. g e (t) is a stochastic process similar to the Ornstein-Uhlenbeck process as used by Destexhe et al. (2001): where e (t) is an independent Gaussian white noise process of unit SD and zero mean, g e,mean (nS) is the average conductance, e (nS) is the noise SD value, and e is the time constant for excitatory synapses.

Network model parameters, rationale
Random connectivity was used throughout, and the probability of connection is given in Table 2, where it is fixed for PYR or PV ϩ cell networks, as estimated in previous work. Network sizes and synaptic time constants are given in Table 2

PYR-PV ϩ cell networks
A full exploration was done for connectivity between PV ϩ and PYR cell networks (Table 2). We use c PV,PYR and c PYR,PV to refer to the probability of connection from PV ϩ to PYR cells or from PYR to PV ϩ cells, respectively. g pv is the maximal inhibitory synaptic conductance between PV ϩ cells, and a parameter value of 3 nS was used based on our previous PV ϩ cell network modeling   Bartos et al. (2002), for most of the simulations. g pyr-pv (nS) is the maximal excitatory synaptic conductance on PV ϩ cells from PYR cells. A value of 1 nS was used for most of the simulations, as estimated from Papp (2013).
Noisy networks Full connectivity explorations were done for g e,mean (nS) ϭ [0, 1, 2], e (nS) ϭ [0, 0.2, 0.4, 0.6] with g pyr ϭ 0.014 nS, all g pyr values given above, and e ϭ [0, 0.2, 0.4, 0.6] with g e,mean ϭ 0 nS. g pv-pyr ϭ 8.7 nS was used, but additional simulations using values of [6, 6.5, 7. . .., 11.5, 12] were done for two parameter sets with g e,mean ϭ 0 nS; [ e , g pyr , c PYR,PV , c PV,PYR ] ϭ [0.2, 0.084, 0.4, 0.5] and [0.6, 0.014, 0.02, 0.3]. Similarly, g pyr-pv ϭ 3 nS was used and scaled for network size (see below). Additional simulations were done for g pyr-pv values of [0.5, Figure 2. Network schematics. Left: Excitatory PYR cell networks. The large circle with PYR represents a population of individual PYR cells, and the black triangles represent excitatory connections. Middle: Inhibitory, fast-firing cell networks. The large circle with PV ϩ represents a population of individual inhibitory, fast-firing PV ϩ cells, and the black circles represent inhibitory connections. The connections are schematized so that there is clearly not all-to-all coupling. "Other" input is either deterministic or noisy, excitatory synaptic input. Right: Excitatory-inhibitory cell networks. PYR and PV ϩ cell populations combined to create networks of 10,500 cells. However, now the excitatory input to PV ϩ cells comes directly from the PYR cell network. Connectivity between PV ϩ and PYR cells was not fixed, but ranged in values as indicated. A resolution of 0.01 was used from 0 to 0.1, and 0.1 upward.

Analyses
For each network simulation, we defined the population activity as the average membrane potential of all model cells. Then, using the fast Fourier transform (fft), the network frequency (f peak in Hz) is defined as the frequency at which there is a spectral peak in the overall population activity. In this analysis, we disregard the initial transient activity (500 ms).
We defined a population burst based on the distribution of spikes of the PYR cell network. To do so, the total number of spikes within a small bin width were summed, where the bin width was dependent on the average peak frequency: bin width ϭ int͑p 1 ϫ round͕͓͑p 2 ϫ exp ͑ Ϫ p 3 ϫ f peak ϩ p 4 ͒ ϩ p 5 ͔/p 1 ͖͒, where p 1 ϭ 2, p 2 ϭ 2.0264, p 3 ϭ 0.2656, p 4 ϭ 2.9288, and p 5 ϭ 5.7907, such that for peak frequencies ranging from f peak ϭ 3-12 Hz, the bin width ranged from Ϸ23 to 7 ms. In this way, the bin would be smaller for higher frequencies. Then, the total number of spikes per bin width was normalized to its maximum (excluding transient activity within the first 500 ms) so that networks with significantly different levels of activity could be compared. A moving threshold capturing approximately five cycles was set to be the mean ϩ 0.35 SDs of the local normalized distribution. Then, the burst was determined to be the midway point between the increase past threshold and the previous decrease past the threshold (with the requirement that these points are at least 1/͑f peak ϫ 2.5͒ ms apart). If the difference between the peak to trough of the burst is Ͻ0.2, it is no longer considered to be a burst. A population burst is considered to be more robust if the power of the fft is larger or the normalized size of the PYR cell spike distribution is larger. We note that if the population burst is reasonably robust, then the burst frequency as determined from the fft is essentially the same as the inverse of the burst width.
We automated the categorization of our network output for the different parameter sets explored. Specifically, nonfiring cases were considered when there were Ͻ300 spikes per burst bin. If network burst frequencies were within theta frequency ranges, they were further examined to determine their stability. Bursts were considered to be stable if there were at least two occurrences of two consecutive amplitudes decreasing by Ͼ79%. For each burst, we determined the burst width, the number of cells that fired in the burst, and the total number of spikes in the burst. In this way, we can not only track these properties for the network as a whole, but also determine how they change over time. This analysis was based on custom code written in Matlab.
For each simulation, we recorded EPSCs and IPSCs from 100 PYR and 50 PV ϩ cells and chose a subset to analyze. Specifically, we used peakfinder in Matlab and ignored any peaks that were below a value thresholded at an order of magnitude less than the main peaks. The first second was not included in the calculations. We computed averages and SDs of 3 PV ϩ and 5 PYR cells and rounded them to give the reported values.
Simulations were run using the Brian simulator (Goodman and Brette, 2009) on the GPC supercomputer at the SciNet High Performance Computing Consortium (Loken et al., 2010). The initial conditions of our membrane potentials (V) were chosen to be uniform random values from -55 to -65 mV. We used the forward Euler method for integration with a time step of 0.02 ms. PYR cell networks were simulated for 10 s, whereas PYR-PV ϩ cell networks were simulated for 4 s. For noisy simulations, simulations were done with second-order explicit Runge-Kutta numerical integration, with a time step of 0.04 ms. A subset of these simulations was also run with the forward Euler method and compared.

Motivating modeling studies Excitatory (PYR cell) networks, deterministic
In previous work, it was considered whether CA1 PYR cell networks on their own could generate theta rhythms (Ferguson et al., 2015a) or, more specifically, given CA1 PYR cell intrinsic properties, connectivity, and cell numbers, can one obtain theta frequency (3-12 Hz) population bursting as observed in the experimental context? This was directly addressed in Ferguson et al. (2015a). Individual PYR cell models were based on whole cell recordings from the whole hippocampus preparation (Ferguson et al., 2015b). Cells exhibited either weak or strong adaptation (determined by how much their frequency changed over the course of a 1 second long input), and postinhibitory rebound (spiking after being released from a hyperpolarizing current), and these models captured these properties. These PYR cell models were connected in a network (see left schematic of Fig. 2 and model details above) and MFT was used to find parameter regimes in which the network exhibited theta frequency population bursts (Ferguson et al., 2015a). Due to a scaling relationship between cell number, connection probability, and g pyr from the MFT, one can use 10,000, rather than 30,000, PYR cells in the network simulations.
Stable theta frequency population bursts emerged from these PYR cell network models without a phasic drive. Similar to what was predicted from the MFT simulations, larger mean excitatory drive, I app , was required to obtain theta frequency bursts as the SD of the drive, app , increased (Ferguson et al., 2015a). Burst frequencies were determined by three factors: the mean excitatory drive, I app (which increases the burst frequency as it increases); the recurrent synaptic strength, g pyr (which decreases burst frequency as it increases); and the SD of the excitatory drive across the cells, app (which increased frequencies as it increased) (Ferguson et al., 2015a;Ferguson, 2015).
From the parameter sets explored, network output was automatically categorized (see specifics in above sections) such that nonfiring, stable theta frequency bursts, unstable bursts, or other groupings were apparent. In Fig.  3, we show three example outputs, one of which exhibits theta rhythms (bottom), another unstable bursts (top), and the third asynchronous behavior (middle). In the theta bursting parameter regimes, we analyzed our networks to determine how many PYR cells were firing (i.e., active) during the population bursts. We found that Ͼ90% of PYR cells are active when theta rhythms are present. In Table 3, we show the minimal number of active PYR cells for each app value. Except for app ϭ 0, these minimal cases all have frequencies that are below theta but ϳ50% of their 10,000 PYR cells are still active. Also, the number of active PYR cells during bursts increased with increasing I app and also with increasing g pyr . These observations are from simulations that were performed using strongly adapting PYR cell models. Weakly adapting PYR cell models were also used, but specifics are not shown; it is already known from previous MFT studies that more input drive is required to achieve theta frequency population bursts with the reduced cellular adaptation (Ferguson et al., 2015a). Thus, these model outputs suggested that an appropriate balance between spike frequency adaptation and excitatory connectivity in CA1 PYR cell networks could provide an essential mechanism for theta population bursts. However, the majority of PYR cells in the models were active during the population bursts, which is not consistent with experiment. This led us to consider how inhibitory cells may contribute to burst dynamics.

Excitatory-inhibitory (PYR-PV ϩ cell) networks, deterministic
To build networks with both excitatory and inhibitory cells, we first took advantage of previous work in which cellular models for PV ϩ fast-spiking cells were developed based on experimental recordings from the whole hippocampus preparation . Given this work, estimates of EPSCs onto PV ϩ cells of ϳ1000 pA, and that we want PV ϩ cells to fire coherent bursts, we set g pv ϭ 3 nS Skinner and Ferguson 2013).
Next, rather than setting the excitatory drive to PV ϩ cell network as deterministic "other" input as has been done previously in Ferguson et al., (2013) (see middle schematic of Fig. 2), we created networks in which the excitatory drive comes directly from the 10,000 PYR cell network. This is shown in the right schematic of Fig. 2. Because our model is designed to explore oscillatory activity intrinsic to the CA1 region of the hippocampus, input from other regions is not specifically included. We chose example PYR cell networks with distinct firing patterns (nonfiring, stable bursts, etc.) and explored how the connectivity between PYR cells and PV ϩ cells affects network activity (Ferguson, 2015). This limited set of deterministic, excitatory-inhibitory network simulations provided a motivating basis for the expanded set of simulations presented in Results.
In PYR cell networks that exhibit stable bursts, introducing PV ϩ cells does not ensure that the bursts are maintained, but instead requires that the connection probability from PV ϩ to PYR cells (c PV,PYR ) surpasses a critical value. This critical connectivity value depended on the PYR-PV ϩ connection probability (c PYR,PV) , as it could not be drastically lower than this critical value. Interestingly, if stable bursts are maintained, the frequency of the population bursts is always higher in the PYR-PV ϩ cell networks relative to the PYR cell networks alone. We note that for each set of simulations, we did not alter the excitatory drive to the PYR cells, and thus the cells did not increase their firing owing to an external change in the amount of excitation. Rather, this increased burst frequency is due to the postinhibitory rebound spiking of the PYR cells as a response to the inhibitory input from the PV ϩ cells.
Alternatively, if the original PYR cell networks were unstable or nonfiring, stable population bursts in the theta frequency range could emerge with the inclusion of the PV ϩ cell population. In all cases, and in contrast with our PYR cell networks, it was possible to simultaneously obtain theta frequency population bursts and sparse firing of the PYR cells. These observations are illustrated in Fig. 4. Thus, even when oscillations did not exist or were not stable in the PYR cell networks, the influence of PV ϩ cells could lead to stable network rhythm generation and sparse firing. However, this depended on an appropriate balance of connection probabilities between the two populations.
Overall, we found that PYR cell networks could exhibit coherent firing at various frequencies (Ϸ0.5-6 Hz), but when theta frequencies were produced, essentially all PYR cells were recruited to fire in every cycle, a behavior that is not consistent with what is seen in the experimental setting, where PYR cells sparsely fire. However, the inclusion of PV ϩ cells allowed sparse firing to emerge while maintaining theta bursts, but of increased frequency, or generating theta bursts. Because we did not change excitation in any other way, a key contributor to this is postinhibitory rebound firing. Thus, population bursts at theta frequency in the PYR-PV ϩ cell networks with sparse firing depend on a number of factors, including the amount of input that the cells receive at any given point in time. As such, the precise connection probabilities in the models is not the essence, but rather the relative balance between excitation and inhibition. More complete exploration and analyses are needed to untangle this.

Results
We use network models that are minimal but at the same time are constrained by experiment. In this way, we reduce the uncertainty in choices for parameters and parameter values. Detailed reasoning and rationale for our choices are provided in Materials and Methods. We examine whether it is possible to capture the experimental observations, and if so, what are the underlying mechanisms that allow this? To start to address this, we have presented previous work and motivating modeling simulations using deterministic networks. Although limited, these simulations showed that there is an intricate tangling of cellular properties and excitatory and inhibitory balances that underlie the generation of theta population bursts with sparse PYR cell firing in the whole hippocampus preparation. That is, there are many interconnected factors. We thus performed an expanded set of simulations using more biologically realistic input to obtain insight into possible underlying mechanisms.

Noisy networks Excitatory PYR cell networks
We performed simulations with noisy, excitatory input. That is, "other" input in Fig. 2 is given by a stochastic rather than a deterministic process (see Materials and Methods). Because we know that PYR cells do not receive large EPSCs during the endogenous theta rhythm (Huh et al., 2016), we can focus on small mean excitatory conductances. That is, an EPSC of 20 pA into a PYR cell would give an excitatory conductance estimate of Ͻ1 nS, given excitatory reversal potential and resting voltage values. Given this, and that PYR cell recurrent connectivity is minimal, we consider excitatory drive (g e,mean ) values of 0 -2 nS to fully encompass the biological situation in the whole hippocampus preparation.
As would be expected, the patterns in these noisy networks are more variable than the deterministic PYR cell network simulations. However, similar to the deterministic simulations, burst frequency increased with increasing excitatory drive (g e,mean ) and decreased with increasing recurrent synaptic strength (g pyr ). There were no theta frequency bursts for g e,mean ϭ 0, but for larger excitatory drives (g e,mean ϭ 1 or 2), there were theta frequency population bursts except for one case that was outside of the theta frequency range: g pyr ϭ 0.094 nS, g e,mean ϭ 1, e ϭ 0.6. Example output is shown in Fig. 5 with (bottom) or without (top) population bursts. Also shown in Fig. 5 is the PYR cell spike distribution that was used to determine whether bursts were present (see Materials and Methods).
The number of active PYR cells always exceeded 50%. Table 3 provides details of the five networks with the most sparse PYR cell firing. For weakly adapting cells, no theta rhythms were obtained for g e,mean ϭ 0 or 1 but could emerge for larger input values. Because this would be beyond EPSCs values observed in experiments, we did not do further detailed explorations using weakly adapting cells. Thus, similar to the deterministic PYR cell network simulations, theta population bursts are present but never with sparse firing of PYR cells in more realistic, noisy PYR cell networks.

Excitatory-inhibitory PYR-PV ϩ cell networks
A full exploration of connectivities (c PV,PYR and c PYR,PV) with g pyr ϭ 0.014 and g e,mean ϭ 0, 1, and 2 nS was done. We first note that, similar to the deterministic simulations, theta bursts could be present in PYR-PV ϩ cell networks even if PYR cell networks did not have any population or not (top). On the left is the raster plot, and on the right is the normalized spike distribution. Red stars delineate the automatic detection of separate bursts. Parameter values: g e,mean ϭ 0, g pyr ϭ 0.074, e ϭ 0.2 (top row: no theta rhythm); g e,mean ϭ 1, g pyr ϭ 0.074, e ϭ 0.2 (bottom row: 4.4 Hz rhythm).
bursts. An example of this is shown in Fig. 6. A summary of the burst frequencies is shown in the top part of Fig. 7 for two excitatory drives (g e,mean ϭ 0 and 1 nS). Theta frequency bursts encompass colors that range from light blue to orange. Given this, it is clear that PYR-PV ϩ cell networks with g e,mean ϭ 1 nS had burst frequencies that exceeded theta, except for when e ϭ 0 nS, meaning noiseless networks. For g e,mean ϭ 2 nS, the burst frequencies far exceeded theta frequencies (not shown). Also, since larger g pyr values result in an increased burst frequency (see Fig. 8, right), it did not make sense to do additional simulations with g e,mean ϭ 1 or 2 nS. This thus led to a focus on simulations with g e,mean ϭ 0 nS. The full range of g pyr values were simulated along with explorations of PYR to PV ϩ and PV ϩ to PYR cell conductance values (see Materials and Methods).
From a computational analysis that consisted of several thousands of simulations, we were able to obtain a lay of the land in terms of required parameter balances for theta rhythms to occur, as well as their characteristics. This is schematized in the bottom of Fig. 7, where the connectivity ranges that refer to the summarized plots above are indicated by the blue ellipse. As can be seen, they are summarized only for smaller c PYR,PV values. However, as c PYR,PV increases, the burst frequency increases further (not shown). As noticed in the motivating deterministic simulations, and similarly here for the noisy runs, if c PV,PYR is too small, theta bursts are not present. That is, there needs to be enough connectivity from PV ϩ to PYR cells to have postinhibitory rebound firing of an appropriate amount in PYR cells for theta bursts to occur. An example of how theta bursts are lost (when c PYR,PV is also small) is shown in Fig. 8, left.

An essence of theta rhythm generation and experimental data matching
From our noisy and motivating deterministic simulations, we are able to distinguish two scenarios, A and B, by which theta rhythms occur. It is because of the wide swaths of simulations and analyses of them that we are able to identify these different scenarios as shown in Fig.  7. As schematized, the transition between scenarios is not meant to be literal but illustrative, to consider the several other parameters that would affect the exact transitions. The difference between these two scenarios lies in how active the PYR and PV ϩ cell populations are, which in turn affects the EPSCs and IPSCs received by them.
So how do theta frequency bursts emerge in PYR-PV ϩ cell networks? We first point out that g e,mean ϭ 0 nS with nonzero e is the appropriate "other" input to use to have theta frequency population bursts. This is because of results from PYR cell networks and knowledge of EPSC values to PYR cells in experiments. Also, as noted above, if c PV,PYR is not large enough, there will be no theta rhythms (Fig. 8, left). This immediately indicates the importance of postinhibitory rebound spiking in PYR cells to generate theta bursts. However, while c PV,PYR cannot be zero (since this would mean that PV ϩ cells are not cou- pled to PYR cells and one would have PYR cell network output, which we know has no theta bursts for g e,mean ϭ 0), the exact nonzero value needed for theta bursts to emerge also depends on c PYR,PV . The example shown in Fig. 8 (left) is for a small value of c PYR,PV , and the PV ϩ cells are not able to fire much or coherently. If c PV,PYR is decreased beyond a critical value, the PYR-PV ϩ cell networks do not generate theta frequency bursts, as the PYR cells do not exhibit enough postinhibitory rebound to enable robust burst firing in the theta range. The model results are consistent with previous experimental findings from , suggesting the importance of postinhibitory rebound firing in the generation of theta rhythms. Once one is within parameter balance regimes with robust theta frequency bursts, the frequency increases with increasing c PYR,PV and less so for increasing c PV,PYR , as determined from an overall examination of the simulations. This observation is particularly apparent in the upper summary plot of Fig. 7 (g e,mean ϭ 0 and e ϭ 0.6). We illustrate this with a slanted box (burst frequency increases) in the schematic of Fig. 7. The changing burst frequency can be seen if g pyr-pv rather than c PYR,PV is modified, as shown in Fig. 9.
The particular parameter balances that allow the generation of theta population bursts affect not only the specific frequency of the population burst but also how robust it is, that is, how easily discernible it is (see Materials and Methods). This interdependence of cellular and synaptic properties affects how much the PYR and PV ϩ cell populations fire. It is clear that the addition of the PV ϩ cell population is what allows the PYR cell population to fire sparsely. However, how sparse the firing of PYR and PV ϩ cells are depends on where the parameter balance lies. From our simulations, we observe that as c PYR,PV decreases, the firing of PV ϩ cells becomes more sparse and PYR cells become less sparse within a theta population burst. This is illustrated by the box in the Fig. 7 bottom schematic. It makes sense that PV ϩ cells would Color plots shows how burst frequency changes for a range of input to the PYR cell population and for a range of connectivities between PYR and PV ϩ cell populations. Note that c PYR,PV is shown for a much smaller range than c PV,PYR . Moving horizontally, one sees that increasing e typically leads to an increase in burst frequency. It is also clear that beyond g e,mean ϭ 0, theta frequency rhythms are not present if e is nonzero. Other parameter values: g pyr (nS), g pyr-pv (nS), g pv-pyr (nS) ϭ [0.014, 3, 8.7]. Bottom: Theta rhythm generation overview. This schematic summarizes the balances in the generation of theta rhythms and their characteristics in the network models. Two scenarios, A and B, can be identified based on connectivity balances (c PV,PYR and c PYR,PV) and are approximately delineated by the squiggly gray dashed lines. Note that these separations are illustrative, as the exact connectivity boundary value will depend on the other parameters in the models. However, it is clear from the many simulations done and analyzed that one is able to differentiate these regions. Network models characteristics (frequency, and PV ϩ and PYR cell firings) are given in boxes with orange arrows. Scenarios A and B are differentiated by their EPSC/IPSC ratios to PV ϩ cells, as given in magenta text. fire less, as they are receiving less excitatory drive with a reduced c PYR,PV . However, that PYR cells would fire less sparsely when the PV ϩ cells are firing less indicates that how much the PYR cells fire is not only dependent on postinhibitory rebound firing. There are clearly different balances going on. It is these different balances as brought forth from our thousands of simulations and analyses that allowed us realize that one could distinguish two scenarios by which theta rhythms emerge. Specifically, in scenario A, the PV ϩ cells fire less sparsely and the PYR cells fire more sparsely than in scenario B.
In Table 4, we present analyses regarding the average number of cells firing per burst as well as the average number of spikes per cell per burst for several parameter sets. From this table, the two different scenarios described above regarding the relative firing of PYR and PV ϩ cells can be appreciated in more detail. The third and fourth columns of Table 4 show examples where the proportion of actively firing PYR cells are larger or smaller, and similarly for PV ϩ cells. That is, it appears that a relative amount of sparseness can be distinguished. In scenario B, the PYR cells are less sparse and are less tightly bound to fire phase-locked with PV ϩ cells due to the postinhibitory rebound, and the PV ϩ cell firing is more sparse. By contrast, in scenario A, the PV ϩ cells fire more, and postinhibitory rebound plays a stronger role to more tightly control PYR cell phase-locking (so more tightly lined up with PV ϩ cells). In this way, scenarios A and B can be differentiated by how much of a role postinhibitory rebound plays in the subsequent theta rhythm. These differences can be seen by comparing the cases shown in Fig. 10 and examining explicit numbers as shown in Table 4. We note that although it is clear that our models are in line with experimental observations in terms of PV ϩ cells spiking more than PYR cells during theta rhythms (population bursts) in the model (compare fifth and sixth columns of Table 4), how much spiking PYR and PV ϩ cells exhibit is not in line with experiment. Specifically, model PV ϩ cells do not spike on each population burst as has been observed experimentally.
Excitatory and inhibitory currents for the chosen parameter sets in Table 4 are shown in Table 5. Specifically, EPSCs and IPSCs to PYR and PV ϩ model cells are measured, and these values along with their ratios are given in Table 5. We note that there can be large EPSCs to PV ϩ cells and small EPSCs to PYR cells as observed experimentally. However, there is not always an appropriate match; IPSCs to PYR cells are too large and IPSCs to PV ϩ cells are too large in some cases. If current ratios rather than currents are compared, then it is always the case that EPSC/IPSC ratios are appropriate for PYR cells relative to experiments, but only in some cases are the EPSC/IPSC ratios for PV ϩ cells somewhat appropriate, that is, close to or greater than 1. As such, we find that the EPSC/IPSC ratio to PV ϩ cells, but not the EPSC/IPSC ratio to PYR cells, allows us to distinguish between scenarios A and B. We then conclude that scenario B, but not scenario A, is consistent with the experimental data, and so is the situation that is appropriate for the biological system; that is, one in which postinhibitory rebound, although required to be present, plays less of a role in theta rhythm generation.

Summary, theta essence, explanation, and predictions
We have developed microcircuit models and obtained an explanation for how theta rhythms can be generated in Figure 9. Changing excitatory conductance g pyr-pv from PYR to PV ϩ cells. Three example raster plots are shown, with the g pyr-pv value shown in red with the network frequency as appropriate. The average number of cells per population burst is 290 (PYR) and 343 (PV ϩ ) for g pyr-pv ϭ 0.5 nS; 81 (PYR) and 288 (PV ϩ ) for g pyr-pv ϭ 2 nS; 32 (PYR) and 238 (PV ϩ ) for g pyr-pv ϭ 5 nS. Other parameter values: g e,mean , e , g pyr , g pv-pyr , c PYR,PV , c PV,PYR ϭ [0, 0.2, 0.084, 6, 0.4, 0.5]. the hippocampus. We used a strategy, as schematized in Fig. 1, in which we took advantage of an experimental context of an intrinsic CA1 theta with developed mathematical models, leveraged theoretical studies, and did extensive parameter variation analyses. This computational analysis allowed us to differentiate between two scenarios of how theta rhythms could be generated, and only one of them is consistent with the experimental data.
We suggest that spike frequency adaptation and postinhibitory rebound in CA1 pyramidal cells are sufficient conditions (building blocks) for the generation of theta rhythms with sparse excitatory cell firing. Moreover, if it is the case that spike frequency adaptation is required, then it is necessary to have postinhibitory rebound. Further, our network simulations predict that theta rhythms are present when the input to the PYR cells has a nonzero   The same examples as in Table 4 but now showing excitatory and inhibitory currents and the delineated scenarios fluctuating input conductance of 0.6 nS or less and a zero mean conductance. Our network models are minimal but were able to capture an essence of the experimental data. As such, we consider our models a foundation on which to build. This building block understanding leads to the following claims and predictions. We first note that if PV ϩ cells are removed from the network, then there are no spontaneously generated theta rhythms. This is in accordance with Amilhon et al. (2015), in which the optogenetic silencing of PV ϩ cells abolished intrinsically generated theta rhythms in a whole hippocampus preparation in vitro. We predict that spike frequency adaptation in PYR cells is required for theta rhythms to be present, as balanced by a large enough PYR cell network with connectivity that is minimal. The amount of spike frequency adaptation controls the existence of theta rhythms and its resulting frequency, as built on our understanding from the theoretical mechanism. Thus, if cellular adaptation in PYR cells is selectively adjustedsay, by modulating potassium and calcium-activated potassium channels-then one should see an effect on the frequency of theta rhythms. If cellular adaptation is reduced enough, there would be no theta rhythms. Further, by selectively altering either the amount of connectivity or the conductance from PYR to PV ϩ cells, the theta frequency would be strongly affected. Alternatively, if the amount of connectivity or the conductance from PV ϩ to PYR cells is selectively adjusted, then the theta frequency would be less affected. Too much reduction, however, would result in no theta rhythms, as postinhibitory rebound would not be present. Finally, we explored a wide range of parameter values in our network models, but only a portion of them are in agreement with the experimental data regarding EPSC/IPSC ratios. As such, we predict that theta rhythms in the hippocampus are generated via scenario B, where the EPSC/IPSC ratio to PV ϩ cells is greater than 1 (Fig. 7). Although postinhibitory rebound is still an important element to have theta rhythms (i.e., c PV,PYR cannot be zero) in scenario B, it plays less of a role relative to scenario A. If scenario B is the mechanism underlying theta rhythms in the hippocampus, we predict that the probability of connections or conductance from PV ϩ to PYR cells is larger than the probability of connection or conductance from PYR to PV ϩ cells (see Fig. 7 schematic). It is of interest to note that large EPSCs onto PV ϩ cells can be present even if there is a minimally connected CA1 PYR cell network that sparsely fires. This is necessarily due to the large size of the PYR cell network.
Our explanation here of theta generation in the hippocampus has Aristotelian elements of efficient, material, and formal causes (Falcon, 2015). That is, we find that spike frequency adaptation and postinhibitory rebound in large, minimally connected networks (efficient cause) consisting of PYR and PV ϩ cells connected by inhibitory (GABA A ) and excitatory (AMPA) synapses (material cause) generates theta rhythms. The models require parameter balances of cellular adaptation in a large PYR cell network with enough connectivity and postinhibitory rebound due to PV ϩ cells, and larger PV ϩ to PYR than PYR to PV ϩ cell connectivity (formal cause). We note that a final cause cannot be considered with the present model, as it uses a whole hippocampus preparation in vitro, so function cannot be directly addressed. There are several ongoing efforts to consider theta rhythm functions (Colgin, 2013), but as emphasized by Colgin (2013), it is important to know what the underlying causes of theta generation are in the first place. In future work, for example, one might consider how our modeling work could be expanded and linked to functional theta modeling studies such as recent efforts by Chadwick et al. (2016) to capture the flexibility of theta sequences by including phase-precessing interneurons in septo-hippocampal circuitry.
It is important to emphasize that the form of theta modeled and explained here is one that is intrinsic to the hippocampus. However, with a baseline model of how the intrinsic CA1 microcircuit can produce theta rhythms in the absence of oscillatory synaptic input, we could gain insight into how external sources influence and build on this CA1 microcircuit. There is always the possibility that in some cases the microcircuit is just perturbed and the mechanisms remain intact, but in other cases it could be that the mechanisms are entirely different. Our models may help to identify which is the case, as it makes clear testable predictions as described above.

Related modeling studies
Previous mathematical modeling studies of theta rhythms have been described (Kopell et al., 2010;Ferguson and Skinner, 2015). Earlier models put forth theta generation mechanism ideas based on coherence between theta frequency firing in oriens-lacunosum/moleculare (O-LM) interneurons, in which a hyperpolarizationactivated inward current (h-current) is critical for the coherence (Rotstein et al., 2005;Wulff et al., 2009). However, this has subsequently been shown to be unlikely, since O-LM cells do not operate as theta pacemakers (Kispersky et al., 2012), an assumption in the earlier models. In another hippocampal modeling study, it was shown that theta rhythms could be generated in a network of basket, O-LM, and pyramidal cells (Neymotin et al., 2011). Although the focus of that study was not on theta generation per se, O-LM cells were strongly implicated in contributing to theta rhythms.
Based on recordings from the hippocampus of behaving rats, models were developed to support the observations that h-currents in pyramidal cells were needed to allow theta frequency spike resonance to occur (Stark et al., 2013). The mathematical models were focused on considering the contribution of h-currents in pyramidal cells and did not specifically consider connectivity between PYR and PV ϩ interneurons or the numbers of cells. The study emphasized the importance of postinhibitory rebound in theta rhythms. Interestingly, this same study found that adaptation was not an important contributing factor. In our modeling study, both cellular adaptation and postinhibitory rebound are needed to bring about theta rhythms in the hippocampus. However, unlike Stark et al. (2013), which had an experimental in vivo context, our models are based on an intrinsically generated hippocampal theta rhythm, thus removing contributions from other brain structures. In this way, we were able to focus on whether and how excitatory-inhibitory networks could produce theta rhythms considering network size, connectivity, and cellular characteristics. Given the complexity of theta rhythms, one would expect that there could be different balances of building blocks (cellular adaptation and postinhibitory rebound), as well as additional ones underlying theta rhythms in vivo. However, it is considerably more challenging to dissect out an understanding of theta rhythm generation in vivo. Using our network models and mechanistic understandings derived from them, this challenge could be reduced.
A full-scale biologically detailed CA1 hippocampal model has recently been developed (Bezaire and Soltesz, 2013;Bezaire et al., 2016). It is loosely based on the whole hippocampus preparation and includes eight inhibitory cell types. This model exhibits theta rhythms phaselocked with gamma oscillations and shows distinct phase relationships for different cell types. From model perturbations, interneuronal diversity, and more specifically, parvalbumin-expressing interneurons and neurogliaform cells, were found to be important in theta generation. Further, although not a particular focus, the generation of theta rhythms in these detailed network models required a particular balance of excitation. Considering Aristotle's four causes (Falcon, 2015), this descriptive understanding of theta generation can be considered as a material cause with elements of an efficient cause. However, due to its very detailed nature, it would be difficult to acquire a formal cause from it.
In a recent study Giovannini et al. (2017) focused on the contribution of a nonspecific cation current in pyramidal cells as critical for the maintenance of theta oscillations in the isolated hippocampus preparation. Interestingly, similar to our work here, they showed that when pyramidal cells were coupled with inhibitory cells, theta oscillations became more robust. This work could be considered as a particular material cause explanation of theta rhythm generation.

Spike frequency adaptation and postinhibitory rebound
Spike frequency adaptation as a mechanism to generate populations bursts has been used before, and previous work examined whether the amount of cellular adaptation expressed in pyramidal cells was appropriate to generate population bursts in pyramidal cell networks (Dur-E-Ahmad et al., 2011;Ferguson et al., 2015a). How best to model and examine adaptation naturally depends on the questions being considered. For example, Benda et al. (2010) examined how adaptation in the signal transmission of sensory systems can be due to either adaptation currents or dynamic thresholds. Coupled oscillator theory was used to examine how adaptation contributed to synchronization (Crook et al., 1998). In another study, a fractional leaky integrate-and-fire model to capture spike frequency adaptation was developed to set a framework to help understand information integration in neocortex (Teka et al., 2014). For our CA1 pyramidal cell models, we used an Izhikevich-type cellular model (Izhikevich, 2006), and adaptation was captured in the d and a parameters when fitted to experimental frequency-current curves (Ferguson et al., 2015b). Further, for postinhibitory rebound to be present in these models, the b parameter needed to be positive. Although the experimental data indicated that both strongly and weakly adapting pyramidal cells exist, this does not of course mean that there are simply two types, as the existence and amount of adaptation depends on the complement of biophysical ion channels in the cells and what is uncovered by the particular experimental protocol. In the work here, we used both types of pyramidal cell models, but networks of weakly adapting cells only were not able to produce theta rhythms given other constraints of the experimental context.
In the intact hippocampus in vitro preparation, CA1 pyramidal cells exhibit postinhibitory rebound inhibition Ferguson et al., 2015b). H-currents in pyramidal cells clearly play a key role in their ability to express postinhibitory rebound. However, both the presence and the distribution of these currents, together with the distribution of other currents, need to be taken into consideration. Specifically, it was shown that postinhibitory rebound is rarely observed in physiologic conditions unless unmasked by the blocking of A-type potassium currents (Ascoli et al., 2010), and rebound and other properties vary along the longitudinal axis of the hippocampus (Malik et al., 2016). Both h-currents and A-type potassium currents are known to have a nonuniform distribution along pyramidal cell dendritic arbors, and putative functional contributions of this to temporal synchrony have been made (Vaidya and Johnston, 2013). It is interesting to note that a difference in the dorsal to ventral patterning of h-currents exists (Dougherty et al., 2013), bringing to light another possibility of rhythm modulation. Further, traveling theta waves have been observed in both rodent (Lubenov and Siapas, 2009;Patel et al., 2012) and human (Zhang and Jacobs, 2015) hippocampus, suggesting a coupled oscillator organizational motif in the hippocampus. Although a "weakly coupled oscillator" terminology has been invoked in describing these waves (Colgin, 2013), this should not be confused with the mathematical theory where the assumption of weakly coupled oscillators is used to reduce the system to a phasecoupled system that is easier to analyze (Schwemmer and Lewis, 2012).

Limitations
Given our highly simplified and minimal network models, we did not expect to find a perfect match to the experimental data. However, it is important to note that given our minimal models, we were able to examine several thousand parameter sets, which in turn enabled us to explore and understand what balances might be important in bringing about theta rhythms. This balance and building block understanding can serve as a basis for how theta rhythm frequency and existence can be modulated by additional inputs from other brain regions, as well as modulation that would affect adaptation and postinhibitory rebound.
Our network models are minimal, but they were able to produce theta rhythms with sparse firing, as represented by population bursts, allowing us to suggest sufficient and necessary conditions for their generation. Although our models took into consideration network size, connectivity, and cellular characteristics in a clear experimental context, architecture was not considered. That is, connectivity used in the network models was random. This is clearly a simplification, especially considering the recent finding of motifs in pyramidal cells in the CA3 region of hippocampus (Guzman et al., 2016). However, it is a reasonable first approximation that allowed us to explore a wide expanse of connectivities.
Variability in intrinsic cell properties and mixing of weakly and strongly adapting pyramidal cells could be considered. Introducing variability immediately raises the question of how it should be done. For example, Harrison et al. (2015) quantified the heterogeneity within and between the neocortical pyramidal cell classes. They examined the class-dependent variance and covariance of electrophysiological parameters using simple models but did not focus directly on network dynamics or rhythms. Given our large explorations of conductances, connectivities, and noisy input, some variability was present in our simulations.
Further, only one type of inhibitory cell was included in our networks, that is, the fast-firing PV ϩ cell type. It is unlikely that only this one type of inhibitory cell contributes to theta rhythms, but we identified this as a good place to start given that Amilhon et al. (2015) found that they were essential. By no means do our models imply that other inhibitory cell types are unimportant. On the contrary, since our models clearly do not fully capture the experimental data (e.g., model PV ϩ cells fire too sparsely relative to experiment), aspects are clearly missing. Because PV ϩ fast-firing cells also include bistratified and axo-axonic cells, an expansion of PV ϩ cell networks along these lines could be considered (Ferguson et al., 2015c). Given the diversity of inhibitory cell types (Chamberland and Topolnik, 2012) and the different types of PV ϩ cells (Baude et al., 2007), we did not specifically try to scale the number of PV ϩ cells as we did for PYR cells; however, since we fully explored connectivity ranges between PV ϩ and PYR cells, this was in effect included. The inclusion of O-LM cells and other inhibitory cell types in the network models is important moving forward to be able to understand how they modulate theta rhythms (Leao et al., 2012;Sekulic and Skinner, 2017).

Conclusions and Future Work
Our network models represent closed, self-consistent, accessible models that can generate theta rhythms in hippocampus CA1. We intend them as a starting point on which to build to understand the mechanisms underlying theta rhythm generation in the hippocampus. Specifically, the network model building block balances need to be fully analyzed so that a more solid formal cause of explanation can be obtained, beyond what was obtained from the computational, parameter variation analyses done here. Further, in combination with full-scale models such as that of Bezaire et al. (2016), it may be possible to obtain an understanding that can fully encompass efficient, material, and formal causes in the Aristotelian sense and could subsequently help in a final cause understanding.
Overall, our models can serve as a backbone on which other cell types as well as details of particular cell types (biophysical channels, dendrites, and spatial considerations), modulatory effects, and input from the medial septum can be incorporated. However, in doing this, it is important to note that interaction and testing with experiment should be designed accordingly, given the strategy used in developing our models (Fig. 1), and that a consideration of the different forms of theta is not lost. Moving forward, it will be important to include biophysical LFP models to allow direct comparisons between models and experiments regarding LFPs and provide further constraints. Interestingly, it has been shown that active subthreshold currents can lead to distinct resonances in the generation of LFPs (Ness et al., 2016).