Theta-Nested Gamma Oscillations in Next Generation Neural Mass Models

Theta-nested gamma oscillations have been reported in many areas of the brain and are believed to represent a fundamental mechanism to transfer information across spatial and temporal scales. In a series of recent experiments in vitro it has been possible to replicate with an optogenetic theta frequency stimulation several features of cross-frequency coupling (CFC) among theta and gamma rhythms observed in behaving animals. In order to reproduce the main findings of these experiments we have considered a new class of neural mass models able to reproduce exactly the macroscopic dynamics of spiking neural networks. In this framework, we have examined two set-ups able to support collective gamma oscillations: namely, the pyramidal interneuronal network gamma (PING) and the interneuronal network gamma (ING). In both set-ups we observe the emergence of theta-nested gamma oscillations by driving the system with a sinusoidal theta-forcing in proximity of a Hopf bifurcation. These mixed rhythms always display phase amplitude coupling. However, two different types of nested oscillations can be identified: one characterized by a perfect phase locking between theta and gamma rhythms, corresponding to an overall periodic behavior; another one where the locking is imperfect and the dynamics is quasi-periodic or even chaotic. From our analysis it emerges that the locked states are more frequent in the ING set-up. In agreement with the experiments, we find theta-nested gamma oscillations for forcing frequencies in the range [1:10] Hz, whose amplitudes grow proportionally to the forcing intensity and which are clearly modulated by the theta phase. Furthermore, analogously to the experiments, the gamma power and the frequency of the gamma-power peak increase with the forcing amplitude. At variance with experimental findings, the gamma-power peak does not shift to higher frequencies by increasing the theta frequency. This effect can be obtained, in our model, only by incrementing, at the same time, also the stimulation power. An effect achieved by increasing the amplitude either of the noise or of the forcing term proportionally to the theta frequency. On the basis of our analysis both the PING and the ING mechanism give rise to theta-nested gamma oscillations with almost identical features.

1. Introduction

Oscillations in the brain, reflecting the underlying dynamics of neural populations, have been measured over a broad frequency range (Buzsaki, 2006). Particularly studied are γ-rhythms (30–120 Hz), due to their ubiquitous presence in many regions of the brain, irrespectively of the species (Buzsáki and Wang, 2012), and for their relevance for cognitive tasks (Fries et al., 2007) and neuronal diseases (Uhlhaas and Singer, 2006; Williams and Boksa, 2010).

Inhibitory networks have been shown to represent a fundamental ingredient for the emergence of γ oscillations (Bartos et al., 2007; Buzsáki and Wang, 2012). Indeed, inhibition is at the basis of the two most known mechanisms: pyramidal interneuronal network gamma (PING) and interneuronal network gamma (ING) (Tiesinga and Sejnowski, 2009). The ING mechanism is observable in purely inhibitory networks in the presence of few ingredients: recurrent connections, a time scale associated with the synaptic GABA_A receptors and an excitatory drive sufficiently strong to lead the neurons supra-threshold (Buzsáki and Wang, 2012). The collective oscillations (COs) emerge when a sufficient number of neurons begins to fire within a short time window and generate almost synchronous inhibitory post-synaptic potentials (IPSPs) in the post-synaptic interneurons. The inhibited neurons fire again when the IPSPs have sufficiently decayed and the cycle will repeat. Thus, the main ingredients dictating the frequency of the COs in the ING set-up are: the kinetics of the IPSPs and the excitatory drive (Whittington et al., 1995). On the other hand the PING mechanism is related to the presence of an excitatory and an inhibitory population, in this case COs emerge whenever the drive on the excitatory neurons is sufficiently strong to induce an almost synchronous excitatory volley that in turn elicits an inhibitory one. The period of the COs is thus determined by the recovery time of the pyramidal neurons from the stimulus received from the inhibitory population (Wilson and Cowan, 1972). A peculiarity of this mechanism, observed both in vivo and in vitro experiments, is that there is a delay between the firing of the pyramidal cells and the interneuronal burst (Buzsáki and Wang, 2012).

In several parts of the brain, one can observe that γ oscillations are modulated by θ oscillations, with θ frequencies corresponding to 4–12 Hz in rodents and to 1–4 Hz in humans. Specific examples have been reported for the hippocampus of rodents in behaving animals and during rapid eye movement (REM) sleep (Lisman, 2005; Colgin et al., 2009; Belluscio et al., 2012; Perńıa-Andrade and Jonas, 2014; Colgin, 2015), for the visual cortex in alert monkeys (Whittingstall and Logothetis, 2009), for the neocortex in humans (Canolty et al., 2006) etc. This is an example of a more general mechanism of cross-frequency coupling (CFC) between a low and a high frequency rhythm, which is believed to have a functional role in the brain (Canolty and Knight, 2010). In particular, low frequency rhythms (such as θ) are usually involving broad brain regions and are entrained to external inputs and/or to cognitive events; on the other hand the high frequency activity (e.g., the γ-rhythm) reflects local computation activity. Thus CFC can represent an effective mechanism to transfer information across spatial and temporal scales (Canolty and Knight, 2010; Lisman and Jensen, 2013). Four different types of CFC of interest for electrophysiology, have been listed in Jensen and Colgin (2007): phase-phase, phase-frequency, phase-amplitude and amplitude-amplitude couplings (PPC, PFC, PAC, and AAC). Two more types of CFCs have later been added as emerging from the analysis of coupled non-linear oscillators (Witte et al., 2008) and coupled neural mass models (Chehelcheraghi et al., 2017): frequency-frequency and amplitude-frequency coupling (FFC and AFC).

In this paper, we will consider θ-nested γ oscillations, where specific features of the γ oscillations are correlated to the θ phase. In particular, we will analyze PPC, PFC, and PAC between θ and γ rhythms. The most studied CFC mechanism is the PAC, which corresponds to the modification of the amplitude (or power) of γ-waves induced by the phase of the θ-oscillations. This mechanism has been reported in the primary visual cortex of anaesthetized macaques subject to naturalistic visual stimulation (Mazzoni et al., 2011), as well as during the formation of new episodic memories in the human hippocampus (Lega et al., 2016). As discussed in Jensen and Colgin (2007), the θ phase can often modulate both amplitude (PAC) and frequency (PFC) of the γ oscillations, therefore these two mechanisms can occur at the same time. PPC, which refers to n:m phase locking between γ and θ phase oscillations (Tass et al., 1998), has been identified in the rodent hippocampus during maze exploration (Belluscio et al., 2012).

Our study is mostly motivated by recent optogenetic experiments revealing PAC in areas CA1 and CA3 of the hippocampus and in the medial enthorinal cortex (MEC) (Akam et al., 2012; Pastoll et al., 2013; Butler et al., 2016, 2018). These experiments have shown that a sinusoidal optogenetic stimulation at θ-frequency of the circuits in vitro is able to reproduce several features of θ-nested γ oscillations, usually observed in behaving rats (Bragin et al., 1995). All these experiments suggest that inhibition has a key role in generating this cross-frequency rhythm; however both ING (Pastoll et al., 2013) and PING (Butler et al., 2016, 2018) mechanisms have been invoked to explain locally generated γ oscillations.

PING and ING oscillation mechanisms have been qualitatively reproduced by employing heuristic neural mass models (Wilson and Cowan, 1972; Gerstner et al., 2014). However, these standard firing rate models do not properly describe the synchronization and desynchronizaton phenomena occurring in neural populations (Devalle et al., 2017; Laing, 2017; Coombes and Byrne, 2019). Recently a new generation of neural mass models has been designed, which are able to exactly reproduce the network dynamics of spiking neurons of class I, for any degree of synchronization among the neurons (Luke et al., 2013; Laing, 2014; So et al., 2014; Montbrió et al., 2015). In particular, for purely inhibitory networks, these mean-field models have been able to reproduce the emergence of COs, observed in the corresponding networks, without the inclusion of an extra time delay (Devalle et al., 2017), as well as the phenomenon of event related synchronization and desynchronization (Coombes and Byrne, 2019).

Our main aim is to understand how θ-nested γ oscillations can emerge when a PING or ING mechanism is responsible for the fast oscillations and which differences can be expected in the population dynamics in the two cases. Therefore we will consider the new class of neural mass models introduced in Montbrió et al. (2015) in two configurations: namely, a purely inhibitory population (ING set-up) and two coupled excitatory-inhibitory populations (PING set-up). In both configurations we will examine the system response to an external sinusoidal θ-drive.

Section 2 is devoted to the introduction of different spiking network configurations of Quadratic Integrate-and-Fire (QIF) neurons able to generate γ COs via PING and ING mechanisms and to the introduction of their corresponding exact neural mass formulations. A detailed bifurcation analysis of the neural mass models for the PING and ING set-ups, in the absence of any external forcing, is reported in section 3. The PAC mechanism is analyzed and discussed in section 4. First, by considering different types of PAC states (namely, phase locked or unlocked) and second, by comparing our numerical results for PAC dynamics with experimental findings reported in Butler et al. (2016, 2018), for the CA1 region of the hippocampus under sinusoidal optogenetic stimulations. Finally, a discussion of our results and of their implications, as well as of possible future developments, will be presented in section 5. The results reported in the paper are mostly devoted to super-critical Hopf bifurcations, however a specific example of a sub-critical Hopf bifurcation leading to COs is discussed for the PING set-up in Appendix A. Further network configurations ensuring the emergence of COs via PING mechanism are presented in Appendix B.

2. Models and Bifurcation Analysis

2.1. Network Models

In this paper we want to compare the two principal mechanisms at the basis of the emergence of collective oscillatory dynamics in neural networks: namely, the PING and ING mechanisms. Therefore we will consider QIF neurons in the two following set-ups: an excitatory and an inhibitory population coupled via instantaneous synapses (PING configuration) and a single inhibitory population interacting via post-synaptic potentials (PSPs) with exponential profile (ING configuration). The corresponding network configurations are shown in Figure 1. Moreover, the neurons are assumed to be fully coupled. As we will show in the following, both these two configurations support the emergence of COs.

FIGURE 1

Figure 1. Network topologies. Two different network configurations have been investigated: on the left side, an excitatory population (E) and an inhibitory population (I) form a circuit that can generate oscillatory output (PING set-up); on the right side one inhibitory population (I) is coupled to itself with an inhibitory coupling (ING set-up). In both cases an external current I^(l) impinging on one single population has been considered.

The dynamics of the membrane potentials of the QIF neurons in the PING configuration is given by

$\begin{array}{l}{\tau }_m^{(e)}{\stackrel{\cdot }{V}}_k^{(e)}={\left(V_k^{(e)}\right)}^2+{\eta }_k^{(e)}+{\tau }_m^{(e)}\left[J^{(ee)}{s}^{(e)}-J^{(ie)}{s}^{(i)}\right]\\ +I^{(e)}(t)k=1,\dots ,N^{(e)}\\ {\tau }_m^{(i)}{\stackrel{\cdot }{V}}_j^{(i)}={\left(V_j^{(i)}\right)}^2+{\eta }_j^{(i)}+{\tau }_m^{(i)}\left[J^{(ei)}{s}^{(e)}-J^{(ii)}{s}^{(i)}\right]\\ +I^{(i)}(t)j=1,\dots ,N^{(i)}\\ s^{(l)}=\frac{1}{N^{(l)}}{\displaystyle \sum _{t_m^{(l)}}}\delta (t-t_m^{(l)})l\in \left\{e,i\right\};& (1)\end{array}$

where the super-scripts e (i) denote the excitatory (inhibitory) population, ${\tau }_m^{(e)}=20$ ms (${\tau }_m^{(i)}=10$ ms) is the excitatory (inhibitory) membrane time constant, ${\eta }_k^{(l)}$ is the excitability of the k-th neuron of population l, J^(ln) is the strength of the synaptic coupling of population l acting on population n. The term I^(l)(t) represents a time-dependent external current applied to the population l; usually we have considered the external drive to be applied to the excitatory population only, i.e., I^(e)(t) ≠ 0 and I⁽ⁱ⁾(t) = 0. The synaptic field s^(l)(t) is the linear super-position of all the pulses p(t) emitted in the past within the l population, p(t) being δ-functions in the present case. Furthermore, since the neurons are fully coupled, each neuron will be subject to the same synaptic field (Olmi et al., 2010). The emission of the m-th spike in the network occurs at time $t_m^{(l)}$ whenever the membrane potential of a generic neuron j reaches infinity, i.e., $V_j^{(l)}\left(t_m^{{(l)}_{-}}\right)\to +\infty$, while the reset mechanism is modeled by setting $V_j^{(l)}\left(t_j^{{(l)}_{+}}\right)\to -\infty$, immediately after the spike emission.

The main part of our analysis of the PING set-up will be devoted to networks with self-activation only (i.e., where J⁽ⁱⁱ⁾ = 0), a configuration which is known to favor the emergence of collective oscillations (Wilson and Cowan, 1972; Kilpatrick, 2015; Onslow et al., 2014). However, as discussed in Appendix B, we have found that COs can arise in different PING set-ups: in the presence of self-inhibition only (i.e., with J⁽ⁱⁱ⁾ ≠ 0 and J^(ee) = 0) and in the absence of both self-activation and inhibition (i.e., with J^(ee) = J⁽ⁱⁱ⁾ = 0).

For what concerns the purely inhibitory network, the membrane potential dynamics of the j-th neuron is ruled by the following equations:

$\begin{array}{l}{\tau }_m^{(i)}{\stackrel{\cdot }{V}}_j^{(i)}={\left(V_j^{(i)}\right)}^2+{\eta }_j^{(i)}-{\tau }_m^{(i)}{J}^{(ii)}{s}^{(i)}+I^{(i)}(t)\\ {\tau }_d{\dot{s}}^{(i)}=-s^{(i)}+\frac{1}{N^{(i)}}{\displaystyle \sum _{t_m^{(l)}}}\delta (t-t_m^{(i)}),& (2)\end{array}$

where ${\tau }_m^{(i)}=10$ ms. In this case the synaptic field s⁽ⁱ⁾(t) is the super-position of the exponential IPSPs $p(t)={\text{e}}^{-t/{\tau }_d}/{\tau }_d$ emitted in the past, where we set τ_d = 10 ms.

For reasons that will become clear in the next paragraph, we assume that the neuron excitabilities ${\eta }_i^{(l)}$ are randomly distributed according to a Lorentzian probability density function (PDF)

$\begin{array}{ll}{g}^{(l)}(\eta )=\frac{1}{\pi }\frac{{\Delta }^{(l)}}{{\left(\eta -H^{(l)}\right)}^2+{\left({\Delta }^{(l)}\right)}^2},& (3)\end{array}$

where H^(l) is the median and Δ^(l) is the half-width half-maximum (HWHM) of the PDF. Therefore each population will be composed of neurons supra- (with ${\eta }_j^{(l)}>0$) and sub-threshold (with ${\eta }_j^{(l)}<0$), the percentage of one group with respect to the other being determined by the Lorentzian parameters. For the PING set-up we fix Δ^(e) = Δ⁽ⁱ⁾ = 1, whereas varying H^(e) and H⁽ⁱ⁾. For the ING set-up we fix Δ⁽ⁱ⁾ = 0.3 and analyze the dynamics by varying H⁽ⁱ⁾.

The dynamical equations are integrated by employing a 4th order Runge-Kutta method in the absence of noise with a time step dt = 0.002 ms (dt = 0.001 ms) for the PING (ING) set-up. Moreover, we define a threshold V_p = 100 and a reset value V_r = −100. Whenever the membrane potential V_j of the j-th neuron overcomes V_p at a time t_p, it is reset to V_r for a refractory period equal to 2/V_j. At the same time the firing time is estimated as t_p + 1/V_j; for more details see Montbrió et al. (2015). The membrane potentials are initialized from a random flat distribution defined over the range [−100:100], while the excitabilities are randomly chosen from the Lorentzian distribution (3).

For instantaneous synapses, we will only employ the following two indicators to characterize the macroscopic dynamics:

$\begin{array}{ll}{r}^{(l)}(t)=\frac{M^{(l)}(\Delta t)}{N^{(l)}\Delta t},v^{(l)}(t)=\frac{1}{N^{(l)}}{\displaystyle \sum _j^{N^{(l)}}}{V}_j^{(l)}(t),& (4)\end{array}$

which represent the average population activity and the average membrane potential of the l-th population, respectively. In particular, the average population activity of the l−network r^(l)(t) is given by the number of spikes M^(l)(Δt) emitted in a time window Δt, divided by the total number of neurons in such population. For finite IPSPs we also consider the synaptic field s^(l)(t). Furthermore, the emergence of COs, corresponding to periodic motions of r^(l)(t) and v^(l)(t), are characterized in terms of their frequencies ν^(l).

We assume that the driving current, mimicking the θ-stimulation in the optogenetic experiments, is a purely sinusoidal excitatory current of the following form

$\begin{array}{ll}{I}_{\theta }(t)=\frac{I_0}{2}\left[1-cos(2\pi {\nu }_{\theta }t)\right]& (5)\end{array}$

where ν_θ is the forcing frequency, usually considered within the θ-range, i.e., ν_θ ∈ [1:10] Hz. In this context a theta phase associated with the forcing field can be defined as θ(t) = mod(2πν_θt, 2π). For the PING configuration we set $I^{(e)}(t)=I_{\theta }(t)$ and I⁽ⁱ⁾(t) ≡ 0 and for the ING set-up $I^{(i)}(t)=I_{\theta }(t)$.

2.2. Neural mass models

As already mentioned, an exact neural mass model has been derived in Montbrió et al. (2015) for a fully coupled network of QIF neurons with instantaneous synapses and with Lorentzian distributed neuronal excitabilities. In this case the macroscopic neural dynamics of a population l is described by two collective variables: the mean field potential v^(l)(t) and the instantaneous firing rate r^(l)(t). In this context, the neural mass model for two coupled E − I populations with instantaneous synapses, corresponding to the microscopic model reported in Equation (1), can be written as

$\begin{array}{l}{\dot{r}}^{(e)}=\frac{{\Delta }^{(e)}}{{\left({\tau }_m^{(e)}\right)}^2\pi }+\frac{2r^{(e)}{v}^{(e)}}{{\tau }_m^{(e)}}\\ {\stackrel{\cdot }{v}}^{(e)}=\frac{{\left(v^{(e)}\right)}^2+H^{(e)}+I^{(e)}(t)}{{\tau }_m^{(e)}}-{\tau }_m^{(e)}{\left(\pi r^{(e)}\right)}^2\\ +J^{(ee)}{r}^{(e)}-J^{(ie)}{r}^{(i)}+A{\xi }^{(e)}\\ {\dot{r}}^{(i)}=\frac{{\Delta }^{(i)}}{{\left({\tau }_m^{(i)}\right)}^2\pi }+\frac{2r^{(i)}{v}^{(i)}}{{\tau }_m^{(i)}}\\ {\stackrel{\cdot }{v}}^{(i)}=\frac{{\left(v^{(i)}\right)}^2+H^{(i)}+I^{(i)}(t)}{{\tau }_m^{(i)}}-{\tau }_m^{(i)}{\left(\pi r^{(i)}\right)}^2\\ +J^{(ei)}{r}^{(e)}-J^{(ii)}{r}^{(i)}+A{\xi }^{(i)}.& (6)\end{array}$

In the equations for the evolution of the average membrane potentials we have also inserted an additive noise term of amplitude A, employed in some of the analysis to mimic the many noise sources present in the brain dynamics. In particular, the noise terms ξ^(e) and ξ⁽ⁱ⁾ are both δ-correlated and uniformly distributed in the interval [−1:1].

In case of finite synapses, the exact derivation of the corresponding neural mass model is still feasible for QIF neurons, but the macroscopic evolution now contains further equations describing the dynamics of the synaptic field characterizing the considered synapses (Devalle et al., 2017; Coombes and Byrne, 2019). In particular, for a single inhibitory population with exponential synapses, the corresponding neural mass model reads as:

$\begin{array}{l}{\dot{r}}^{(i)}=\frac{{\Delta }^{(i)}}{{({\tau }_m^{(i)})}^2\pi }+\frac{2r^{(i)}{v}^{(i)}}{{\tau }_m^{(i)}}\\ {\stackrel{\cdot }{v}}^{(i)}=\frac{{(v^{(i)})}^2+H^{(i)}+I^{(i)}(t)}{{\tau }_m^{(i)}}-{\tau }_m^{(i)}{(\pi r^{(i)})}^2\\ -J^{(ii)}{s}^{(i)}+A{\xi }^{(i)}\\ {\dot{s}}^{(i)}=\frac{1}{{\tau }_d}\left[-s^{(i)}+r^{(i)}\right].& (7)\end{array}$

In the present case the equation for the average membrane potential contains, as already shown before in Equation (6), an additive noise term of amplitude A.

It should be noticed that in Equations (6) and (7) the noise has been added in an effective manner and not with a consistent procedure, that would amount to take into account the effect of microscopic noise on the mean-field formulation. This can be achieved by considering a Fokker-Planck description for the distribution of the membrane potentials, e.g., as done in Brunel and Hakim (1999), or by considering a reduced approach in terms of circular cumulants (Tyulkina et al., 2018; Goldobin et al., 2018). However, all these formulations will lead to much more complicated evolution equations for the macroscopic quantities.

To analyse the stability of the macroscopic solutions of Equations (6) and (7), one should estimate the corresponding Lyapunov spectrum (LS) (Pikovsky and Politi, 2016). This can be done by considering the time evolution of the tangent vector, which for the PING set-up turns out to be four dimensional, i.e., δ = {δr^(e), δv^(e), δr⁽ⁱ⁾δv⁽ⁱ⁾}. The dynamics of the tangent vector is ruled by the linearization of the Equation (6), namely

$\begin{array}{l}\delta {\dot{r}}^{(e)}=\frac{2\left(r^{(e)}\delta v^{(e)}+v^{(e)}\delta r^{(e)}\right)}{{\tau }_m^{(e)}}\\ \delta {\stackrel{\cdot }{v}}^{(e)}=\frac{2v^{(e)}\delta v^{(e)}}{{\tau }_m^{(e)}}-2{\tau }_m^{(e)}{\pi }^2{r}^{(e)}\delta r^{(e)}+J^{(ee)}\delta r^{(e)}-J^{(ie)}\delta r^{(i)}\\ \delta {\dot{r}}^{(i)}=\frac{2\left(r^{(i)}\delta v^{(i)}+v^{(i)}\delta r^{(i)}\right)}{{\tau }_m^{(i)}}\\ \delta {\stackrel{\cdot }{v}}^{(i)}=\frac{2v^{(i)}\delta v^{(i)}}{{\tau }_m^{(i)}}-2{\tau }_m^{(i)}{\pi }^2{r}^{(i)}\delta r^{(i)}+J^{(ei)}\delta r^{(e)}-J^{(ii)}\delta r^{(i)}.& (8)\end{array}$

For the ING set-up the tangent vector is three dimensional, δ = {δr⁽ⁱ⁾, δv⁽ⁱ⁾, δs⁽ⁱ⁾}, and its time evolution can be obtained by the linearization of Equation (7), which reads as

$\begin{array}{l}\delta {\dot{r}}^{(i)}=\frac{2\left(r^{(i)}\delta v^{(i)}+v^{(i)}\delta r^{(i)}\right)}{{\tau }_m^{(i)}}\\ \delta {\stackrel{\cdot }{v}}^{(i)}=\frac{2v^{(i)}\delta v^{(i)}}{{\tau }_m^{(i)}}-2{\tau }_m^{(i)}{\pi }^2{r}^{(i)}\delta r^{(i)}-J^{(ii)}\delta s^{(i)}\\ \delta {ṡ}^{(i)}=\frac{1}{{\tau }_d}\left[-\delta s^{(i)}+\delta r^{(i)}\right].& (9)\end{array}$

Please notice that the presence of additive external noise or of forcing terms in Equations (6) and (7) does not modify the evolution equations in the tangent space Equations (8) and (9).

The LS is composed by 4 (3) Lyapunov exponents (LEs) {λ_i} for the PING (ING) set-ups, which quantify the average growth rates of infinitesimal perturbations along the orthogonal manifolds. In details, LEs are estimated as follows

$\begin{array}{ll}{\lambda }_i={\displaystyle \underset{t\to \infty }{lim}}\frac{1}{t}log\frac{\left|{\delta }_i(t)\right|}{\left|{\delta }_i(0)\right|},& (10)\end{array}$

where the technique described in Benettin et al. (1980) to maintain the tangent vectors δ_i orthonormal during the evolution is employed. The autonomous system will be chaotic for λ₁ > 0, while a periodic (quasi-periodic) dynamics will be characterized by λ₁ = 0 (λ₁ = λ₂ = 0) and a fixed point by λ₁ < 0. In a non-autonomous system in the presence of an external forcing, one Lyapunov exponent will be necessarily zero, therefore a periodic behavior corresponds to λ₁ < 0 and a quasi-periodic dynamics to λ₁ = 0 (Pikovsky and Politi, 2016).

In the absence of noise, neural mass models have been directly integrated by employing a Runge-Kutta 4th order integration scheme, while in the presence of additive noise with a Heun scheme. In both cases the time step has been set to dt = 0.01 ms. In order to estimate the Lyapunov spectra we have integrated the direct and tangent space evolution with a Runge-Kutta 4th order integration scheme with dt = 0.001 ms, for a duration of 200 s, after discarding a transient of 10 s.

Besides LEs, in order to characterize the macroscopic dynamics of the model, we have estimated the frequency power spectra $P_S^{(e)}(F)$ ($P_S^{(i)}(F)$) of the mean excitatory (inhibitory) membrane potential v^(e)(t) (v⁽ⁱ⁾(t)) for the PING (ING) set-up. The power spectra have been obtained by calculating the temporal Fourier transform of the mean membrane potentials sampled at time intervals of 2 ms. In the deterministic (noisy) case, time traces composed of 2048 (1024) consecutive intervals have been considered to estimate the spectra, which are obtained at a frequency resolution of ΔF = 0.244 Hz (ΔF = 0.488 Hz). Finally, the power spectra have been averaged over 12 (488) independent realizations for the deterministic (noisy) dynamics. To compare our numerical findings with the experimental results reported in Butler et al. (2016), as a measure of the power of the γ oscillations, we have estimated the area of the power spectrum P_γ in an interval ±15 Hz around the main peak position F_r of the corresponding power spectrum.

3. Dynamics in the Absence of Forcing

Due to the low dimensionality of the neural mass models we have been able to obtain the corresponding bifurcation diagrams by employing the software MATCONT developed for orbit continuation (Govaerts et al., 2006).

In particular, we have derived the bifurcation diagrams in the absence of forcing [I^(e) = I⁽ⁱ⁾ ≡ 0] as a function of the medians H^(e) and H⁽ⁱ⁾ of the excitability distributions for the PING and ING configuration. In general, we observe either asynchronous dynamics, corresponding to a stable fixed point (a focus) of the neural mass equations, or COs, corresponding to stable limit cycles for the same set of equations.

3.1. PING set-up

For the excitatory-inhibitory set-up, as already mentioned, we usually fix H⁽ⁱ⁾ = −5 and we vary H^(e). In this case the inhibitory neurons are mostly below threshold (apart from 6 to 7% of them) and they can be driven supra-threshold from the activity of the excitatory population for sufficiently large values of H^(e). COs emerge when a sufficient number of neurons is supra-threshold, i.e., when H^(e) becomes positive enough, Indeed, as shown in Figure 2A, at negative or low values of H^(e), one can observe asynchronous dynamics, where the neurons fire independently and without any collective behavior (as an example see Figure 2C). By increasing H^(e), a supercritical Hopf bifurcation occurs at $H_c^{(e)}\simeq 1.5$ leading to the emergence of COs. The COs regime is characterized in the network by almost periodic population bursts, where the neurons in one population partially synchronize over a short time window in the order of a few milliseconds. An example for H^(e) = 5 is shown in Figure 2D, where one can observe two salient characteristics of the oscillatory dynamics. Firstly, the excitatory anticipates always the inhibitory burst by a certain time interval T_a (in this case T_a ≃ 5 ms), as usually observed for the PING mechanism (Tiesinga and Sejnowski, 2009). Secondly, the bursts of the excitatory population have a temporal width (≃8 ms) which is two or three times larger than those of the inhibitory ones (≃2−3 ms). This is also due to the fact that a large part of the inhibitory neurons is sub-threshold, therefore most of them fire within a short time window, irrespective of their excitabilities, due to the arrival of the synaptic stimulation from the excitatory population. Instead, the excitatory neurons, which are mostly supra-threshold, recover from silence, due to the inhibitory stimulation received during the inhibitory burst, over a wider time interval, driven by their own excitabilities. It is evident that the CO frequency of the excitatory and inhibitory population coincide in this set-up.

FIGURE 2

Figure 2. (PING set-up) (A) Bifurcation diagram of the average membrane potential v^(e) as a function of H^(e), for H⁽ⁱ⁾ = −5.0. The black continuous (dashed) line identifies the stable (unstable) fixed point. The red lines denote the maxima and minima of the limit cycles. The supercritical Hopf bifurcation occurs at $H_c^{(e)}=1.5$. The inset shows the frequency ν^(e) of the COs vs. H^(e). (B) Bifurcation diagram of the average membrane potential v^(e) vs. H⁽ⁱ⁾ for H^(e) = 10. The Hopf bifurcations are located at $H_{c1}^{(i)}=-8.4$ and $H_{c2}^{(i)}=0.20$, while the saddle-node bifurcation of limit cycles occurs at $H_{SN}^{(i)}=-10.0$. The inset show the frequency ν⁽ⁱ⁾ ≡ ν^(e) of the COs vs. H⁽ⁱ⁾. (C,D) Raster plots of the excitatory (green dots) and inhibitory (blue dots) networks are calculated in correspondence with the stable fixed point for H^(e) = −5.0 (C) and with the limit cycle for H^(e) = +5.0 (d) for the case analyzed in (A). For a better visualization, the activity of only 500 neurons of each population is shown. (E) Delay T_a as a function of H^(e). The red dashed line denotes $H_c^{(e)}$. Here we have used the same parameters as in (A). In the inset is reported the dependence of T_a vs. H⁽ⁱ⁾ for the parameters in (B). The other parameters of the system are J^(ee) = 8, J^(ie) = J^(ei) = 10, J⁽ⁱⁱ⁾ = 0 and the sizes of the networks are N^(e) = 5,000, N⁽ⁱ⁾ = 5,000.

Moreover, it is important to investigate the bifurcation diagram of the system at fixed median excitatory drive by varying H⁽ⁱ⁾. The corresponding bifurcation diagram is displayed in Figure 2B for H^(e) = 10. By increasing H⁽ⁱ⁾, COs emerge from the asynchronous state via a sub-critical Hopf bifurcation at $H_{c1}^{(i)}\simeq -8.4$ and they disappear via a super-critical Hopf bifurcation at $H_{c2}^{(i)}\simeq 0.20$. Since the first transition is hysteretical, COs disappear via a saddle-node of the limit cycles at a value $H_{SN}^{(i)}\simeq -10.00$ lower than $H_{c1}^{(i)}$. Indeed, in the interval $\left[H_{SN}^{(i)};H_{c1}^{(i)}\right]$ we have the coexistence of a stable focus with a stable limit cycle. In summary, COs are clearly observable as long as H⁽ⁱ⁾ is negative or sufficiently small. If the inhibitory neurons become mostly supra-threshold, this destroys the collective behavior associated with the PING mechanism.

It is worth noticing that the frequencies of the COs are in the γ-range, namely ν^(e) ∈ [22:71] Hz (as shown in the inset of Figure 2A): in this set-up the maximum achievable frequency ≃100 Hz, since the decay time of inhibition is dictated by ${\tau }_m^{(i)}=10$ ms (Tiesinga and Sejnowski, 2009). On the other hand, the influence of H⁽ⁱ⁾ on the frequency of the COs is quite limited. As shown in the inset of Figure 2B for a specific case corresponding to H^(e) = 10.0, ν⁽ⁱ⁾ ≡ ν^(e) varies by few Hz (namely, from 42.8 to 46.9 Hz), when H⁽ⁱ⁾ is varied by an order of magnitude.

For what concerns the delay T_a between the excitatory and inhibitory bursts, we observe a decrease of T_a with the increase of the excitatory drive H^(e), from T_a ≃ 10 ms at the Hopf bifurcation, toward 2 ms for large H^(e) value, see Figure 2E. The largest value of T_a is of the order of ${\tau }_m^{(i)}$. This can be explained by the fact that the excitatory stimulations should reach the inhibitory population within a time interval of (at most) $\simeq {\tau }_m^{(i)}$ to be able to sum up in an effective manner and to ignite the inhibitory burst. As shown in the inset of Figure 2E, the increase of H⁽ⁱ⁾ has in general the effect to reduce T_a; this should be expected since for larger excitabilities [larger H⁽ⁱ⁾], the inhibitory neurons are faster in responding to the excitatory stimulations. However, this is not the case in proximity of the saddle-node bifurcation at $H_{SN}^{(i)}$ and for positive H⁽ⁱ⁾, where the effect is reversed and T_a increases with H⁽ⁱ⁾. It is worth noticing that the same parameters as in Figure 2A are used for the main panel Figure 2E, while in the inset of Figure 2E, the data shown are calculated for the same parameters as in Figure 2B.

For the PING set-up we can also observe sub-critical Hopf bifurcations. A specific example is discussed in some detail in Appendix A.

3.2. ING set-up

As shown in Devalle et al. (2017), in order to observe COs in globally coupled inhibitory QIF networks and in the corresponding neural mass models, it is sufficient to include a finite synaptic time scale τ_d. On the other hand, in sparse balanced QIF networks, COs are observable even for instantaneous synapses (di Volo and Torcini, 2018). Indeed, for the set of parameters here employed, by varying the median of the inhibitory excitabilities H⁽ⁱ⁾, we observe a super-critical Hopf bifurcation at $H_c^{(i)}\simeq 2.4$, from an asynchronous state to COs (see Figure 3A). Analogously to the PING set-up, the frequencies of the COs observable in the ING set-up are within the γ-range, namely ν⁽ⁱ⁾ ∈ [26:83] Hz. In particular, we observe an almost linear increase of ν⁽ⁱ⁾ with H⁽ⁱ⁾. An example of the observed dynamics is shown in Figure 3B (Figure 3C) where the raster plot of the inhibitory network is calculated for H⁽ⁱ⁾ = 0(H⁽ⁱ⁾ = 10).

FIGURE 3

Figure 3. (ING set-up) (A) Bifurcation diagram of the average membrane potential v⁽ⁱ⁾ as a function of H⁽ⁱ⁾. The black continuous (dashed) line identifies the stable (unstable) fixed point. The red lines denote the maxima and minima of the limit cycles. The supercritical Hopf bifurcation occurs at $H_c^{(i)}\simeq 2.4$. The inset shows the COs' frequency ν⁽ⁱ⁾ of the inhibitory population as a function of H⁽ⁱ⁾. (B,C) Raster plots of the inhibitory network (blue dots) are calculated in correspondence with the stable fixed point at H⁽ⁱ⁾ = 0.0 (B) and with the limit cycle at H⁽ⁱ⁾ = +10.0 (C). Only the firing activity of 1,000 neurons is displayed. Parameters of the system: J⁽ⁱⁱ⁾ = 21.0, ${\bar{\eta }}^{(i)}=2.0$, Δ⁽ⁱ⁾ = 0.3, ${\tau }_m^{(i)}=10.0$ ms, τ_d = 10.0 ms, A = 0. The system size for the purely inhibitory network is N⁽ⁱ⁾ = 10, 000.

Therefore, the PING and ING set-ups considered here are ideal candidates to analyse the influence of θ-forcing on γ-oscillatory populations, which represents the main focus of this paper. In particular, the response of the system to the excitatory θ-forcing current (5) can be interpreted in terms of the bifurcation diagrams for the model in the absence of forcing shown, respectively, in Figure 2A for the PING set-up and in Figure 3A for the ING set-up. The interpretation is possible due to the fact that the response of the system to the sinusoidal current (5) can be considered as almost adiabatic, because the forcing frequencies ν_θ ∈ [1:10] Hz are definitely slower than those of the COs [ν^(e) and ν⁽ⁱ⁾], which lie in the γ-range.

4. Dynamics Under θ-Forcing

As a first step, we have verified that the reduced mean-field models are able to reproduce the macroscopic evolution of the spiking network in both considered set-ups, under the external forcing (5). In particular, we set the unforced systems in the asynchronous regime in proximity of a super-critical Hopf bifurcation, by choosing $H_0^{(e)}=1.3<H_c^{(e)}$ and $H_0^{(i)}=-5$ ($H_0^{(i)}=2.0<H_c^{(i)}$) and considered a forcing term with frequency ν_θ = 5 Hz and amplitude I₀ = 10 (I₀ = 9) for the PING (ING) set-up.

The comparisons, reported in Figures 4A,C, reveal a very good agreement in both set-ups between the network and the neural mass simulations, for the mean membrane voltages and the instantaneous firing rates. Furthermore, in both cases, we clearly observe COs, whose amplitudes are modulated by the amplitude of the θ-forcing term (5), suggesting that we are in the presence of a Phase-Amplitude Coupling (PAC) mechanism (Hyafil et al., 2015). The corresponding spectrograms shown in Figures 4B,D reveal that the frequencies of the COs are in the γ-range with the maximum power localized around 50–60 Hz. Moreover, the spectrograms indicate that the process is stationary and due to the external stimulation. The gamma oscillations repeat during each θ-cycle and they arrest when the external stimulation is stopped. The characteristics of these COs resemble θ-nested γ-oscillations reported in many experiments for neural systems in vitro under optogenetic stimulation (Akam et al., 2012; Pastoll et al., 2013; Butler et al., 2016, 2018) as well as in behaving animals (Chrobak and Buzsáki, 1998).

FIGURE 4

Figure 4. Theta-nested gamma oscillations (PING set-up). (A) From top to bottom: temporal traces of r^(e), v^(e), r⁽ⁱ⁾, v⁽ⁱ⁾, for the spiking network (red curves) and the neural mass model (black curves). I_θ, reported in the bottom panel in blue, is the external current (5). For the neural mass model the average rates and membrane potentials are solutions of Equation (6), while for the network they are calculated according to Equation (4). (B) Spectrogram of the mean membrane potential v^(e) (top) as a function of the external forcing (bottom). The amplitude of the forcing is I₀ = 10 and its frequency is ν_θ = 5 Hz. Parameters of the system: J^(ee) = 8, J^(ie) = J^(ei) = 10, J⁽ⁱⁱ⁾ = 0, $H_0^{(e)}=1.3$, $H_0^{(i)}=-5.0$, Δ^(e) = 1, ${\tau }_m^{(e)}=20$, Δ⁽ⁱ⁾ = 1, ${\tau }_m^{(i)}=10.0$, A = 0, network size N^(e) = N⁽ⁱ⁾ = 5, 000. The average firing rates are ${\bar{R}}^{(e)}\simeq 37$ Hz, ${\bar{R}}^{(i)}\simeq 36$ Hz. (ING set-up) (C) From top to bottom: temporal traces of r⁽ⁱ⁾, v⁽ⁱ⁾ where the line colors have the same meaning as in (A). For the neural mass model, average rates and membrane potentials are solutions of Equation (7). (D) Spectrogram of the mean membrane potential v⁽ⁱ⁾ (top) as a function of the external forcing (bottom). The amplitude of the forcing is I₀ = 9 and its frequency is ν_θ = 5 Hz. Parameters of the system: J⁽ⁱⁱ⁾ = 21.0, $H_0^{(i)}=2.0$, Δ⁽ⁱ⁾ = 0.3, ${\tau }_m^{(i)}=10.0$ ms, τ_d = 10.0 ms, A = 0, system size for the purely inhibitory network N⁽ⁱ⁾ = 10, 000. The corresponding average firing rate is ${\bar{R}}^{(i)}\simeq 28$ Hz.

4.1. Wavelet Analysis

To get a deeper insight into these dynamics we have estimated the continuous wavelet transform of the average membrane potential on each θ-cycle. As an example, we report in Figure 5 the wavelet spectrogram of the mean potential within a single θ-cycle for the previously examined PING (Figure 5A) and ING (Figure 5B) set-ups. Indeed, from the comparison of Figures 5A,B, we practically do not observe any difference: the system responds with COs in the range [40, 80] Hz and it exhibits alternating maxima and minima in the wavelet spectrogram as a function of the θ-phase. Similar results have been reported in Figure 4G in Butler et al. (2016) for the CA1-region of rat hippocampus under optogenetic sinusoidal θ-stimulation.

FIGURE 5

Figure 5. Wavelet analysis. Continuous wavelet transform over a single θ-cycle of the mean membrane potentials v^(e) and v⁽ⁱ⁾ appearing in the neural mass models for PING (A) and ING (B) set-up, respectively. This analysis allows for accurate automated detection and extraction of γ activity without the need for bandpass filtering. Parameters as in Figure 4.

Differences among the two cases appear when one considers the wavelet spectrograms averaged over many θ-periods: for the PING case the spectrogram remains unchanged, instead for the ING set-up the spectrogram smears out and it does not present anymore the clear oscillations reported in Figure 5B. This difference indicates that, in the PING case, the observed pattern repeats exactly over each cycle: γ-oscillations and θ-oscillations are perfectly phase locked. This is not the case for the ING set-up: although the PAC patterns appear quite similar in successive cycles, as shown in Figure 4C, indeed they do not repeat exactly. From the point of view of non-linear dynamics, the PING case would correspond to a perfectly periodic case, while the other case could be quasi-periodic or even chaotic. Therefore, we can observe PAC with an associated phase locking, but also in the absence of phase locking.

Furthermore, according to the data shown in Figure 5, this can also represent an example of PFC, since COs with frequencies ≃40 Hz occur at small and large θ-phases, while in the middle range π/2 < θ <3π/2 one observes similar oscillations with F ≃ 60 Hz.

For what concerns the wavelet analysis obtained from optogenetic experiments and shown in Figure 4G in Butler et al. (2016), we should stress two important aspects: (i) the wavelet spectrogram, averaged over several θ cycles (namely 30), displays clear correlations among the θ-phase and the γ-oscillations; (ii) the spectrogram is highly asymmetric indicating that γ-oscillations emerge in proximity of θ-phase ≃π and disappear ≃3/2π. The former aspect reveals that θ and γ oscillations were perfectly locked in the experiment, while the latter suggests that the bifurcation associated with the emergence of COs in the experiment is probably hysteretic. This would explain the asymmetry that we do not observe here for super-critical Hopf bifurcations in Figure 5, but that emerges for sub-critical Hopf bifurcations, as discussed in Appendix A.

4.2. Phase-Amplitude Locked and Unlocked States

To better examine the dynamical regimes emerging in our set-ups, we have first estimated the maximal Lyapunov exponent λ₁ associated with the neural mass models, for the same parameters considered in Figure 4, over a wide range of forcing amplitudes, that is 0 ≤ I₀ ≤ 20. From the results reported in Figures 6A,B, it is clear that λ₁ is almost always zero, apart from some limited intervals where it is negative and a few values of I₀ for the ING set-up, where it can be even positive. This means that the dynamics is usually quasi-periodic, apart from some Arnold tongues where there is perfect locking between the external forcing and the forced system.

FIGURE 6

Figure 6. Maximal Lyapunov exponent λ₁ estimated for the neural mass models as a function of the forcing amplitude I₀, for the PING (A) and ING (B) set-ups. In both cases the system is subject to a forcing frequency ν_θ = 5 Hz. Insets in (A,B) report the instantaneous firing rate r^(e)(t) (r⁽ⁱ⁾(t)) vs. time for the PING (ING) set-up, respectively. The three cases shown are representative of the states identified by circles in the main panels. The color code is the same, i.e., the color used in the inset identifies the corresponding circle in the main panel. The black continuous lines in the inset correspond to I_θ in arbitrary units. Parameters are the same as in Figure 4.

We notice that for small amplitudes the forcing entrains the system in a 1:1 periodic locking, therefore the instantaneous firing rate displays one peak for each θ-period with the same frequency as the forcing ν_θ. This locking is present in a wider region in the ING case (namely, I₀ < 1.70) with respect to the PING set-up (namely, I₀ < 0.40). More interesting locking regimes, where the forced populations oscillate in the γ-range, emerge at larger I₀. These locking regimes can be considered as θ-nested γ-oscillations; most of them are of the type m:1, with m ∈ [5:10], which means that, for each θ-period, the firing rate of the forced populations has m maxima (for specific examples see the insets of Figures 6A,B). In extremely narrow parameter intervals other, more complex, kinds of locking of the type m:n emerge, where exactly m maxima in the population activity appear for every n θ-oscillations. In the examined cases we have identified locked patterns with n up to four. Moreover, for the ING case, we have even observed a chaotic region (see Figure 6B), which emerges at quite large forcing amplitude I₀ ≃ 19. On the basis of our analysis we cannot exclude that chaos could emerge also in the PING set-up, for sufficiently strong forcing.

Let us now focus on the m:1 perfectly locked states with m > 1, which are worth investigating due to their relevance for θ-γ mixed oscillations, as well as to their relative large frequency of occurrence with respect to more complex m:n locked states. In particular, we have examined the response of the system to different forcing amplitudes I₀ ∈ [0:20] and frequencies ν_θ ∈ [1:10] Hz. The m:1 locked oscillations are reported in Figures 7A,B and characterized by the number m of oscillations displayed within a single θ-cycle.

FIGURE 7

Figure 7. Phase locked m:1 states Locked states for the neural mass models are displayed in (A,B) for the PING and ING set-ups, respectively. The color code identifies the locked states according to the value of m, from 3 to 15. (C,D) Minimal (red circles), average (blue circles), and maximal (black circles) frequencies of the COs as a function of the forcing amplitude I₀ for PING (C) and ING (D) set-ups. These values are obtained by considering all possible m:1 locked states corresponding to the examined I₀. The frequencies ν^(e) (ν⁽ⁱ⁾) (green solid lines) of the COs obtained from the bifurcation analysis in the adiabatic set-up are reported as a function of $H^{(e)}-H_0^{(e)}$ ($H^{(i)}-H_0^{(i)}$) for the PING (ING). Parameters are the same as in Figure 4.

These locked states appear only for ν_θ > 2−3 Hz. Moreover, the states with equal m are arranged in stripes in the (ν_θ, I₀)-plane. Locked states in the PING configuration occur in separated stripes whose order m increases for increasing I₀; in particular, states with 3 ≤ m ≤ 10 are clearly identifiable. In the ING set-up, for sufficiently large ν_θ and I₀, we have a continuum of locked states, thus indicating that, for the ING set-up, phase locking to the forcing frequency is easier to achieve. In this case the order of occurrence of m-order states is not clearly related to the forcing amplitude; however locked states with order m and 2m are often nested within each other as shown in Figure 7B.

To examine which frequencies are excited in these states we have measured for each amplitude I₀ the minimal, the maximal and the average frequency of the COs associated with m:1 locked states over the whole range of examined forcing frequencies ν_θ. These frequencies are reported in Figures 7C,D. The analysis clearly reveals that the minimal CO frequency is essentially independent from I₀ and its value is around 20 Hz, while the maximal and the average grow with I₀. However all these frequencies stay within the γ-range for the examined forcing amplitudes.

To better understand the mechanism underlying the emergence of θ-nested γ oscillations, we have reported in Figures 7C,D the COs frequencies ν^(e) (ν⁽ⁱ⁾) (green solid lines) obtained from the adiabatic bifurcation analysis of the neural mass models (these frequencies are also shown in the insets of Figures 2A, 3A). The very good agreement between ν^(e) and ν⁽ⁱ⁾ and the maximal frequency measured for the locked states suggests that the nested COs are induced by the crossing of the super-critical Hopf bifurcation during the periodic stimulation. In particular, during forcing, the maximal achievable γ-frequency is the one corresponding to the maximal stimulation current $I_0+H_0^{(e)}$ ($I_0+H_0^{(i)}$) for the unforced PING (ING) set-ups. Furthermore, under sinusoidal forcing, the system spends a longer time in proximity of the maximal stimulation value, since it is a turning point. This explains why this frequency is always present in the response of the driven system for the considered locked states.

4.3. Comparison with Experimental Findings

In a series of recent optogenetic experiments on the mouse enthorinal-hippocampal system, clear evidence has been reported that phase-amplitude coupled γ-rhythms can be generated locally in brain slices ex vivo in the CA1-region, as well as in the CA3 and MEC, under sinusoidal θ stimulations (Akam et al., 2012; Pastoll et al., 2013; Butler et al., 2016, 2018). In particular, in Butler et al. (2018) the authors reported evidence that, for all the regions CA1, CA3 and MEC, the generation of the γ-rhythms, under θ-rhythmic activation of pyramidal neurons, is due to a PING mechanism.

However, due to the fact that pyramidal neurons are directly activated during experiments, their result cannot exclude that tonic activation of interneurons contributes to θ–γ oscillations in vivo. Furthermore, in Pastoll et al. (2013) the authors affirm that θ-nested γ-oscillations due to the optogenetic θ-frequency drive, are generated, in MEC, by local feedback inhibition without recurrent excitation, therefore by a ING mechanism. In this section, we try to reproduce some of the analyses reported in these experimental studies by employing both the PING and ING set-ups, in order to understand if these two set-ups give rise to different dynamical behaviors.

By following the analysis performed in Butler et al. (2016, 2018), we have considered the response of the two set-ups to forcing of different frequencies ν_θ and amplitudes I₀. The results reported in Figure 8 reveal that the phenomenon of PAC is present for all the considered frequencies ν_θ ∈ [1, 10] Hz and amplitudes I₀ ∈ [1, 20] in both set-ups. Moreover, analogously to what was reported in Butler et al. (2016, 2018), the amplitude of the γ-oscillations increases proportionally to I₀, while the number of nested oscillations in each cycle increases for decreasing ν_θ. On the basis of this comparison, the forced PING and ING set-ups display essentially the same dynamics.

FIGURE 8

Figure 8. Theta-nested gamma COs for PING (A,B) and ING set-up (C,D). Left column: dependence of the mean membrane potential of the excitatory (inhibitory) population v^(e) (v⁽ⁱ⁾) on the frequency ν_θ of the external forcing $I^{(e)}=I_{\theta }$ ($I^{(i)}=I_{\theta }$) with I₀ = 10 (I₀ = 9) for the PING (ING) set-up. The current profiles (blue lines) are displayed immediately below the corresponding membrane potential evolution. From top to bottom, the frequency ν_θ is 1, 5, and 10 Hz. Right column: dependence of the mean membrane potential v^(e) (v⁽ⁱ⁾) on the amplitude I₀ of the external current. Here the forcing frequency is kept constant at the value ν_θ = 5Hz. The amplitude is changed from 100% of maximum (top) to 20% of maximum (bottom) in 20% increments, the maximum being given by I₀ = 10. The data refer to the evolution of neural mass models, the parameters are the same as in Figure 4.

To get a more detailed information about the dynamics in the two set-ups, we will now consider the features of the power spectra $P_S^{(e)}$ ($P_S^{(i)}$) of the mean excitatory (inhibitory) potential for the PING (ING) set-up. These features are obtained for different forcing amplitudes and frequencies, somehow similar to the analysis performed for the power spectra of the Local Field Potential (LFP) in Butler et al. (2016, 2018).

Let us first consider, as an example of the obtained power spectra, the case corresponding to the PING set-up with a forcing characterized by ν_θ = 5 Hz and amplitude I₀ = 10, shown in Figure 9A. In the spectrum we observe very well defined spectral lines located at frequencies which can be obtained as a linear combination of the forcing frequency ν_θ = 5 Hz and of the response frequency F_r = 45 Hz. In particular F_r is associated with the main peak and should correspond to the intrinsic frequency of the forced system. In the present case, the adiabatic bifurcation diagram reported in Figure 2A tells us that the maximal achievable frequency is ${\nu }_{max}^{(e)}\simeq 49.3$ Hz, corresponding to $H^{(e)}=I_0+H_0^{(e)}=11.3$. Indeed $F_r<{\nu }_{max}^{(e)}$ due to the interaction with the forcing current that eventually induces a locking phenomenon at a frequency that is exactly a multiple of ν_θ, as it happens in the present case. However, in general, a spectrum as the one shown in Figure 9A, is the emblem of a quasi-periodic motion characterized by two incommensurate frequencies. This can be easily observable in most cases in our system, where ν_θ and F_r are usually incommensurate.

FIGURE 9

Figure 9. Power spectra for the PING set-up. Spectra $P_S^{(e)}$ of the mean membrane potential v^(e) of the neural mass model estimated when the excitatory population is subject to an external drive with frequency ν_θ = 5 Hz and amplitude I₀ = 10, in the absence of noise (A) and for additive noise with amplitude A = 1.4 (B). The data refer to the evolution of neural mass models, the parameters are as in Figure 4.

The spectra obtained from optogenetic stimulation, reported in Butler et al. (2016, 2018), do not resemble the one shown in Figure 9A; indeed they present only two peaks: one corresponding to the stimulation frequency and one, quite broad, associated with the γ-oscillations. We can expect that the difference is due to the multiple noise sources that are always present in an experimental analysis (in particular for neurophysiological data), but that are absent in our model. Indeed, by considering the neural mass model for the PING set-up with additive noise on the membrane potentials of suitable amplitude, that is A = 1.4, we get a power spectrum resembling the experimental one, as shown in Figure 9B. The presence of noise induces the merging of the principal peaks in a unique broad one and the shift of the position of the main peak toward some larger values (F_r = 54 Hz in the present case) with respect to the fully deterministic case.

Let us now consider the power spectra obtained for different forcing frequencies ν_θ ∈ [1:10] Hz in the θ-range, in case of fixed forcing amplitude and in the absence of noise. The position of the main and auxiliary peaks are shown in Figure 10A (Figure 10C) for the PING (ING) set-up and compared with the experimental results (red circles) obtained for the CA1 region of the hippocampus in Butler et al. (2016). It is clear that, for both set-ups, the position of the main peak F_r (green squares) has a value ≃50 Hz and it does not show any clear dependence on ν_θ. This is in contrast with the experimental data, which reveal an increase proportional to ν_θ from 49 to 60 Hz. The same trend is displayed in our simulation from the subsidiary peak located at F_r + ν_θ (black stars), showing an increase with ν_θ.

FIGURE 10

Figure 10. Power spectra features (PING set-up). (A) Frequencies of the peaks of the power spectrum $P_S^{(e)}$ as a function of the stimulation frequency ν_θ. Green squares correspond to the main peak frequency F_r, while the black stars to F_r + ν_θ and the blue diamonds to F_r − ν_θ. The red circles are the experimental data extrapolated from Figure 4C of Butler et al. (2016). The amplitude of the forcing is I₀ = 10. (B) Normalized power of the γ oscillations $P_{\gamma }/P_{\gamma }^{max}$ associated with the signal v^(e) as a function of the amplitude stimulation, where we set $I_0^{max}=20$ and the frequency of stimulation at ν_θ = 5 Hz. In the inset we report the same quantity as a function of the frequency stimulation ν_θ for I₀ = 10. The black stars correspond to our simulations, while the red circles to experimental data extrapolated from Figure 4E (Figure 4B for the inset) of Butler et al. (2016) (filled circles) and from Figure 4C of Butler et al. (2018) (empty circles). The other parameters are as in Figure 4. (ING set-up) (C) Same as in (A) for the power spectrum $P_S^{(i)}$ with I₀ = 9. (D) Same as (B) for the signal v⁽ⁱ⁾ with $I_0^{max}=40$. The data refer to the evolution of neural mass models. For the inset we set I₀ = 9, other parameters as in Figure 4.

Let us now take into account the power of the γ oscillations P_γ as defined in section 2.2. As shown in the insets of Figures 10B,D, this quantity remains essentially constant for low frequencies (namely, for ν_θ ≤ 5 Hz in the PING and for ν_θ ≤ 7 Hz in the ING), while it drops to smaller values at larger frequencies. On the other hand, the experimental results (red circles) reveal a similar decrease at frequencies ν_θ>5 Hz, but they also reveal an increase at low frequencies, not present in our numerical data, thus suggesting a sort of resonance at 5 Hz. For what concerns the dependence of P_γ on the forcing amplitude, we have fixed ν_θ = 5 Hz and varied I₀ in the range [4:10] ([8:20]) for the PING (ING) set-up. In both cases and analogously to experimental data, P_γ increases proportionally to I₀, see Figures 10B,D.

In both set-ups, our model is unable to reproduce, in the absence of noise and for fixed forcing amplitude I₀, the steady increase of F_r with ν_θ reported in the experiments for the mice CA1 in Butler et al. (2016). Therefore, in order to cope with this problem, we will now investigate how a similar trend can emerge in our data. In particular, in the remaining part of the paper we consider noisy dynamics, to have a better match with experiments where is unavoidable. In Figure 11A we report, for the PING set-up, the estimated power spectra for different noise levels, under constant external sinusoidal forcing. The effect of noise is to render the spectrum more flat and to shift the position of the peak in the γ-range toward higher frequencies. As shown in the inset of Figure 11A, the frequency F_r is almost insensitive to the noise up to amplitudes A ≃ 1.0, then it increases steadily with A from ≃45 Hz to ≃62 Hz. The effect of varying the forcing amplitude I₀, for constant forcing frequency ν_θ = 5 Hz and noise amplitude A = 1.4, is shown in Figure 11B. In this case the amplitude increase of the forcing leads to more defined peaks in the γ-range and to an almost linear increase with I₀ of F_r, as reported in the inset. In the same inset we also have reported the results related to two optogenetic experiments for the CA1-region of the mice hippocampus. In particular, the data-sets refer to two successive experiments performed by the same group: namely, red filled circles refer to Butler et al. (2016) and red open circles to Butler et al. (2018). While in one experiment (red open circles) a constant increase of F_r with the forcing amplitude is observable from 60 to 70 Hz, in the other one (red filled circles) the frequency initially increases with I₀ and then decreases with it. As a matter of fact in the latter case, F_r remains around 45–50 Hz for a variation of I₀ from 40 to 100 % of the maximal amplitude $I_0^{max}$. From the comparison with our results, we can affirm that our data reproduce the correct range of frequencies in both experiments and also the dependence on the forcing amplitude for $I_0/I_0^{max}\ge 60\%$ reported in Butler et al. (2018). The decrease of F_r for $I_0/I_0^{max}$ larger than the 50% reported in Butler et al. (2016) is inconsistent with our data, but also with the experimental results of the same group published in Butler et al. (2018).

FIGURE 11

Figure 11. Power spectra dependency on noise and forcing amplitudes (PING set-up). Power spectra $P_S^{(e)}$ for different noise level A (A) and different amplitude of the external input I₀ (B), for a fixed forcing frequency ν_θ = 5 Hz. In the insets are reported the frequencies F_r of the main peak as a function of the noise level (A) and of the amplitude of the external drive I₀ (B). In the inset of (B) are also reported experimental data extracted from Figure 4F of Butler et al. (2016) (filled red circles) and from Figure 4D of Butler et al. (2018) (open red circles). The curves in (A) are obtained by varying the noise amplitude A ∈ [0.9:3.0] with a step of 0.3, while keeping I₀ = 10 fixed. On the other hand the curves in (B) refer to different forcing amplitudes 2 ≤ I₀ ≤ 20, varied in steps of 0.2, with fixed noise amplitude A = 1.4. The other parameters are as in Figure 4. Data have been obtained by the integration of neural mass models.

From this last analysis we have understood that, for constant forcing frequency, the γ-peak shifts toward higher frequencies by increasing the forcing amplitude or the noise level, i.e., by increasing the stimulation power.

Therefore, to obtain an increase of F_r with the forcing frequency ν_θ, analogously to the results reported in Butler et al. (2016) (and displayed as filled red circles in Figures 10A,C), we perform numerical experiments where ν_θ increases together with A or I₀. The simplest protocol is to assume that A (I₀) will increase linearly with ν_θ. The results obtained for the PING (ING) set-up are reported in Figure 12A (Figure 12B). As evident from the figures, in both set-ups and for both protocols we obtain results in reasonable agreement with the experiments. In the present framework, we have also analyzed the dependence of the γ-power P_γ on ν_θ. In particular, this quantity increases almost linearly with the forcing frequency, at variance with the experimental results in Butler et al. (2016) which revealed a sort of resonance with an associated maximal γ-power around ν_θ = 5 Hz (the experimental data are displayed as red circles in the insets of Figures 10B,D).

FIGURE 12

Figure 12. Influence of the theta frequency on the gamma oscillations. Frequency F_r of the main peak of the power spectrum $P_S^{(e)}$ vs. ν_θ for the PING (A) and ING (B) set-ups. Red filled circles represent the experimental data extrapolated from Figure 4C in Butler et al. (2016). Black stars (magenta triangles) refer to numerical data obtained by varying linearly the noise amplitude A (the forcing amplitude I₀) as a function of ν_θ and maintaining the forcing amplitude I₀ (the noise amplitude A) constant. The data shown as black stars for the PING (ING) set-up in (A) (panel B) are obtained by adding white noise to the evolution of the mean membrane potentials and by varying linearly its amplitude in the interval A ∈ [1.4:2.9] as a function of ν_θ with I₀ = 10 (I₀ = 9). The magenta triangles refer to data obtained by keeping fixed the noise amplitude at the value A = 1.4 and by varying linearly with ν_θ the forcing amplitude I₀ in the range [9.5:18] ([8:14]) for the PING (ING) set-up in (a) (panel B). Other parameters for as in Figure 4. Data are obtained from neural mass models integration.

5. Discussion and Conclusions

In this paper we have analyzed the dynamics of a new class of neural mass models arranged in two different set-ups: an excitatory-inhibitory network (or PING set-up) and a purely inhibitory network (or ING set-up). These neural mass models are extremely relevant to mimick neural dynamics for two reasons. On one side, because they are not derived heuristically, since they reproduce exactly the dynamics of excitatory and inhibitory networks of spiking neurons for any degree of synchronization (Montbrió et al., 2015; Devalle et al., 2017; Ceni et al., 2019). On another side, these neural masses reproduce the macroscopic dynamics of quadratic integrate-and-fire neurons, which are normal forms of class I neurons, therefore they are expected to represent the dynamics of this large class of neurons (Ermentrout and Kopell, 1986).

In this present work we have shown that θ-nested γ oscillations can emerge both in the PING and ING set-up under an external excitatory θ-drive whenever the system, in the absence of forcing, is in a regime of asynchronous dynamics, but in proximity of a Hopf bifurcation toward collective γ oscillations. The external forcing drives the system across the bifurcation inside the oscillatory regime, thus leading to the emergence of γ oscillations. The amplitude of these collective oscillations is related to the distance from the bifurcation point, therefore it depends on the phase of the θ-forcing term. These nested oscillations can arise in proximity of a super-critical and also a sub-critical Hopf bifurcations. As shown in Appendix A, in the latter case the amplitudes are no more symmetric with respect to the maximum value of the theta stimulation, analogously to the experimental findings reported in Butler et al. (2016).

Equivalent results have been reported for an excitatory-inhibitory network with a recurrent coupling among the excitatory neurons, by considering the Wilson-Cowan rate model (Onslow et al., 2014). However, at variance with our neural mass model, the Wilson-Cowan model fails to reproduce the emergence of γ-oscillations, displayed by the corresponding spiking networks, in several other set-ups. In particular, the Wilson-Cowan model is unable to display COs for purely inhibitory populations (the ING set-up), without the addition of a delay in the IPSPs transmission, delay that is not required in the network model. Moreover, the Wilson-Cowan model is unable to display COs even for excitatory-inhibitory coupled populations in the absence of a recurrent excitation (Onslow et al., 2014; Devalle et al., 2017). As shown in Appendix B, the considered neural mass model in the PING set-up displays clear θ-nested γ-oscillations in the absence of any recurrent coupling or with recurrent couplings only among the inhibitory neurons.

Furthermore, we have identified two different types of phase amplitude couplings. One characterized by a perfect locking between θ and γ-rhythms, corresponding to an overall periodic behavior dictated by the slow forcing. The other one where the locking is imperfect and the dynamics is quasi-periodic or even chaotic. The perfectly locked θ-nested γ oscillations display in turn two types of CFC: phase-phase and phase amplitude coupling (Hyafil et al., 2015). These states arise for ν_θ larger than 2–3 Hz and for sufficiently large forcing amplitudes. From the results reported in Butler et al. (2016) for the CA1-region of the hippocampus under sinusoidal forcing in vitro, it is evident that perfectly phase locked PACs have been observed in each single slice. However, in vivo this perfect phase-phase locking cannot be expected, see the detailed discussion of phase-phase coupling reported in Scheffer-Teixeira and Tort (2016), where the authors clarify that phase locking is indeed observable, but only over a limited number of successive θ-cycles. Therefore, PAC with an underlying chaotic (or noisy) dynamics is the scenario usually expected in behaving animals.

From our analysis it also emerges that locked states are more frequent in the ING set-up. The purely inhibitory population is more easily entrained by the forcing with respect to the coupled excitatory-inhibitory population system, where the forcing is applied to the excitatory population. This result is somehow in agreement with recent findings based on the analysis of phase response curves, which suggest that stimulating the inhibitory population facilitates the entrainment of the gamma-bands with an almost resonant frequency (Akao et al., 2018; Dumont and Gutkin, 2019). However, these analyses do not consider θ-γ entrainment: this will be a subject of future studies based on exact macroscopic phase response curves (Dumont et al., 2017; Dumont and Gutkin, 2019).

Our modelization of the PAC mechanism induced by an external θ-forcing is able to reproduce several experimental features reported for optogenetic experiments concerning the region CA1, CA3 of the hippocampus, as well as MEC (Akam et al., 2012; Pastoll et al., 2013; Butler et al., 2016, 2018). In agreement with the experiments, we observe nested γ COs for forcing frequencies in the range [1:10] Hz, whose amplitude grows proportionally to the forcing one. Furthermore, the γ-power and the frequency of the γ peak increase almost linearly with the forcing amplitude, i.e., with the input θ-power. Moreover these findings are consistent with recent results for behaving rats, where it has been shown that hippocampal γ-frequency and the associated power increase proportionally to the animal speed (Ahmed and Mehta, 2012; Sheremet et al., 2019). In addition, in Richard et al. (2013), the authors have clearly demonstrated that the hippocampal θ-power and the mouse speed are positively correlated. This proportionality between the θ-power and the mouse speed has been recently employed to develop a computational model able to successfully reproduce CA1 network activity (Haimerl et al., 2019).

However, the neural mass model in all the examined PING and ING set-ups is unable to reproduce the increase in frequency of the γ-power peak with ν_θ reported in Butler et al. (2016). Indeed, such effect was expected by the observation that during movement, both the frequencies of hippocampal θ oscillations (Sławińska and Kasicki, 1998) and γ oscillations (Ahmed and Mehta, 2012) increase with the running speed of the animal. However, the variation of the γ frequency reported in Ahmed and Mehta (2012) for behaving animals amounts to 40–60 Hz, while in the optogenetic experiment by Butler et al. (2016), the increase was limited to ≃10 Hz. In order to get a similar increase in the neural mass model, we have been obliged to assume that the stimulation power (namely, the noise or the forcing amplitude) increases proportionally to ν_θ. On one side, further experiments are required to clarify if, during optogenetic experiments, the forcing (or noise amplitude) affecting the neural dynamics is indeed dependent on ν_θ. This could be due to a reinforcement of the synaptic strengths for increasing forcing frequencies, or to the fact that higher θ frequencies can favor neural discharges in regions different from CA1, thus being assimilated to external noise. On another side it should be analyzed if other bifurcation mechanisms, beside the Hopf one, here considered, can give rise to such a dependence of γ power on θ forcing.

Finally, experiments on behaving rodents report clear evidence that θ-power and ν_θ, as well as the power of the γ-peak and the corresponding frequency, increase all proportionally to the animal speed (Sławińska and Kasicki, 1998; Ahmed and Mehta, 2012; Richard et al., 2013; Sheremet et al., 2019). Furthermore, in Sheremet et al. (2019) the authors report evidence of the increase of the phase-amplitude coupling with the speed. This scenario is consistent with the results reported in our analysis, where we have shown that an increase of ν_θ and of the stimulation power leads to an increase of P_γ and of the frequency of the γ peak as well as of the PAC.

Data Availability Statement

The datasets generated for this study are available on request to the corresponding author.

Author Contributions

MS, HB, and SO performed the simulations and data analysis. SO and AT were responsible for the state-of-the-art review and the paper write-up. All the authors conceived and planned the research.

Funding

AT received financial support by the Excellence Initiative I-Site Paris Seine (Grant No ANR-16-IDEX-008) (together with HB and MS), by the Labex MME-DII (Grant No ANR-11-LBX-0023-01) (together with SO) and by the ANR Project ERMUNDY (Grant No ANR-18-CE37-0014), all part of the French programme Investissements d'Avenir.

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

This work has been initiated during the Advanced Study Group 2016/17 From Microscopic to Collective Dynamics in Neural Circuits at the Max-Planck Institute for Physics of Complex System in Dresden (Germany), where SO and AT had extremely useful interactions with E. Montbrió and H. P. Robinson. We also acknowledge fruitful discussions with D. Angulo-Garcia, F. Devalle, G. Dumont, M. di Volo, A. Pikovsky and we thank T. Kreuz for a carefully reading of the final manuscript. We would also like to thank A.P. Maureer for pointing out recent experimental findings on the Hippocampus of behaving rodents.

Supplementary Material

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

References

Ahmed, O. J., and Mehta, M. R. (2012). Running speed alters the frequency of hippocampal gamma oscillations. J. Neurosci. 32, 7373–7383. doi: 10.1523/JNEUROSCI.5110-11.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Akam, T., Oren, I., Mantoan, L., Ferenczi, E., and Kullmann, D. M. (2012). Oscillatory dynamics in the hippocampus support dentate gyrus-ca3 coupling. Nat. Neurosci. 15:763. doi: 10.1038/nn.3081

PubMed Abstract | CrossRef Full Text | Google Scholar

Akao, A., Ogawa, Y., Jimbo, Y., Ermentrout, G. B., and Kotani, K. (2018). Relationship between the mechanisms of gamma rhythm generation and the magnitude of the macroscopic phase response function in a population of excitatory and inhibitory modified quadratic integrate-and-fire neurons. Phys. Rev. E 97:012209. doi: 10.1103/PhysRevE.97.012209

PubMed Abstract | CrossRef Full Text | Google Scholar

Bartos, M., Vida, I., and Jonas, P. (2007). Synaptic mechanisms of synchronized gamma oscillations in inhibitory interneuron networks. Nat. Rev. Neurosci. 8:45. doi: 10.1038/nrn2044

PubMed Abstract | CrossRef Full Text | Google Scholar

Belluscio, M. A., Mizuseki, K., Schmidt, R., Kempter, R., and Buzsáki, G. (2012). Cross-frequency phase-phase coupling between theta and gamma oscillations in the hippocampus. J. Neurosci. 32, 423–435. doi: 10.1523/JNEUROSCI.4122-11.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Benettin, G., Galgani, L., Giorgilli, A., and Strelcyn, J.-M. (1980). Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them. Part 1: Theory. Meccanica 15, 9–20. doi: 10.1007/BF02128236

CrossRef Full Text | Google Scholar

Bragin, A., Jandó, G., Nádasdy, Z., Hetke, J., Wise, K., and Buzsáki, G. (1995). Gamma (40-100 hz) oscillation in the hippocampus of the behaving rat. J. Neurosci. 15, 47–60. doi: 10.1523/JNEUROSCI.15-01-00047.1995

PubMed Abstract | CrossRef Full Text | Google Scholar

Brunel, N., and Hakim, V. (1999). Fast global oscillations in networks of integrate-and-fire neurons with low firing rates. Neural Comput. 11, 1621–1671. doi: 10.1162/089976699300016179

PubMed Abstract | CrossRef Full Text | Google Scholar

Butler, J. L., Hay, Y. A., and Paulsen, O. (2018). Comparison of three gamma oscillations in the mouse entorhinal-hippocampal system. Eur. J. Neurosci. 48, 2795–2806. doi: 10.1111/ejn.13831

PubMed Abstract | CrossRef Full Text | Google Scholar

Butler, J. L., Mendonça, P. R., Robinson, H. P., and Paulsen, O. (2016). Intrinsic cornu ammonis area 1 theta-nested gamma oscillations induced by optogenetic theta frequency stimulation. J. Neurosci. 36, 4155–4169. doi: 10.1523/JNEUROSCI.3150-15.2016

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsaki, G. (2006). Rhythms of the Brain. New York, NY: Oxford University Press. doi: 10.1093/acprof:oso/9780195301069.001.0001

CrossRef Full Text | Google Scholar

Buzsáki, G., and Wang, X.-J. (2012). Mechanisms of gamma oscillations. Annu. Rev. Neurosci. 35, 203–225. doi: 10.1146/annurev-neuro-062111-150444

PubMed Abstract | CrossRef Full Text | Google Scholar

Canolty, R. T., Edwards, E., Dalal, S. S., Soltani, M., Nagarajan, S. S., Kirsch, H. E., et al. (2006). High gamma power is phase-locked to theta oscillations in human neocortex. Science 313, 1626–1628. doi: 10.1126/science.1128115

PubMed Abstract | CrossRef Full Text | Google Scholar

Canolty, R. T., and Knight, R. T. (2010). The functional role of cross-frequency coupling. Trends Cogn. Sci. 14, 506–515. doi: 10.1016/j.tics.2010.09.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Ceni, A., Olmi, S., Torcini, A., and Angulo-Garcia, D. (2019). Cross frequency coupling in next generation inhibitory neural mass models. arXiv preprint arXiv:1908.07954. doi: 10.1101/745828

CrossRef Full Text | Google Scholar

Chehelcheraghi, M., van Leeuwen, C., Steur, E., and Nakatani, C. (2017). A neural mass model of cross frequency coupling. PLoS ONE 12:e0173776. doi: 10.1371/journal.pone.0173776

PubMed Abstract | CrossRef Full Text | Google Scholar

Chrobak, J., and Buzsáki, G. (1998). Gamma oscillations in the entorhinal cortex of the freely behaving rat. J. Neurosci. 18, 388–398. doi: 10.1523/JNEUROSCI.18-01-00388.1998

PubMed Abstract | CrossRef Full Text | Google Scholar

Colgin, L. L. (2015). Theta-gamma coupling in the entorhinal-hippocampal system. Curr. Opin. Neurobiol. 31, 45–50. doi: 10.1016/j.conb.2014.08.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Colgin, L. L., Denninger, T., Fyhn, M., Hafting, T., Bonnevie, T., Jensen, O., et al. (2009). Frequency of gamma oscillations routes flow of information in the hippocampus. Nature 462:353. doi: 10.1038/nature08573

PubMed Abstract | CrossRef Full Text | Google Scholar

Coombes, S., and Byrne, Á. (2019). “Next generation neural mass models,” in Nonlinear Dynamics in Computational Neuroscience, eds F. Corinto and A. Torcini (Cham: Springer), 1–16. doi: 10.1007/978-3-319-71048-8_1

PubMed Abstract | CrossRef Full Text | Google Scholar

Devalle, F., Roxin, A., and Montbrió, E. (2017). Firing rate equations require a spike synchrony mechanism to correctly describe fast oscillations in inhibitory networks. PLoS Comput. Biol. 13:e1005881. doi: 10.1371/journal.pcbi.1005881

PubMed Abstract | CrossRef Full Text | Google Scholar

di Volo, M., and Torcini, A. (2018). Transition from asynchronous to oscillatory dynamics in balanced spiking networks with instantaneous synapses. Phys. Rev. Lett. 121:128301. doi: 10.1103/PhysRevLett.121.128301

PubMed Abstract | CrossRef Full Text | Google Scholar

Dumont, G., Ermentrout, G. B., and Gutkin, B. (2017). Macroscopic phase-resetting curves for spiking neural networks. Phys. Rev. E 96:042311. doi: 10.1103/PhysRevE.96.042311

PubMed Abstract | CrossRef Full Text | Google Scholar

Dumont, G., and Gutkin, B. (2019). Macroscopic phase resetting-curves determine oscillatory coherence and signal transfer in inter-coupled neural circuits. PLoS Comput. Biol. 15:e1007019. doi: 10.1371/journal.pcbi.1007019

PubMed Abstract | CrossRef Full Text | Google Scholar

Ermentrout, G. B., and Kopell, N. (1986). Parabolic bursting in an excitable system coupled with a slow oscillation. SIAM J. Appl. Math. 46, 233–253. doi: 10.1137/0146017

CrossRef Full Text | Google Scholar

Fries, P., Nikolić, D., and Singer, W. (2007). The gamma cycle. Trends Neurosci. 30, 309–316. doi: 10.1016/j.tins.2007.05.005

CrossRef Full Text | Google Scholar

Gerstner, W., Kistler, W. M., Naud, R., and Paninski, L. (2014). Neuronal Dynamics: From Single Neurons to Networks and Models of Cognition. Cambridge: Cambridge University Press. doi: 10.1017/CBO9781107447615

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldobin, D. S., Tyulkina, I. V., Klimenko, L. S., and Pikovsky, A. (2018). Collective mode reductions for populations of coupled noisy oscillators. Chaos 28:101101. doi: 10.1063/1.5053576

PubMed Abstract | CrossRef Full Text | Google Scholar

Govaerts, W., Kuznetsov, Y. A., and Sautois, B. (2006). Matcont. Scholarpedia 1:1375. doi: 10.4249/scholarpedia.1375

CrossRef Full Text | Google Scholar

Haimerl, C., Angulo-Garcia, D., Villette, V., Reichinnek, S., Torcini, A., Cossart, R., and Malvache, A. (2019). Internal representation of hippocampal neuronal population spans a time-distance continuum. Proc. Natl. Acad. Sci. U.S.A. 116, 7477–7482. doi: 10.1073/pnas.1718518116

PubMed Abstract | CrossRef Full Text | Google Scholar

Hyafil, A., Giraud, A.-L., Fontolan, L., and Gutkin, B. (2015). Neural cross-frequency coupling: connecting architectures, mechanisms, and functions. Trends Neurosci. 38, 725–740. doi: 10.1016/j.tins.2015.09.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Jensen, O., and Colgin, L. L. (2007). Cross-frequency coupling between neuronal oscillations. Trends Cogn. Sci. 11, 267–269. doi: 10.1016/j.tics.2007.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Kilpatrick, Z. P. (2015). “Wilson-Cowan model,” in Encyclopedia of Computational Neuroscience, eds D. Jaeger and R. Jung (New York, NY: Springer), 3159–3163. doi: 10.1007/978-1-4614-6675-8_80

CrossRef Full Text | Google Scholar

Laing, C. R. (2014). Derivation of a neural field model from a network of theta neurons. Phys. Rev. E 90:010901. doi: 10.1103/PhysRevE.90.010901

PubMed Abstract | CrossRef Full Text | Google Scholar

Laing, C. R. (2017). “Phase oscillator network models of brain dynamics,” in Computational Models of Brain and Behavior, 1st ed, ed A. Moustafa (Hoboken, NJ: WILEY Blackwell), 505–517. doi: 10.1002/9781119159193.ch37

CrossRef Full Text | Google Scholar

Lega, B., Burke, J., Jacobs, J., and Kahana, M. J. (2016). Slow-theta-to-gamma phase-amplitude coupling in human hippocampus supports the formation of new episodic memories. Cereb. Cortex 26, 268–278. doi: 10.1093/cercor/bhu232

PubMed Abstract | CrossRef Full Text | Google Scholar

Lisman, J. (2005). The theta/gamma discrete phase code occuring during the hippocampal phase precession may be a more general brain coding scheme. Hippocampus 15, 913–922. doi: 10.1002/hipo.20121

PubMed Abstract | CrossRef Full Text | Google Scholar

Lisman, J. E., and Jensen, O. (2013). The theta-gamma neural code. Neuron 77, 1002–1016. doi: 10.1016/j.neuron.2013.03.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Luke, T. B., Barreto, E., and So, P. (2013). Complete classification of the macroscopic behavior of a heterogeneous network of theta neurons. Neural Comput. 25, 3207–3234. doi: 10.1162/NECO_a_00525

PubMed Abstract | CrossRef Full Text | Google Scholar

Mazzoni, A., Brunel, N., Cavallari, S., Logothetis, N. K., and Panzeri, S. (2011). Cortical dynamics during naturalistic sensory stimulations: experiments and models. J. Physiol. Paris 105, 2–15. doi: 10.1016/j.jphysparis.2011.07.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Montbrió, E., Pazó, D., and Roxin, A. (2015). Macroscopic description for networks of spiking neurons. Phys. Rev. X 5:021028. doi: 10.1103/PhysRevX.5.021028

CrossRef Full Text | Google Scholar

Olmi, S., Livi, R., Politi, A., and Torcini, A. (2010). Collective oscillations in disordered neural networks. Phys. Rev. E 81:046119. doi: 10.1103/PhysRevE.81.046119

PubMed Abstract | CrossRef Full Text | Google Scholar

Onslow, A. C., Jones, M. W., and Bogacz, R. (2014). A canonical circuit for generating phase-amplitude coupling. PLoS ONE 9:e102591. doi: 10.1371/journal.pone.0102591

PubMed Abstract | CrossRef Full Text | Google Scholar

Pastoll, H., Solanka, L., van Rossum, M. C., and Nolan, M. F. (2013). Feedback inhibition enables theta-nested gamma oscillations and grid firing fields. Neuron 77, 141–154. doi: 10.1016/j.neuron.2012.11.032

PubMed Abstract | CrossRef Full Text | Google Scholar

Pernía-Andrade, A. J., and Jonas, P. (2014). Theta-gamma-modulated synaptic currents in hippocampal granule cells in vivo define a mechanism for network oscillations. Neuron 81, 140–152. doi: 10.1016/j.neuron.2013.09.046

PubMed Abstract | CrossRef Full Text | Google Scholar

Pikovsky, A., and Politi, A. (2016). Lyapunov Exponents: A Tool to Explore Complex Dynamics. Cambridge: Cambridge University Press. doi: 10.1017/CBO9781139343473

CrossRef Full Text | Google Scholar

Richard, G. R., Titiz, A., Tyler, A., Holmes, G. L., Scott, R. C., and Lenck-Santini, P.-P. (2013). Speed modulation of hippocampal theta frequency correlates with spatial memory performance. Hippocampus 23, 1269–1279. doi: 10.1002/hipo.22164

PubMed Abstract | CrossRef Full Text | Google Scholar

Scheffer-Teixeira, R., and Tort, A. B. (2016). On cross-frequency phase-phase coupling between theta and gamma oscillations in the hippocampus. Elife 5:e20515. doi: 10.7554/eLife.20515

PubMed Abstract | CrossRef Full Text | Google Scholar

Sheremet, A., Kennedy, J., Qin, Y., Zhou, Y., Lovett, S. D., Burke, S. N., et al. (2019). Theta-gamma cascades and running speed. J. Neurophysiol. 121, 444–458. doi: 10.1152/jn.00636.2018

PubMed Abstract | CrossRef Full Text | Google Scholar

Sławińska, U., and Kasicki, S. (1998). The frequency of rat's hippocampal theta rhythm is related to the speed of locomotion. Brain Res. 796, 327–331. doi: 10.1016/S0006-8993(98)00390-4

PubMed Abstract | CrossRef Full Text | Google Scholar

So, P., Luke, T. B., and Barreto, E. (2014). Networks of theta neurons with time-varying excitability: macroscopic chaos, multistability, and final-state uncertainty. Phys. D Nonlinear Phenomena 267, 16–26. doi: 10.1016/j.physd.2013.04.009

CrossRef Full Text | Google Scholar

Tass, P., Rosenblum, M. G., Weule, J., Kurths, J., Pikovsky, A., Volkmann, J., et al. (1998). Detection of n:m phase locking from noisy data: application to magnetoencephalography. Phys. Rev. Lett. 81, 3291–3294. doi: 10.1103/PhysRevLett.81.3291

CrossRef Full Text | Google Scholar

Tiesinga, P., and Sejnowski, T. J. (2009). Cortical enlightenment: are attentional gamma oscillations driven by ING or PING? Neuron 63, 727–732. doi: 10.1016/j.neuron.2009.09.009

CrossRef Full Text | Google Scholar

Tyulkina, I. V., Goldobin, D. S., Klimenko, L. S., and Pikovsky, A. (2018). Dynamics of noisy oscillator populations beyond the ott-antonsen ansatz. Phys. Rev. Lett. 120:264101. doi: 10.1103/PhysRevLett.120.264101

PubMed Abstract | CrossRef Full Text | Google Scholar

Uhlhaas, P. J., and Singer, W. (2006). Neural synchrony in brain disorders: relevance for cognitive dysfunctions and pathophysiology. Neuron 52, 155–168. doi: 10.1016/j.neuron.2006.09.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Whittingstall, K., and Logothetis, N. K. (2009). Frequency-band coupling in surface EEG reflects spiking activity in monkey visual cortex. Neuron 64, 281–289. doi: 10.1016/j.neuron.2009.08.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Whittington, M. A., Traub, R. D., and Jefferys, J. G. (1995). Synchronized oscillations in interneuron networks driven by metabotropic glutamate receptor activation. Nature 373, 612–615. doi: 10.1038/373612a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Williams, S., and Boksa, P. (2010). Gamma oscillations and schizophrenia. J. Psychiatry Neurosci. 35:75. doi: 10.1503/jpn.100021

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilson, H. R., and Cowan, J. D. (1972). Excitatory and inhibitory interactions in localized populations of model neurons. Biophys. J. 12, 1–24. doi: 10.1016/S0006-3495(72)86068-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Witte, H., Putsche, P., Hemmelmann, C., Schelenz, C., and Leistritz, L. (2008). Analysis and modeling of time-variant amplitude-frequency couplings of and between oscillations of EEG bursts. Biol. Cybernet. 99, 139–157. doi: 10.1007/s00422-008-0245-x

PubMed Abstract | CrossRef Full Text | Google Scholar