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.

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.

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 suprathreshold (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;Perní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 crossfrequency 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 amplitudeamplitude 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., 2016Butler et al., , 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(Butler et al., , 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 excitatoryinhibitory 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 setups, 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. (2016Butler et al. ( , 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 setup in Appendix A. Further network configurations ensuring the emergence of COs via PING mechanism are presented in Appendix B.

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.
The dynamics of the membrane potentials of the QIF neurons in the PING configuration is given by 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.
where the super-scripts e (i) denote the excitatory (inhibitory) population, τ (e) m = 20 ms (τ (i) m = 10 ms) is the excitatory (inhibitory) membrane time constant, η (l) k 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 (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 (l) m whenever the membrane potential of a generic neuron j reaches infinity, i.e., 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 (ii) = 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 (ii) = 0 and J (ee) = 0) and in the absence of both self-activation and inhibition (i.e., with For what concerns the purely inhibitory network, the membrane potential dynamics of the j-th neuron is ruled by the following equations: where τ (i) m = 10 ms. In this case the synaptic field s (i) (t) is the super-position of the exponential IPSPs p(t) = e −t/τ d /τ 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 η (l) i are randomly distributed according to a Lorentzian probability density function (PDF) where H (l) is the median and (l) is the half-width halfmaximum (HWHM) of the PDF. Therefore each population will be composed of neurons supra-(with η (l) j > 0) and sub-threshold (with η (l) j < 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) = (i) = 1, whereas varying H (e) and H (i) . For the ING set-up we fix (i) = 0.3 and analyze the dynamics by varying H (i) .
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: 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 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 θ (t) and I (i) (t) ≡ 0 and for the ING set-up I (i) (t) = I θ (t).

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 aṡ 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 ξ (i) 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:ṙ 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 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 (i) δv (i) . The dynamics of the tangent vector is ruled by the linearization of the Equation (6), namely For the ING set-up the tangent vector is three dimensional, δ = δr (i) , δv (i) , δs (i) , and its time evolution can be obtained by the linearization of Equation (7), which reads as 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 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 λ 1 > 0, while a periodic (quasi-periodic) dynamics will be characterized by λ 1 = 0 (λ 1 = λ 2 = 0) and a fixed point by λ 1 < 0. In a nonautonomous system in the presence of an external forcing, one Lyapunov exponent will be necessarily zero, therefore a periodic behavior corresponds to λ 1 < 0 and a quasi-periodic dynamics to λ 1 = 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 (e) S (F) (P (i) S (F)) of the mean excitatory (inhibitory) membrane potential v (e) (t) (v (i) (t)) for the PING (ING) setup. 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.

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 (i) ≡ 0] as a function of the medians H (e) and H (i) 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.

PING set-up
For the excitatory-inhibitory set-up, as already mentioned, we usually fix H (i) = −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 (e) c ≃ 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.
Moreover, it is important to investigate the bifurcation diagram of the system at fixed median excitatory drive by varying H (i) . The corresponding bifurcation diagram is displayed in Figure 2B for H (e) = 10. By increasing H (i) , COs emerge from the asynchronous state via a sub-critical Hopf bifurcation at H (i) c1 ≃ −8.4 and they disappear via a super-critical Hopf bifurcation at H (i) c2 ≃ 0.20. Since the first transition is hysteretical, COs disappear via a saddle-node of the limit cycles at a value we have the coexistence of a stable focus with a stable limit cycle. In summary, COs are clearly observable as long as H (i) 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 τ (i) m = 10 ms (Tiesinga and Sejnowski, 2009). On the other hand, the influence of H (i) 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, ν (i) ≡ ν (e) varies by few Hz (namely, from 42.8 to 46.9 Hz), when H (i) 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 τ (i) m . This can be explained by the fact that the excitatory stimulations should reach the inhibitory population within a time interval of (at most) ≃ τ (i) m 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 (i) has in general the effect to reduce T a ; this should be expected since for larger excitabilities [larger H (i) ], 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 (i) SN and for positive H (i) , where the effect is reversed and T a increases with H (i) . 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.

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 (i) , we observe a super-critical Hopf bifurcation at H (i) c ≃ 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 ν (i) ∈ [26 : 83] Hz. In particular, we observe an almost linear increase of ν (i) with H (i) . An example of the observed dynamics is shown in Figure 3B (Figure 3C) where the raster plot of the inhibitory network is calculated for 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 ν (i) ], which lie in the γ -range.

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 (e) 0 = 1.3 < H (e) c and H (i) 0 = −5 (H (i) 0 = 2.0 < H (i) c ) and considered a forcing term with frequency ν θ = 5 Hz and amplitude I 0 = 10 (I 0 = 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., 2016Butler et al., , 2018 as well as in behaving animals (Chrobak and Buzsáki, 1998).

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

set-ups. Indeed, from the comparison of
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.

Phase-Amplitude Locked and Unlocked States
To better examine the dynamical regimes emerging in our setups, we have first estimated the maximal Lyapunov exponent λ 1 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 0 ≤ 20. From the results reported in Figures 6A,B, it is clear that λ 1 is almost always zero, apart from some limited intervals where it is negative and a few values of I 0 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.
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 0 < 1.70) with respect to the PING set-up (namely, I 0 < 0.40). More interesting locking regimes, where the forced populations oscillate in the γ -range, emerge at larger I 0 . 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 0 ≃ 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 FIGURE 6 | Maximal Lyapunov exponent λ 1 estimated for the neural mass models as a function of the forcing amplitude I 0 , 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 (i) (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.
θ -γ 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 ∈ [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.
These locked states appear only for ν θ > 2 − 3 Hz. Moreover, the states with equal m are arranged in stripes in the (ν θ , I 0 )plane. Locked states in the PING configuration occur in separated stripes whose order m increases for increasing I 0 ; in particular, states with 3 ≤ m ≤ 10 are clearly identifiable. In the ING setup, for sufficiently large ν θ and I 0 , 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 0 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 0 and its value is around 20 Hz, while the maximal and the average grow with I 0 . 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) (ν (i) ) (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 ν (i) 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 (e) 0 (I 0 + H (i) 0 ) 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.

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., 2016Butler et al., , 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. (2016Butler et al. ( , 2018, we have considered the response of the two set-ups to forcing of different frequencies ν θ and amplitudes I 0 . 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 0 ∈ [1, 20] in both set-ups. Moreover, analogously to what was reported in Butler et al. (2016Butler et al. ( , 2018, the amplitude of the γ -oscillations increases proportionally to I 0 , 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.
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 (e) S (P (i) S ) 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. (2016Butler et al. ( , 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 0 = 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 ν (e) max ≃ 49.3 Hz, corresponding to H (e) = I 0 + H (e) 0 = 11.3. Indeed F r < ν (e) max 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. The spectra obtained from optogenetic stimulation, reported in Butler et al. (2016Butler et al. ( , 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. population v (e) (v (i) ) on the frequency ν θ of the external forcing I (e) = I θ (I (i) = I θ ) with I 0 = 10 (I 0 = 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 (i) ) on the amplitude I 0 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 0 = 10. The data refer to the evolution of neural mass models, the parameters are the same as in Figure 4.
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 ν θ .
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 In both set-ups, our model is unable to reproduce, in the absence of noise and for fixed forcing amplitude I 0 , 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 0 , 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 0 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 CA1region 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 0 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 0 from 40 to 100 % of the maximal amplitude I max 0 . 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 max 0 ≥ 60% reported in Butler et al. (2018). The decrease of F r for I 0 /I max 0 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).
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 0 . The simplest protocol is to assume that A (I 0 ) will increase linearly with ν θ . The results obtained for the PING (ING) setup 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).

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 S 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 0 = 10. (B) Normalized power of the γ oscillations P γ /P max γ associated with the signal v (e) as a function of the amplitude stimulation, where we set I max 0 = 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 0 = 10. The black stars correspond to our simulations, while the red circles to experimental data extrapolated from Figure 4E ( Figure 4B   S for different noise level A (A) and different amplitude of the external input I 0 (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 0 (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   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 excitatoryinhibitory 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., 2016Butler et al., , 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.