Theta-Nested Gamma Oscillations in Next Generation Neural Mass Models

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.

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

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:

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)

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:

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_{\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

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:

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

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

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 λ₁ > 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.

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

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

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.

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.

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.

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.

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.

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 ν_θ.

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

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

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