A Hypothesis for Theta Rhythm Frequency Control in CA1 Microcircuits

Computational models of neural circuits with varying levels of biophysical detail have been generated in pursuit of an underlying mechanism explaining the ubiquitous hippocampal theta rhythm. However, within the theta rhythm are at least two types with distinct frequencies associated with different behavioral states, an aspect that must be considered in pursuit of these mechanistic explanations. Here, using our previously developed excitatory-inhibitory network models that generate theta rhythms, we investigate the robustness of theta generation to intrinsic neuronal variability by building a database of heterogeneous excitatory cells and implementing them in our microcircuit model. We specifically investigate the impact of three key “building block” features of the excitatory cell model that underlie our model design: these cells' rheobase, their capacity for post-inhibitory rebound, and their spike-frequency adaptation. We show that theta rhythms at various frequencies can arise dependent upon the combination of these building block features, and we find that the speed of these oscillations are dependent upon the excitatory cells' response to inhibitory drive, as encapsulated by their phase response curves. Taken together, these findings support a hypothesis for theta frequency control that includes two aspects: (i) an internal mechanism that stems from the building block features of excitatory cell dynamics; (ii) an external mechanism that we describe as “inhibition-based tuning” of excitatory cell firing. We propose that these mechanisms control theta rhythm frequencies and underlie their robustness.

1. Introduction

Hippocampal theta rhythms (≈3–12 Hz) observed in local field potential (LFP) recordings are associated with cognitive processes of memory formation and spatial navigation (Colgin, 2013, 2016; Hinman et al., 2018). Exactly how theta rhythms emerge is a complicated and multi-layered problem, but it is known that there are two types of theta in the hippocampus, denoted type 1 and type 2, that have high (7–12 Hz) or low (4–7 Hz) frequencies, respectively. Type 2, but not type 1, rhythms are dependent on cholinergic drive (Kramis et al., 1975; Bland, 1986; Buzsáki, 2002). That is, type 2 rhythms are atropine-sensitive and type 1 are not. In rodents, it has been shown that social stimuli elicit high theta, and fearful stimuli elicit low theta (Tendler and Wagner, 2015), and type 2 theta oscillations have been shown to be associated with increased risk-taking behavior (Mikulovic et al., 2018). In humans, theta frequencies are lower overall (Jacobs, 2014), but it is still possible to distinguish high and low theta frequencies, with low theta supporting encoding and retrieval of memories (Kota et al., 2020). Clearly, what the particular theta frequency is has functional significance, and so how theta frequencies are controlled is functionally important to consider.

It is now well-documented that theta rhythms can be generated intra-hippocampally, emerging spontaneously from an isolated whole hippocampus preparation in vitro (Goutagny et al., 2009). In other words, hippocampal microcircuitry is able to produce theta rhythms on its own. Simultaneous access to cellular and population output in this in vitro preparation presented an opportunity to untangle cellular and population dynamics of how theta rhythms are generated. We took advantage of this and built cellular and microcircuit models that could generate theta rhythms with parameters directly constrained by experimental data from the whole hippocampus preparation and the experimental literature (Ferguson et al., 2013, 2015a, 2017). Motivated by the perspective presented by Gjorgjieva et al. (2016), we considered a “building blocks for circuit dynamics” analysis approach in our microcircuit model design (Ferguson et al., 2017). In this perspective, biologically known cellular, synaptic, and connectivity characteristics are considered as building blocks for circuit dynamics. For example, one such cellular “building block” is post-inhibitory rebound (PIR), which has previously been invoked as a contributor to the generation of cortical oscillations (McCormick et al., 2015).

In this paper we use our theta-generating microcircuit model to develop a hypothesis of how the theta frequencies could be controlled. We first describe the model microcircuit design and then assess the robustness of theta generation in the model by considering heterogeneous pyramidal (PYR) cell populations. From this, we use phase response curves (PRCs) and show that inhibitory inputs affect the theta frequency. We thus propose a hypothesis for theta frequency control in CA1 microcircuits that is dependent on internal features of PYR cells and “inhibition-based tuning” of PYR cell firing. We summarize our study in schematic form in Figure 1.

FIGURE 1

Figure 1. Schematic showing aspects involved in the hypothesis developed in this study. Theta rhythms are generated intrinsically in a whole hippocampus preparation of Goutagny et al. (2009) (“Experiment”). Their generation is captured in a microcircuit model design by Ferguson et al. (2017) (“Model Networks”). In the present paper we assess the robustness of this model design and develop a hypothesis for theta frequency control (“Hypothesis Development”).

2. A Design of Microcircuit Models That Produce Theta Rhythms

We have built cellular-based excitatory-inhibitory (E-I) network models (Ferguson et al., 2017) to understand how the theta rhythms observed in a whole hippocampus preparation by Goutagny et al. (2009) could be generated. The model networks (see Figure 1 schematic) are designed to represent a “piece” of the CA1 region of the hippocampus – ~1 mm³ that was determined to be enough to self-generate theta rhythms (Goutagny et al., 2009; Ferguson et al., 2017). It includes only two distinct cell types, pyramidal (PYR) cells and fast-firing parvalbumin-positive (PV+) cells, as represented by a single compartment model with an Izhikevich mathematical model structure (Izhikevich, 2006). These cellular models have parameters determined from fits to electrophysiological recordings from the whole hippocampus preparation to match frequency-current (f-I) curves of PYR cells (Ferguson et al., 2015a) and fast-firing PV+ cells (Ferguson et al., 2013). The individual PYR cell model does not generate bursting output. The model network consists of 10,500 cells (10,000 PYR cells and 500 fast-firing PV+ cells) (Ferguson et al., 2013, 2015b), with connectivity parameters estimated from the experimental literature. In particular, there is sparse coupling (<1%) between PYR cells. We note that we took advantage of a scaling relationship between cell number, connection probability and excitatory synaptic weight that allowed us to use 10,000 PYR cells rather than the 30,000 cell number size as estimated for the “piece” of tissue.

We examined our models from a “building block for circuit dynamics” perspective (Gjorgjieva et al., 2016) to determine if theta rhythms (i.e., theta frequency population bursts) could be generated according to experimental constraints. We first found that experimentally constrained PYR cell network models (E-cell networks alone) could generate population bursts of theta frequency (Ferguson et al., 2015b), suggesting that a cellular “building block” feature of spike frequency adaptation (SFA) present in the constrained PYR cell models could be an important contributor to theta rhythm generation. A sample of output from a E-cell only network is shown in the Supplementary Material. However, we also found that in these E-cell only networks the PYR cells do not fire sparsely as was observed experimentally (Huh et al., 2016). When we included PV+ cells to create E-I model networks, population bursts of theta frequency were still possible and were now associated with sparse PYR cell firing in accordance with the experimental data. As the addition of PV+ cells allows PIR to be possible in the PYR cells, we consider PIR as another building block feature of importance in generating these intrinsic theta rhythms. Along with SFA and PIR features, the PYR neurons have an inherent rheobase (Rheo) feature, which is the amount of current required to make the PYR cell spike (derived from fitting to the experimental data in Ferguson et al., 2015a). We consider this to be a third building block feature for theta rhythm generation. Further, for the model output to be consistent with experimental observations of excitatory postsynaptic current (EPSC) and inhibitory postsynaptic current (IPSC) amplitude ratios, we found that the connection probability from PV+ to PYR cells was required to be larger than from PYR to PV+ cells—a particular prediction that has been examined and found to be consistent with empirically derived connectivities (Chatzikalymniou et al., 2020).

3. An Assessment of the Model Design for Robust Theta Rhythms

In our previous work, we did not specifically examine the sensitivity of theta rhythms to SFA, PIR, or Rheo features. To address this here, we create a model database of 10,000 PYR cell models.

3.1. PYR Cell Model Database and SFA/PIR/Rheo Features

While there are various ways in which a model database could be created, we do this by simply varying specific parameter values of the PYR cell model (Equation 1) in a regular fashion.

$\begin{array}{l}{C}_m\stackrel{\bullet }{V}=k(V-v_r)(V-v_t)-u+I_{other}\\ \stackrel{\bullet }{u}=a\left[b(V-v_r)-u\right]\\ \\ \text{if }V\ge v_{peak},\text{then }V\leftarrow c,u\leftarrow u+d\\ \text{where }k=k_{low}\text{if }V\le v_t,k=k_{high}\text{if }V>v_t& (1)\end{array}$

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 cut-off value, a (ms⁻¹) 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. I_other is applied to determine values for PYR cell features, as described below. Model parameter values (units above) for the PYR cell are: v_r = −61.8; v_t = −57; v_peak = 22.6; c = −65.8; k_high = 3.3; C_m = 115; a = 0.0012; b = 3; d = 10; k_low = 0.1. These parameters are as previously determined for strongly adapting cells (Ferguson et al., 2015a). We refer to them as default parameter values, and specifically, the a, b, d, k_low parameters are varied in creating the model database. Further details about the model database are provided in the Supplementary Material.

From the created model database of PYR cell models, we obtain a range of SFA, PIR, and Rheo features that are quantified in the following fashion. For Rheo: Starting from v_r, each PYR cell model is given a constant current from −25 to 25 pA in 0.5 pA increments. If a spike is generated within the first 500 ms, then that constant current value is considered as the Rheo quantified value. For PIR: Starting from v_r, each PYR cell model is subjected to a one second hyperpolarizing step current for current values from 0 to −25 pA with a resolution of 0.5 pA. If a spike occurs upon termination of a given hyperpolarization step (i.e., a PIR spike) but not at the previous step value, then that step value is considered as the PIR quantified value. For SFA: Starting from v_r, each PYR cell model is subjected to input currents for one second, from 0 to 98 pA (inclusive) in 2 pA increments. For each input current, the number of spikes is recorded, and the interspike interval is calculated between the first and second spikes, and the last and second from last spike. The inverse is taken and defined as the initial and final frequency at that current. The initial and final frequencies as a function of the current steps creates a smooth, approximately linear relationship, so lines are fitted to the initial and final frequency plots. The slopes of those lines are subtracted from one another to produce the SFA quantified value. Therefore, the larger the quantified SFA value is, the stronger is the amount of the PYR cell adaptation, i.e., we get more reduction in the PYR cell spike frequency for a fixed amount of input current; the more negative the quantified PIR value is, the larger is the hyperpolarizing step required to generate a spike at the end of the step; the larger the quantified Rheo value is, the more input is required to cause the cell to spike. For the PYR cell model with default parameter values, the quantified values for the building block features are: SFA = 0.46 Hz/pA, Rheo = 4.0 pA, and PIR = −5.0 pA. We refer to these as base feature values. With the created database of PYR cell models, we obtain a range of building block feature values distributed as shown in Figure 2. Further details about SFA, PIR, and Rheo features are provided in the Supplementary Material.

FIGURE 2

Figure 2. Distributions of PYR cell features from created model database. A heterogeneous set of PYR cells was created and their “building block” features of SFA, Rheo, and PIR were quantified. Details of this quantification are provided in the Supplementary Material. Histograms show the number of occurrences of SFA [=] Hz/pA, Rheo [=] pA, PIR [=] pA values, and vertical black arrows indicate [SFA,Rheo,PIR] base values. Also shown are narrow (N) and broad (B) subsets of heterogeneous PYR cell populations and low (L), medium (M), or high (H) subsets of heterogeneous PYR cell populations that do or do not include base building block values. SFA histogram has a bin resolution of 0.05, and Rheo, PIR histograms have a bin resolution of 0.5.

3.2. Observations

In the extensive E-I network simulations of Ferguson et al. (2017), the PYR cell models used were homogeneous, and all had default model parameter values. However, the networks themselves were not homogeneous because of the noisy external drives to the PYR cell models. To examine the robustness of the theta-generating mechanism in the E-I network models to variability in the SFA, PIR, and Rheo features, we create heterogeneous PYR cell populations from the model database and examine whether the presence of theta rhythms in E-I networks is affected by varying these building block features. We carry out our examination such that the heterogeneous PYR cell population in the E-I networks either does or does not include PYR cells that have base values. As a brief aside, we note that when we examine E-I networks that have homogeneous PYR cell models with parameter values different from the default ones, but that have similar SFA, PIR, and Rheo base values, the resulting networks produce clear population bursts, but with a bit of variation in frequency and power. Specific examples are provided in the Supplementary Material.

For E-I networks with heterogeneous PYR cell populations that have PYR cells that do include SFA, Rheo, and PIR base values, theta rhythms continue to be expressed. We also find that the network theta power is larger when there is a narrow rather than a broad range of values encompassing base ones. Figure 2 shows the narrow and broad ranges of values in our created database. Further details are provided in the Supplementary Material. This observation of theta power difference suggests that particular quantified feature values affect the robustness of theta rhythms since the power is larger when it more narrowly encompasses base values.

For heterogeneous E-I networks that have PYR cells that do not include base values for all features, we build E-I networks that have a low (L), medium (M), or high (H) range of values for SFA, Rheo, and PIR features in different combinations. Thus, a given heterogeneous E-I network has a triplet of [SFA,Rheo,PIR] features that have a L, M, or H range of values. These values are shown in Figure 2. In Figure 3, we show the frequency (left) and power (right) of the output of these heterogeneous E-I networks designated by dots of a given color. The red circled dot is the only E-I network that does have base values for all of the building block features, i.e., [SFA,Rheo,PIR] = HML. We observe the following for the network frequency: Networks with Rheo = L do not produce theta rhythms when PIR and SFA = M or H; There are no theta rhythms when Rheo = M values and SFA and PIR = H; As Rheo increases, the network frequency increases, and there appears to be a stronger control of frequency by the Rheo feature relative to SFA and PIR features. For the theta power, we find that it is lowest when Rheo = L and increases as Rheo increases, but decreases as SFA or PIR increase. However, when Rheo = M, the power increases as SFA increases and as PIR decreases. From these trends, it would appear that the Rheo feature controls the theta frequency and power more than SFA or PIR. As larger values of Rheo refer to larger depolarizing currents being required for the PYR cell to fire, our observations imply that the amount of current needed for a PYR cell to fire is an essential controller of theta frequency and power, assuming that other features allow rhythms to exist in the first place. Further details from this examination are provided in the Supplementary Material.

FIGURE 3

Figure 3. Frequency and power of theta rhythms in heterogeneous E-I networks. Each dot represents the frequency (left) or power (right) of the output of the network that has [SFA,Rheo,PIR] features with a L, M, or H range of values as plotted, with the dot color representing the specific frequency or power value given in the color bar. The red circled dot is the network that has feature values that include base values for all of the features, i.e., [SFA,Rheo,PIR] = HML. The dark blue circles do not produce a rhythmic output, and the vertices that do not have any dots are where there were no individual PYR cell models to generate the particular heterogeneous network. Further details are provided in the Supplementary Material.

In summary, the exploration of our microcircuit model of theta rhythm generation in the whole hippocampus preparation leads us to the following conclusions regarding the influence of the three “building blocks” on this dynamic: (i) a larger theta power occurs in E-I networks with heterogeneous PYR cells that include their base values and are narrowly distributed around them, and (ii) particular rheobase current values control the frequency and power of network rhythms more than the ability of the PYR cell to spike on inhibitory rebound or the particular amount of spike frequency adaptation. Thus, these simulations of E-I networks with heterogeneous PYR cell populations have allowed us to gauge the contributions of the different features and have helped us to confirm the robustness to cellular heterogeneity of the theta-generating rhythm mechanism in our microcircuit model design.

4. Using the Assessment and Design to Develop a Hypothesis for Theta Frequency Control

As described above, we find that large, minimally connected recurrent networks with fast-firing PV+ cells and PYR cells can produce theta frequency population rhythms consistent with experiment, driven, and controlled in part by the building block features of SFA, PIR, and Rheo in PYR cells. In our previous I-cell only network models of PV+ cells, coherent network output was possible with experimentally constrained PV+ cellular models and synaptic connectivities (Ferguson et al., 2013). A sample of output from an I-cell only network is shown in the Supplementary Material. In creating the E-I network model setup, the PV+ cell network was “designed” to be in a coherent state—a function of the appropriate excitatory drive being received and the connectivity of PV+ cells. Specifically, we chose the synaptic weight (between PV+ cells) to be such that it could be at the “edge” of firing coherently (high frequency) or not (see Figure 3 in Ferguson et al., 2013), and as such, given an appropriate excitatory drive from the PYR cells, the PV+ cell network could be in a high frequency coherent regime and be considered to be producing an inhibitory “bolus” to the PYR cells. This is an important consideration for our phase response curve (PRC) considerations below.

From the several model sets of heterogeneous E-I model network outputs described in the previous section, we choose three that exhibit strong population rhythms of different frequencies. Details on these three chosen networks (specifically the heterogeneous PYR population as well as the classification of their rhythms as “strong”) can be found in the Supplementary Material. Raster plot outputs of the PYR cells in these chosen heterogeneous E-I networks are shown in Figure 4 where the different rhythms are referred to as “slow,” “medium,” and “fast.” Given the minimal nature of the microcircuit model, the frequencies of these rhythms fall a bit outside theta ranges (higher) for some networks, although the underlying theta generation mechanism and the model design is the same.

FIGURE 4

Figure 4. Raster plot outputs of PYR cells from three heterogeneous E-I simulations. These three model sets generating population burst rhythmic output exhibit three different frequencies that we refer to as “slow” (9.6 Hz), “medium” (13 Hz), and “fast” (15 Hz) from their respective model sets. For all three sets, the heterogeneous PYR cells include those with Rheo base values, whereas only the model set producing the “medium” output has PYR cells with SFA base values. Except for the model set producing “slow” output, PYR cells have PIR base values. That is, the triplet [SFA,Rheo,PIR] feature for the slow, medium and fast networks are MMH, HML, and LML, respectively.

Let us now take advantage of our microcircuit design to examine how these frequencies are controlled by turning to PRC considerations (Schultheiss et al., 2011). We note that PRCs are commonly calculated using a brief, strong, excitatory current pulse as a perturbation. We slightly modify that paradigm here and instead use a negative pulse whose amplitude and duration is motivated by the type of synaptic inputs generated during an “inhibitory bolus” in our network model (see Figure 5). We know that the PYR cell network can generate theta population bursts on its own given its cellular adaptation characteristics (SFA feature; Ferguson et al., 2015b). While on their own the PYR cells do not fire sparsely as in experiment, they do when a PV+ cell population is included (Ferguson et al., 2017). We consider that the resulting frequency of the E-I network's population bursts is due to a combination of the individual PYR cell's firing frequency and how much an inhibitory input could advance or delay the PYR cell spiking (as quantified by PRCs). The setup to consider this is schematized in Figure 5 and consists of the following: Each PYR cell in the heterogeneous population receives excitatory input from other PYR cells as well as a noisy drive (other input). The amount of input a PYR cell receives would of course fluctuate over time, but under reasonable approximation the PYR cell receives a mean excitatory input of about 20–30 pA. This approximation is based on the fact that in our E-I network models (see Figure 1), theta population bursts occur when PYR cells receive a zero mean excitatory drive with fluctuations of ≈10–30 pA (Ferguson et al., 2017). We then calculate PRCs as described above. The inhibitory pulse can advance or delay the subsequent PYR cell's spike as quantified by the PRC, which in turn is dependent on the PYR cell's intrinsic properties. All of these aspects are schematized in Figure 5.

FIGURE 5

Figure 5. Motivation for the setup in phase response curve (PRC) calculations. Assuming a theta-generating mechanism based on model design, PRCs are generated based on an inhibitory input (“bolus”) coming from the PV+ cell network to an individual PYR cell. Each PYR cell is receiving a noisy drive shown as “Other Input” as well as from other PYR cells. An illustrative f-I curve is shown for an individual PYR cell. An illustration of a computed PRC based on the inhibitory input to a particular PYR cell is also shown. It would be dependent on the particular PYR cell's model parameter values that dictates its f-I curve.

We consider the three cases of heterogeneous E-I networks exhibiting different population burst frequencies shown in Figure 4 and described as having a “slow,” “medium,” or “fast” population burst frequency output. We generate PRCs for the several PYR cell models in the population for each of these model sets that produce the different frequency population burst outputs. Each PYR cell model in the heterogeneous population has particular PRC characteristics due to its given model parameter values, and thus exhibits a specific intrinsic frequency for a given input.

4.1. PRC Calculations

These proceed as follows: A set input current (20:2:30 pA) is tonically applied to the model cell, and the period (defined λ) of the cell's firing is calculated as the time between the 9th and 10th cell spike. The inverse of the period represents the firing frequency of the cell, reported as averages and standard deviations for entire model sets. We compute the phase response of a model neuron to a perturbation at 100 equidistant times in its normal firing cycle, where the perturbation is a 1 ms current pulse with −500 pA amplitude (as mentioned previously, considered an approximation of the synaptic input received by these cells following an “inhibitory bolus”). For 1 ≤ i ≤ 100, we define $\Delta p=\frac{\lambda }{100}$ and deliver the perturbation at i*Δp ms after the 10th cell spike. We then measure the time between the 10th and 11th cell spike as the “perturbed period” (defined λ_p). We calculate the difference between this and the previously calculated period (in the absence of any perturbation) and normalize this by the normal firing period, meaning that in the PRC plots the y-axis is $\frac{\lambda -{\lambda }_p}{\lambda }$. This means that negative values plotted in the PRC correspond with a phase-delay, i.e., the perturbed period was longer than the unperturbed period, and vice-versa. The x-axis in the PRC plots are the normalized time at which the perturbation was delivered, simply calculated as $\frac{i}{100}$. We note that we perform this calculation separately for each i, i.e., we re-initialize the cell and let it respond naturally to a tonic input until the 10th spike for each value of i, rather than perform these perturbations sequentially and risk confounding the responses.

In Figures 6B,C, we quantify aspects of the PRC curves. In Figure 6B, we simply extract the value of the normalized phase difference from the mean PRC curve for a perturbation delivered at a normalized phase of 0.3 (denoted by the arrows overlaid on Figure 6A). In Figure 6C, we quantify one aspect of the mean PRC curve's rate of change, specifically the variability of the difference quotient calculated at each phase step, in the following straightforward way: first, this difference quotient is calculated for all but the last value of the normalized phase; second, the variance of these data is calculated simply using the var function in MATLAB.

FIGURE 6

Figure 6. Theta rhythm frequency is influenced by inhibitory drive as quantified via PRCs and firing frequencies of individual PYR cells. (A) Mean PRC (solid line) for the heterogeneous PYR cell population involved in “slow” (left), “medium” (middle), and “fast” (right) theta oscillations, calculated with an input current of 20 pA and with the shading representing ± the standard deviation. There are 25, 556, and 74 different PYR cell models in the 10,000 PYR cell populations of slow, medium, and fast cases, respectively. More details are provided in the Supplementary Material. The mean and standard deviation of the firing frequencies of the PYR cells at this input level are included in the inset of each panel. (B,C) After calculating both the mean PRC and mean intrinsic firing frequency for the PYR cell populations associated with our “slow” (red), “medium” (blue), and “fast” (green) theta oscillations for six input currents (20:2:30 pA), we extract a particular feature of the mean PRC (the mean phase shift caused by a perturbation delivered at a phase of 0.3 in (B) and the variance of the mean PRC's derivative in (C)) and plot it against the mean intrinsic firing frequency. In neither case is there a linear relationship between either axis and the theta rhythm frequency, indicating that it is a more complex combination that determines the population frequency. Note that, given the monotonic relationship between the input current and firing frequency in this range, the leftmost point for each color represents an input current of 20 pA, with each subsequent point moving rightwards representing the next input current step.

The code for generating and plotting these PRCs can be found at https://github.com/sbrich/Theta_PRCs. PRCs for input currents other than 20 pA that is shown in Figure 6A can be found at https://osf.io/yrkfv/.

4.2. Observations

In Figure 6, we first show an example of PRCs calculated for an input current of 20 pA (Figure 6A). PRCs are calculated for each model in a particular model set of heterogeneous PYR cell models, with the averaged curve presented along with a range of ± one standard deviation (shown by the shading around the curve in each plot of Figure 6A). These PRCs showcase distinct features: for instance, the PYR cells in the medium case uniquely exhibit a region of phase-advance, while the PYR cells in the fast case have the largest phase delay for perturbations delivered at all but the latest phases. Clear distinctions between the PRCs for each model set persist for all the input currents used.

To better visualize the influence of the intrinsic properties of the PYR neurons on theta rhythm frequency, we plot an extracted feature of the mean PRC against the mean firing frequency of these model sets for each of our computed input currents in Figures 6B,C, with the corresponding theta rhythm frequencies associated with each model set denoted by the data point's color, with the extracted PRC features in each case described in the previous section. These visualizations clearly illustrate that both the PRC and the mean intrinsic firing frequency of the PYR neurons in a given model set contribute to the overall theta rhythm frequency; otherwise, these points would be “flat” with respect to either the x or y axis. Furthermore, the relationship between the extracted PRC feature of interest and the mean intrinsic firing frequency varies notably depending on the output theta rhythm frequency: for instance, in Figure 6B both the “slow” and “medium” model sets show a monotonically decreasing relationship between the extracted PRC value and the mean intrinsic firing frequency, while the “fast” model set shows a monotonically increasing relationship. Taken together, these results show that it is a combination of the inhibitory drive and the PYR cell's excitability that contributes to the overall theta rhythm frequency.

The PRCs measure of how the PYR cells respond to perturbation help articulate potential mechanisms by which these differing theta rhythm frequencies arise. For instance, while the PYR cells in the fast case have the fastest individual firing frequencies (notably faster than what is seen in population models), their PRCs may be illustrative of how the inhibitory “bolus” decreases this firing frequency toward the theta range. Meanwhile, the PYR cells in the medium case have the slowest individual firing frequencies, although they participate in “medium” theta rhythm frequencies. The PRC in this case, particularly the region of phase-advance, may elucidate how inhibitory synaptic input actually accelerates PYR cell activity. These particular examples rely upon the PRC feature extracted and plotted in Figure 6B.

Thus, in the context of having stable population burst output given by an E-I model network design with appropriate SFA, PIR, and Rheo features, the PRC calculations here show that an appropriate inhibitory input contributes to the resulting population burst frequency. In essence, this analysis of the PRC features of our model sets supports a hypothesis that intrinsic properties of PYR cells, including those encapsulated by the PRC measure, play a critical role in controlling the frequency of population theta rhythms.

5. Discussion

We have used a microcircuit model, as designed to generate theta rhythms representing those observed in a whole hippocampus preparation, to develop a hypothesis for theta frequency control. Our work has allowed us to propose a hypothesis that encompasses two aspects: (i) an internal mechanism that stems from SFA, PIR, and Rheo building block features of PYR cells, and (ii) an external mechanism that involves an “inhibition-based tuning” of PYR cell firing. From previous work, we already knew that our sparsely coupled CA1 PYR cell E-only networks could produce theta frequency population bursts on their own, but the majority of the PYR cells would be firing, which is unlike the experimental observations of sparse PYR cell firing. With the inclusion of PV+ cells to create E-I networks, the population of PYR cells fired sparsely in accordance with experiment. It makes sense that the addition of inhibitory cells leads to less firing of PYR cells due to potential silencing from the inhibition. That theta rhythms of strong power can still emerge despite the participation of fewer PYR cells in the rhythm is likely due to the PV+ cells tuning the otherwise diverse frequencies of the PYR cells to similar frequencies, enabling this smaller group of cells to produce strong rhythms. This constitutes a main part of our proposed hypothesis. Relatedly, it has been shown that feedforward inhibition plays a role in maintaining low levels of correlated variability of spiking activity (Middleton et al., 2012).

It is important to highlight two key points that underlie our proposed hypothesis. First, the PYR cell population needs to be large enough so that it can collectively generate a strong excitatory drive to the inhibitory PV+ cells, and in turn the PV+ cell population should be able to fire enough (and coherently) to create a strong inhibitory “bolus” to tune the PYR cell population output. Second, the net input (recurrent excitation, excitatory drive, incoming inhibition) thus received by the PYR cells leads to the generation of theta rhythms and its resultant frequency. Some similarities exist between these points and the “PING mechanism” underlying the generation of gamma rhythms in E-I networks (Kopell et al., 2010; ter Wal and Tiesinga, 2013). In PING, the interneurons are recruited by high firing PYR cells which feed back to the PYR cells that subsequently fire when inhibition has sufficiently decayed, and it is interesting to note that recent research shows that rhythms with frequencies approaching the theta range can arise in PING-motivated networks (Rich et al., 2017). For the theta rhythms here, the large population of PYR cells is able to create a large enough excitatory drive to the inhibitory population that is able to fire coherently and feedback to tune the PYR cell firing to robust theta frequency rhythms. It is important to acknowledge that this mechanism and the subsequent frequency control hypothesis emerged from using E-I networks that have parameter values with direct links to the whole hippocampus preparation. That is, although the model is minimal, it did not have arbitrary model parameter values.

We do not know whether a clear relationship between PYR cell inputs and network frequency as described in the second key point above actually exists, and it would be highly challenging to directly examine this experimentally. We further note that the minimal nature of our E-I model networks regarding different inhibitory cell types and biophysical characteristics limits experimental design considerations. However, since a detailed, biophysical CA1 network model that includes PYR cells and eight different inhibitory cell types has been created (Bezaire et al., 2016), one can consider using it as a proxy for the actual biological system to start to explore this. We have done this by bringing together the described microcircuit model used herein and the detailed, full-scale CA1 microcircuit model, and examined how the theta network frequency produced by the detailed model depends on the net input received by the PYR cells (Chatzikalymniou et al., 2020). We found that the biologically detailed models strongly support this dependence and thus our proposed hypothesis for theta rhythm frequency control. Thus, this indicates that theta frequencies in the biological system may be controlled in such a fashion. In turn, this suggests that directly modulating intrinsic properties of PYR cells could strongly affect population theta rhythms.

In the previous work of Ferguson et al. (2015a), we had created PYR cell models that were either strongly adapting based on fits to the experimental data, or weakly adapting based on another experimental dataset. In Ferguson et al. (2015b), when either PYR cell models were used in E-cell only networks, that could produce theta frequency population bursts. As discussed in Ferguson et al. (2015a), it is unlikely that there are distinct types of biological PYR cells that are strongly or weakly adapting, but rather a continuum of adaptation amount dependent on the underlying balances of biophysical ion channel currents. Our explorations of the robustness of the theta generation mechanism in the microcircuit model here revealed that the frequency and power of theta rhythms were not strongly controlled by SFA feature values relative to Rheo feature values. Thus, although we created the model database starting from the strongly adapting PYR cell model parameter basis, it likely would not have mattered if the robustness examination of theta rhythm generation had been undertaken using weakly adapting PYR cell models instead.

It is perhaps not surprising that Rheo feature values are the main controller of the existence of theta rhythms and their frequency and power, as the particular Rheo value dictates whether a PYR cell would spike or not. We note that the experimental findings of Goutagny et al. (2009) had already suggested the importance of PIR in the generation of theta rhythms. In actual CA1 PYR cells, it has been shown that PIR spiking does occur, mediated by h-channels, and is locally controlled by biophysical ion channel balances (Ascoli et al., 2010). Whether PYR cells actually fire due to PIR during ongoing theta rhythms may or may not be the case, and one could potentially disentangle this in the models. However, the hypothesis developed in this work points to a confluence of features that culminate in the net current to individual PYR cells being a focus of theta rhythm frequency control. Thus, changes in the net drive to PYR cells or changes to the PYR cell's intrinsic properties such as h-currents that would affect PIR would be expected to affect the resulting theta rhythm frequency.

Previously, in vitro slices from hippocampus have been shown to generate theta rhythms (Fellous and Sejnowski, 2000; Gillies et al., 2002), and modeling studies have been built on these data to examine how theta rhythms are generated (e.g., see Gloveli et al., 2005; Orbán et al., 2006). The difference between such studies and the work here is that the in vitro whole hippocampus preparation spontaneously produces theta rhythms without any pharmacological or stimulatory influences, indicating that the theta rhythm generating mechanism in the whole hippocampus preparation is likely relevant in the behaving animal. Indeed, direct correlates between the in vitro preparation and in vivo recordings were shown regarding a reversal of theta rhythms (Jackson et al., 2014). Our E-I network models were built in direct correspondence with the whole hippocampus preparation, including the large PYR cell network size, and the specific cellular model parameters derived from fits to recordings from the whole hippocampus preparation. In this way, they are fundamentally different from the several models of theta rhythms that have previously been developed (Kopell et al., 2010; Ferguson and Skinner, 2018).

PRC theory has been used in a variety of ways in the Neuroscience field (Schultheiss et al., 2011), and particularly in consideration of network dynamics. For example, Hansel et al. (1995) used PRCs to explain the differential capacity for excitatory drive to synchronize networks of Type I or Type II neurons (these types are differentiated by their bifurcation type; Izhikevich, 2006), Rich et al. (2016) analyzed synchronization features in purely inhibitory networks using PRCs, and Achuthan and Canavier (2009) used PRCs to understand clustering in networks. We took advantage of PRC theory by considering phase-resetting of the PYR cells due to incoming inhibitory input. In this way, we were able to hypothesize an inhibition-based tuning mechanism for control of the theta rhythm frequency based on the PRC shape (amount of advance or delay) and the PYR cell's intrinsic firing frequency. Our use of PRCs relied on our observations of the effect of different PRC shapes on the resulting theta rhythm. For example, such a consideration was used by Rich et al. (2016) to explain differential synchrony patterns in inhibitory networks of Type I vs Type II neurons, and Pervouchine et al. (2006) used spike time response curve methods to analytically explore theta-generating circuits dependent on biophysical characteristics of oriens-lacunosum molecular inhibitory cell types. While theory is strongly needed and analytical considerations should be exploited whenever possible, at the same time, this should not limit the integration of experimental constraints (Skinner, 2013). We note that the work here did not start from theoretical considerations such as Type I/II cell types which can yield insights, but rather, due to the opportunity of being able to obtain simultaneous cellular and population output, multi-level experimental constraints were used in building and designing the E-I microcircuit models from the beginning.

In conclusion, we have developed a hypothesis for how theta rhythm frequencies are controlled in the CA1 hippocampus. This hypothesis is built on the theta-generating mechanism of the microcircuit model design. Even though it does not include all of the known inhibitory cell types, it perhaps captures essential elements in play in biological circuits and may apply more widely in the brain regarding the generation and control of theta rhythm frequencies.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author Contributions

FS, JL, and AC contributed to conception and supervision of the study. SR and AL performed computations and analyses. FS wrote the first draft of the manuscript. SR wrote sections of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

Funding

This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery Grant RGPIN-2016-06182 (FS).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

Parts of this work have been released as a preprint (Chatzikalymniou et al., 2020).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncir.2021.643360/full#supplementary-material

References

Achuthan, S., and Canavier, C. C. (2009). Phase-resetting curves determine synchronization, phase locking, and clustering in networks of neural oscillators. J. Neurosci. 29, 5218–5233. doi: 10.1523/JNEUROSCI.0426-09.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Ascoli, G. A., Gasparini, S., Medinilla, V., and Migliore, M. (2010). Local control of postinhibitory rebound spiking in CA1 pyramidal neuron dendrites. J. Neurosci. 30, 6434–6442. doi: 10.1523/JNEUROSCI.4066-09.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Bezaire, M. J., Raikov, I., Burk, K., Vyas, D., and Soltesz, I. (2016). Interneuronal mechanisms of hippocampal theta oscillation in a full-scale model of the rodent CA1 circuit. eLife 5:e18566. doi: 10.7554/eLife.18566

PubMed Abstract | CrossRef Full Text | Google Scholar

Bland, B. H. (1986). The physiology and pharmacology of hippocampal formation theta rhythms. Prog. Neurobiol. 26, 1–54. doi: 10.1016/0301-0082(86)90019-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsáki, G. (2002). Theta oscillations in the hippocampus. Neuron 33, 325–340. doi: 10.1016/S0896-6273(02)00586-X

CrossRef Full Text | Google Scholar

Chatzikalymniou, A. P., Gumus, M., Lunyov, A. R., Rich, S., Lefebvre, J., and Skinner, F. K. (2020). Linking minimal and detailed models of CA1 microcircuits reveals how theta rhythms emerge and how their frequencies are controlled. bioRxiv. 2020.07.28.225557. doi: 10.1101/2020.07.28.225557

CrossRef Full Text | Google Scholar

Colgin, L. L. (2013). Mechanisms and functions of theta rhythms. Annu. Rev. Neurosci. 36, 295–312. doi: 10.1146/annurev-neuro-062012-170330

CrossRef Full Text | Google Scholar

Colgin, L. L. (2016). Rhythms of the hippocampal network. Nat. Rev. Neurosci. 17, 239–249. doi: 10.1038/nrn.2016.21

CrossRef Full Text | Google Scholar

Fellous, J. M., and Sejnowski, T. J. (2000). Cholinergic induction of oscillations in the hippocampal slice in the slow (0.5–2 Hz), theta (5–12 Hz), and gamma (35–70 Hz) bands. Hippocampus 10, 187–197. doi: 10.1002/(SICI)1098-1063(2000)10:2<187::AID-HIPO8>3.0.CO;2-M

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferguson, K. A., Chatzikalymniou, A. P., and Skinner, F. K. (2017). Combining theory, model, and experiment to explain how intrinsic theta rhythms are generated in an in vitro whole hippocampus preparation without oscillatory inputs. eNeuro 4. doi: 10.1523/ENEURO.0131-17.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferguson, K. A., Huh, C. Y. L., Amilhon, B., Williams, S., and Skinner, F. K. (2013). Experimentally constrained CA1 fast-firing parvalbumin-positive interneuron network models exhibit sharp transitions into coherent high frequency rhythms. Front. Comput. Neurosci. 7:144. doi: 10.3389/fncom.2013.00144

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferguson, K. A., Huh, C. Y. L., Amilhon, B., Williams, S., and Skinner, F. K. (2015a). Simple, biologically-constrained CA1 pyramidal cell models using an intact, whole hippocampus context. F1000Research 3:104. doi: 10.12688/f1000research.3894.2

PubMed Abstract | CrossRef Full Text

Ferguson, K. A., Njap, F., Nicola, W., Skinner, F. K., and Campbell, S. A. (2015b). Examining the limits of cellular adaptation bursting mechanisms in biologically-based excitatory networks of the hippocampus. J. Comput. Neurosci. 39, 289–309. doi: 10.1007/s10827-015-0577-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferguson, K. A., and Skinner, F. K. (2018). “Hippocampal Theta, gamma, and theta/gamma network models,” in Encyclopedia of Computational Neuroscience, eds D. Jaeger and R. Jung (New York, NY: Springer New York), 1–14. doi: 10.1007/978-1-4614-7320-6_27-2

CrossRef Full Text

Gillies, M. J., Traub, R. D., LeBeau, F. E. N., Davies, C. H., Gloveli, T., Buhl, E. H., et al. (2002). A model of atropine-resistant theta oscillations in rat hippocampal area CA1. J. Physiol. 543, 779–793. doi: 10.1113/jphysiol.2002.024588

PubMed Abstract | CrossRef Full Text | Google Scholar

Gjorgjieva, J., Drion, G., and Marder, E. (2016). Computational implications of biophysical diversity and multiple timescales in neurons and synapses for circuit performance. Curr. Opin. Neurobiol. 37, 44–52. doi: 10.1016/j.conb.2015.12.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Gloveli, T., Dugladze, T., Rotstein, H. G., Traub, R. D., Monyer, H., Heinemann, U., et al. (2005). Orthogonal arrangement of rhythm-generating microcircuits in the hippocampus. Proc. Natl. Acad. Sci. U.S.A. 102, 13295–13300. doi: 10.1073/pnas.0506259102

PubMed Abstract | CrossRef Full Text | Google Scholar

Goutagny, R., Jackson, J., and Williams, S. (2009). Self-generated theta oscillations in the hippocampus. Nat. Neurosci. 12, 1491–1493. doi: 10.1038/nn.2440

PubMed Abstract | CrossRef Full Text | Google Scholar

Hansel, D., Mato, G., and Meunier, C. (1995). Synchrony in excitatory neural networks. Neural Comput. 7, 307–337. doi: 10.1162/neco.1995.7.2.307

CrossRef Full Text | Google Scholar

Hinman, J. R., Dannenberg, H., Alexander, A. S., and Hasselmo, M. E. (2018). Neural mechanisms of navigation involving interactions of cortical and subcortical structures. J. Neurophysiol. 119, 2007–2029. doi: 10.1152/jn.00498.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

Huh, C. Y. L., Amilhon, B., Ferguson, K. A., Manseau, F., Torres-Platas, S. G., Peach, J. P., et al. (2016). Excitatory inputs determine phase-locking strength and spike-timing of CA1 stratum oriens/alveus parvalbumin and somatostatin interneurons during intrinsically generated hippocampal theta rhythm. J. Neurosci. 36, 6605–6622. doi: 10.1523/JNEUROSCI.3951-13.2016

PubMed Abstract | CrossRef Full Text | Google Scholar

Izhikevich, E. M. (2006). Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting, 1st Edn. Cambridge, MA: The MIT Press. doi: 10.7551/mitpress/2526.001.0001

CrossRef Full Text | Google Scholar

Jackson, J., Amilhon, B., Goutagny, R., Bott, J.-B., Manseau, F., Kortleven, C., et al. (2014). Reversal of theta rhythm flow through intact hippocampal circuits. Nat. Neurosci. 17, 1362–1370. doi: 10.1038/nn.3803

PubMed Abstract | CrossRef Full Text | Google Scholar

Jacobs, J. (2014). Hippocampal theta oscillations are slower in humans than in rodents: implications for models of spatial navigation and memory. Philos. Trans. R. Soc. Lond. B 369:20130304. doi: 10.1098/rstb.2013.0304

PubMed Abstract | CrossRef Full Text | Google Scholar

Kopell, N., Börgers, C., Pervouchine, D., Malerba, P., and Tort, A. (2010). “Gamma and theta rhythms in biophysical models of hippocampal circuits,” in Hippocampal Microcircuits, Springer Series in Computational Neuroscience, eds V. Cutsuridis, B. Graham, S. Cobb, and I. Vida (New York, NY: Springer), 423–457. doi: 10.1007/978-1-4419-0996-1_15

CrossRef Full Text | Google Scholar

Kota, S., Rugg, M. D., and Lega, B. C. (2020). Hippocampal theta oscillations support successful associative memory formation. J. Neurosci. 40, 9507–9518. doi: 10.1523/JNEUROSCI.0767-20.2020

PubMed Abstract | CrossRef Full Text | Google Scholar

Kramis, R., Vanderwolf, C. H., and Bland, B. H. (1975). Two types of hippocampal rhythmical slow activity in both the rabbit and the rat: relations to behavior and effects of atropine, diethyl ether, urethane, and pentobarbital. Exp. Neurol. 49, 58–85. doi: 10.1016/0014-4886(75)90195-8

PubMed Abstract | CrossRef Full Text | Google Scholar

McCormick, D. A., McGinley, M. J., and Salkoff, D. B. (2015). Brain state dependent activity in the cortex and thalamus. Curr. Opin. Neurobiol. 31, 133–140. doi: 10.1016/j.conb.2014.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Middleton, J. W., Omar, C., Doiron, B., and Simons, D. J. (2012). Neural correlation is stimulus modulated by feedforward inhibitory circuitry. J. Neurosci. 32, 506–518. doi: 10.1523/JNEUROSCI.3474-11.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Mikulovic, S., Restrepo, C. E., Siwani, S., Bauer, P., Pupe, S., Tort, A. B. L., et al. (2018). Ventral hippocampal OLM cells control type 2 theta oscillations and response to predator odor. Nat. Commun. 9:3638. doi: 10.1038/s41467-018-05907-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Orbán, G., Kiss, T., and Érdi, P. (2006). Intrinsic and synaptic mechanisms determining the timing of neuron population activity during hippocampal theta oscillation. J. Neurophysiol. 96, 2889–2904. doi: 10.1152/jn.01233.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Pervouchine, D. D., Netoff, T. I., Rotstein, H. G., White, J. A., Cunningham, M. O., Whittington, M. A., et al. (2006). Low-dimensional maps encoding dynamics in entorhinal cortex and hippocampus. Neural Comput. 18, 2617–2650. doi: 10.1162/neco.2006.18.11.2617

PubMed Abstract | CrossRef Full Text | Google Scholar

Rich, S., Booth, V., and Zochowski, M. (2016). Intrinsic cellular properties and connectivity density determine variable clustering patterns in randomly connected inhibitory neural networks. Front. Neural Circ. 10:82. doi: 10.3389/fncir.2016.00082

PubMed Abstract | CrossRef Full Text | Google Scholar

Rich, S., Zochowski, M., and Booth, V. (2017). Dichotomous dynamics in EI networks with strongly and weakly intra-connected inhibitory neurons. Front. Neural Circ. 11:104. doi: 10.3389/fncir.2017.00104

CrossRef Full Text | Google Scholar

Schultheiss, N., Butera, R., and Prinz, A. (2011). Phase Response Curves in Neuroscience: Theory, Experiment, and Analysis. Springer. doi: 10.1007/978-1-4614-0739-3

CrossRef Full Text | Google Scholar

Skinner, F. K. (2013). Moving beyond type I and type II neuron types. F1000Research 2:19. doi: 10.12688/f1000research.2-19.v1

PubMed Abstract | CrossRef Full Text | Google Scholar

Tendler, A., and Wagner, S. (2015). Different types of theta rhythmicity are induced by social and fearful stimuli in a network associated with social memory. eLife 4:e03614. doi: 10.7554/eLife.03614

PubMed Abstract | CrossRef Full Text | Google Scholar

ter Wal, M., and Tiesinga, P. (2013). “Hippocampal oscillations, mechanisms (PING, ING, sparse),” in Encyclopedia of Computational Neuroscience, eds D. Jaeger and R. Jung (New York, NY: Springer New York), 1–14. doi: 10.1007/978-1-4614-7320-6_475-3

CrossRef Full Text