Theta-gamma coupling: nonlinearity as a universal cross-frequency coupling mechanism

The Cross Frequency Coupling (CFC) phenomenon is defined as a statistical correlation between characteristic parameters neural oscillations. This study demonstrates and analyzes the nonlinear mechanism of the CFC, with a focus on the coupling between slow and fast oscillations, as a model for theta-gamma coupling. We first discuss the usage of the spectrum/bispectrum CFC measure using experimental data. As a physical paradigm, we propose the concept of a Class II neural population at low activity: neurons fire intermittently, and the time spent in the subthreshold regime is much larger that the duration of an action potential. We verify the emergence of fast oscillations (gamma) using a direct numerical simulations (DNS) of a population of Hodgkin-Huxley neurons forced by a slow theta oscillation. To deconstruct the mechanism, we derive a mean field approximation based on a reduction of the Hodgkin-Huxley model to a two-equation leaky-integrate-and-fire (LIF) model. Under theta forcing, mean field model generates gamma oscillations; the solutions exhibit spectrum/bispectrum CFC patterns that agree qualitatively with both the DNS model and experimental data. For the theta-gamma coupling problem, the mean field model may be linearized using an asymptotic expansion. The analytical solution of the linear system describes theta-gamma interaction as a gamma stabilization/destabilization cycle and provides explicit expressions of the gamma amplitude and frequency modulation by theta. The results demonstrate that nonlinearity as a universal/unifying mechanism of all CFC types.

1 Introduction

The Cross Frequency Coupling (CFC) phenomenon may be defined as a statistical correlation between some characteristic parameters of two LFP “rhythms” (Jensen and Colgin, 2007; Jirsa and Müller, 2013; Hyafil et al., 2015). The coupling of theta and gamma rhythms is a well known example, identified in rat hippocampus (e.g., Soltesz and Deschênes, 1993; Bragin et al., 1995; Rosenblum et al., 2000; Colgin et al., 2009; Lisman, 2005; Belluscio et al., 2012; Colgin, 2014; Pernía-Andrade and Jonas, 2014), in the visual cortex of primates (Whittingstall and Logothetis, 2009; Mazzoni et al., 2011), and in humans (Canolty et al., 2006), but CFC phenomena also include sharp waves/ripples (Buzsáki et al., 1992; Skaggs et al., 2007; Bragin et al., 1999; Buzsáki, 2006; Clemens et al., 2011; Buzsáki and da Silva, 2012), and slow oscillations/spindles (Marshall et al., 2006; Rasch and Born, 2013; Staresina et al., 2015; Ladenbauer et al., 2017). CFC has been associated with a long list of cognitive activity aspects, such as local processing and communication between thalamus and neocortex (Fitzgerald et al., 2013); formation of short term memories (Colgin et al., 2009; Lisman and Jensen, 2013; Pernía-Andrade and Jonas, 2014; Lega et al., 2014; Bergmann and Born, 2018); information extraction for decoding motor imagery types (Gwon and Ahn, 2021); as a mechanism for the orienting response (Isler et al., 2008); and many others.

However, despite decades of research, the physical foundation of the CFC process is not well understood (e.g., Jensen and Colgin, 2007; Jirsa and Müller, 2013; Hyafil et al., 2015). The CFC classification into PP, PA, AA, PF, and AF subtypes, basically all two-letter combinations of A, F, and P (amplitude, frequency, and phase) seems taxonomic than rather founded on some mechanistic distinction. Only three (PP, PA, and AA) types have been identified experimentally (Cohen, 2008; Belluscio et al., 2012; De Hemptinne et al., 2013; Munia and Aviyente, 2019; Qin et al., 2021; Chehelcheraghi et al., 2017). Of these, only the PP type, also known as n:m phase synchronization (e.g., Ermentrout, 1981; Jirsa and Müller, 2013), has found a mechanistic description in the limit cycle oscillator model (e.g., Kuramoto, 1975, 1984; Ermentrout, 1991; Strogatz, 2000; Stiefel and Ermentrout, 2016; Nakao, 2016). The origin of other CFC subtypes remains unclear.

This situation is surprising, because all CFC subtypes are in fact aspects of a single overarching physical phenomenon, the coupling of activity across scales—if we think of frequencies as time scales. It is interesting to note that a unifying statistical CFC test does exists: the bispectrum (Jirsa and Müller, 2013; Sheremet et al., 2016, 2018b) subsumes other statistical CFC measures (Kovach et al., 2018). The bispectrum also arises naturally as a statistical characteristic of nonlinear dynamics, and as such it is a universal foundation of CFC physics (Sheremet et al., 2020, 2019). Nonlinearity is of course also the PP coupling mechanism.

Examining the “solved” PP case may shed some light. The limit cycle oscillator model is a highly simplified mathematical description of a population of independent oscillators in a state of sustained periodic oscillations. Oscillators have fixed amplitudes and frequencies, but are allowed weak interaction, which causes mutual phase shifts (Stiefel and Ermentrout, 2016). Because the physical state of the population is completely described by their phases, the model also known as a “phase equation” representation.¹

Because it assumes that amplitudes are constant, a phase model cannot describe AA or AF coupling. To understand amplitude related CFC we need a model that accounts for significant amplitude variations. However, we cannot simply go back to the “parent” model of the phase equation and remove the constant amplitude restriction, because a “parent” model does not exist for neural populations. We need to construct it.

For this we need a physical paradigm of neural activity. We reason that understanding amplitude effects requires a model capable of describing arbitrarily low neural activity, in other words, amplitude emergence. But amplitude emergence means that 1) new frequency emerge; and 2) there is an active energy exchange between frequencies. Transfers of energy across frequencies generate phase correlations, in other words, CFC. These two points are in fact the very definition of nonlinearity.

Therefore, we focus here on neural populations at low activity. In this state, neurons fire intermittently (as opposed to periodically, as limit cycle oscillators) and spend most of the time in the subthreshold regime (e.g., u < −55 mV, where u is the membrane potential; see Figure 1). Action potentials are random and sparse, separated by time intervals at least an order of magnitude larger than the duration of an action potential, and caused by the superposition of presynaptic excitation received from the large number of neighborhood neurons. These states are relevant for populations of recurrently connected Class II neurons (e.g., Hodgkin, 1948; Ermentrout and Terman, 2010; Gerstner et al., 2014) where each neuron has of many random and weak connections (possibly thousands, e.g., Hoppensteadt and Izhikevich, 1997; Muller et al., 2018).

Figure 1

Figure 1. (a) Example of membrane potential evolution of a pyramidal neuron in CA1; 1-s segment intracellular recording data (blue; courtesy of D. English, please refer to English et al., 2014 for details). The red line represents subthreshold fluctuations, with action potentials removed using a simple thresholding algorithm applied to the voltage time derivative). Circles mark the initiation of an action potential (arbitrary voltage reference). (b) PDF function of the threshold for action potential (gray rectangle in (a). In this time series: total time: 3,059 s; total action potentials: 326 seconds; number of action potentials: 25,907; average firing rate of 8.49 Hz.

In these states, the population may transition from a chaotic firing to some degree of oscillatory organization into spontaneous neural assemblies. This organization is dynamically distinct from synchronization of fixed oscillatory circuitry. In this context, the concept of “amplitude” is related to collective neural activity, and is measured by the LFP spectrum. Remarkably, these states exhibit significant CFC patterns, measurable using the general bispectrum test (e.g., Sheremet et al., 2016, 2019).

The goal of this study is to demonstrate the nonlinear foundation of CFC. We use two modeling approaches: direct numerical simulations (DNS) of a population of Hodgkin-Huxley neurons (e.g., Hodgkin and Huxley, 1952), and a simplified mean field approximation based on a two-equation leaky integrate-and-fire (LIF2) model derived from the Hodgkin-Huxley model. The models represent dual, microscopic-macroscopic descriptions. We use the DNS model as a baseline for describing nonlinearity in a neural population. It is straightforward to implement, but its results are difficult to analyze and interpret. We use the mean field approximation to deconstruct the nonlinear mechanism and obtain a simple interpretation of the CFC process. Both models are applied to study the emergence of gamma-like fast oscillations under a slow theta forcing. To identify CFC, we use the spectrum/bicoherence (Jirsa and Müller, 2013; Sheremet et al., 2016, 2018b).

Some experimental evidence is briefly discussed in section (2), where we also introduce and describe briefly the use of the spectrum/bicoherence test for CFC. Section 3 formulates the framework and scope of this study. In Section 4 we discuss direct numerical simulations of a population of fully nonlinear Class II neurons (e.g., Hodgkin and Huxley, 1952; Gerstner et al., 2014) that illustrate the emergence of oscillatory behavior. In Section 5 we derive the simplified model. Slow-fast oscillation coupling in the mean field approximation is analyzed analytically and numerically in Section 6. We summarize and discuss the results in Section 7. The Appendix provides an informal short introduction to bispectral analysis.

2 Some experimental evidence: the spectrum/bicoherence pair

In rat hippocampus, power spectra of local field potential (LFP) recordings exhibit two well studied rhythms: theta, with frequency f≈8 Hz, and gamma, a weaker family of oscillations that cover wide frequency band, approximately 70 Hz. < f < 120. With spatial scales in the order of 10⁻² m and frequency in the order of 100 Hz (speeds of ~ 1 m/s), gamma classifies as mesoscopic activity (Muller et al., 2018).

A useful method to detect rhythms in LFP recordings and estimate their statistical coupling uses the Fourier spectrum and bispectrum (or its modulus, the bicoherence; see Appendix for a brief account of the bispectral estimate). This method has the advantage of requiring minimal data processing, and being comparatively free of assumptions (less prone to a priori fallacies). The spectrum allows for identifying the dominant LFP frequencies, as peaks of power density. The bicoherence is a measure of statistical coupling across the spectrum. If the bicoherence is statistically zero at coordinates (f₁, f₂), the two frequencies are statistically independent; alternatively, a statistically significant (e.g., >0.15, Sheremet et al., 2016) bicoherence peak at coordinates (f₁, f₂) indicates CFC between frequencies f₁, f₂, and f₁+f₂.

Figure 2 illustrates these ideas on using the transformation of LFP with the running speed (Sheremet et al., 2016, 2018a, 2020). As running speed increases, the original theta peak at f≈8 Hz grows and sharpens, and develops prominent secondary peaks at 16 Hz and 24 Hz, clearly visible in the spectrum corresponding to v>15 cm/s (red lines in Figure 2a). The spectrum also shows power accumulating in a broad gamma band spanning frequencies between approximately 70 Hz and 120 Hz (gamma oscillations). In the bicoherence map (Figures 2b, c), the growth of the theta and gamma spectral peaks is accompanied by the development of significant peaks at matching frequencies. For example, the two large peaks at coordinates (8,8) and (16,8) Hz in the bicoherence map (Figure 2c), indicate CFC between theta (8 Hz), and 16 Hz and 24 Hz spectral peaks; this identifies the 16 Hz and 24 Hz oscillations as theta harmonics. Weaker, but statistically significant bicoherence peaks match theta harmonic frequencies (domain marked red triangles in Figures 2b, c) up to the 6th, 48 Hz harmonic. This is remarkable, because there are no detectable corresponding peaks in the spectrum: the bicoherence is able to detect CFC at nascent forced frequencies with very weak amplitudes. The two wide horizontal bands parallel to the x-axis and extending approximately between 75 Hz and 120 Hz indicate CFC between theta and its second harmonic, on the one hand, and gamma on the other.

Figure 2

Figure 2. Example of CFC involving both amplitudes and phases of rat hippocampal LFP rhythms, as a function of rat running speed v (Sheremet et al., 2019, 2020). (a) LFP power density as a function of speed; the spectrum is normalized by f⁻², where f is the frequency. (b, c) Bicoherence maps; peaks indicate CFC at the coordinates of the peak; triangles mark the domains of theta-theta coupling (red) and theta-gamma (blue); see text for discussion.

This brief spectrum/bicoherence analysis shows that emergence of new frequencies (theta harmonics and gamma) and their amplitude growth are intimately connected to CFC. The spectrum/bicoherence test can detect CFC at low activity levels. Below, we use this approach to test for CFC coupling.

3 Formulation of the problem

Here, we consider low activity levels of populations of Class II neurons. To demonstrate the nonlinear mechanism of CFC, we formulate the problem in its simplest form, retaining only the elements strictly necessary to support the emergence of new frequencies: nonlinearity, and the ability to support collective oscillations. The former is an intrinsic property of neural activity; the latter requires some negative feedback, such as inhibitory activity, adaptation, etc (e.g., Palkar et al., 2023, and references therein). For simplicity, we use here only Hodgkin-Huxley neurons (Hodgkin and Huxley, 1952); their intrinsic refractoriness will be used as negative feedback.

3.1 Neuron model

The Hodgkin-Huxley model (Hodgkin and Huxley, 1952; Hodgkin, 1958), describes explicitly the dynamics of ion gates. Because the model is well known, (e.g., FitzHugh, 1955; Abbott and Kepler, 1990; Kepler et al., 1992, and many others), here we only summarize its basic characteristics.

The microscopic model describes the evolution of the unit surface membrane potential u and three gating variables m, n, and h

$\begin{array}{llll}\frac{du}{dt}& =& -j_L-j_{\text{Na}}-j_{\text{K}}+j_{Nt},& (1)\end{array}$

$\begin{array}{l}{j}_L=g_L\left(u-E_L\right);j_{\text{Na}}=g_{\text{Na}}{m}^{3}h\left(u-E_{\text{Na}}\right);\\ j_{\text{K}}=g_{\text{K}}{n}^4\left(u-E_{\text{K}}\right),& (2)\end{array}$

$\begin{array}{llll}\frac{d\xi (u,t)}{dt}& =& -{\kappa }_{\xi }(u)\left[\xi (u,t)-{\xi }_0(u)\right],\xi =m,n,h,& (3)\end{array}$

where j_L, j_Na, j_K, and j_Nt denote ionic currents due to: sodium-potassium pump and passive transport (leaky, L), Na+, and K+ ion channels, and neurotransmitters, respectively. The currents are normalized by the capacity of the cell; E_L, E_Na, E_K are the corresponding reversal potentials; g_L, g_Na, and g_K are constant conductivity values. The evolution of gating variables m, n, and h (Equation 3) is described as relaxation process to equilibrium values ξ₀(u) with characteristic relaxation rate κ_ξ, or relaxation time τ_ξ = 1/κ_ξ. The parameters used in the calculations below are based on Mainen et al. (1995) for a cortex excitatory neuron. The characteristic behavior of the neuron is shown in Figure 3.

Figure 3

Figure 3. Characteristic evolution during an action potential: (a) membrane potential u; (b) gating variables m, h, n; (c) ion currents. The action potential event is marked by a gray rectangle. (d) Voltage-dependent equilibrium values of gating variables m, h and n; (e) Voltage-dependent time constants of the gating variables. vertical lines mark reversal potentials. Circles mark the limits of the subthreshold band E_L ≤ u≾E_T, where E_L = −65 mV is the leaky reversal potential, and E_T≈−55 mV is a characteristic value of the threshold potential.

3.2 Numerical models

The activity of a neural population activity is modeled here using direct numerical simulations (DNS) of neural populations (e.g., Kadmon and Sompolinsky, 2015; Hasegawa, 2000) and a mean field approximation (e.g., Wilson and Cowan, 1972; Jirsa and Haken, 1996, 1997; Coombes et al., 2003; Palkar et al., 2023). The former describes explicitly the dynamics each of the neurons, and can incorporate arbitrarily detailed information of microscopic brain physiology, but are numerically expensive, and produce a wealth of information that is difficult to process and interpret. The latter describe the evolution of population averages; it has direct relevance to macroscopic measurements, it is fast and relatively inexpensive numerically, but involves significant simplifications, and the relationship between its parameter space and physiological quantities is not explicit or simple to obtain.

We use direct numerical simulations of a population of Hodgkin-Huxley neurons under oscillatory forcing to demonstrate the emergence new frequencies, and to test for CFC. We use the mean field approximation to deconstruct the role of nonlinearity in the CFC generation.

3.3 Cautionary note about the interpretation of results

To study the nonlinear foundation of CF we leverage here three means on investigation: experimental data, and two modeling approaches, a DNS and a mean field approximation. These tools provide significantly radically different perspectives on the activity of neural populations.

The numerical models represent a highly simplified neural system, a neural mass of Class II point neurons with refractoriness, under a simple oscillatory forcing. Experimental data is infinitely more complicated that that.

The numerical models themselves are not just different, they represent two opposite perspectives: microscopic in the DNS case, macroscopic in the mean field case. In principle, they could be calibrated to produce matching results, but this is an effort beyond the scope of this study. Getting the models to agree with the experimental data is obviously impossible in the simplified setup used her.

However, our purpose is to demonstrate and investigate the nonlinear mechanism of CFC in populations of neurons at low activity levels. For this, a qualitative agreement between patterns in the CFC spectral/bicoherence test suffices.

4 The DNS model

The DNS model uses the Brian2 numerical simulator for spiking neural networks (Stimberg et al., 2019) to integrate equations (Equations 1–3) for a population of 1,000 interconnected neurons. Neurons receive input through 20 connections on average, and have on average 30 recurrent connections. The external forcing is modeled as Poisson time series corresponding to a firing rate of the form F(t) = F₀+Acos2πft, where F₀ = 100 Hz, A = 10 Hz, and f = 8 Hz. The 8 Hz modulation emulates the theta oscillation.

The simulations confirm that a population of Hodgkin-Huxley neurons exhibits emergent oscillations under periodic forcing. Figure 4 summarizes some essential aspects of the population activity that illustrate our concept of low activity. The population responds to the modulated input with an oscillation on the same frequency, clearly detectable in the mean firing rate, membrane potential u, and gating variable h (Figures 4a, c, d). However, individual neurons fire intermittently (Figure 4b), with action potentials separated by large time intervals spent in subthreshold regime. The oscillations represent a spontaneous organization of collective activity oscillation, driven by input, rather than synchronization of limit cycle activity of individual neurons. (Synchronization may be achieved if the input is strong enough to drive all neurons into limit cycles).

Figure 4

Figure 4. Response of a population of Hodgkin-Huxley neurons to slow oscillatory forcing. (a) Mean input action potential rate received by a population neuron (blue) and mean firing rate of a population neuron. (b) Time series of the membrane potential of neuron; the subthreshold regime is marked by a gray band. (c, d) Population mean membrane potential u and gating variable h. (e) Bicoherence of the mean membrane potential; (f) normalized spectral density of mean membrane potential and input firing rate.

The results also show the hallmark signs of nonlinear CFC. While the mean input is Gaussian (not skewed, and symmetric), the response of the population shows significant skewness (positive for mean u, negative for mean h), and asymmetry (negative for the mean firing rate and mean u; positive for the mean h; Figures 4a, c, d). This deformation is associated with the appearance of higher harmonics of the theta frequency. The spectral density shows four, maybe five peaks at theta harmonics, with matching peaks are prominent in the bicoherence map (Figures 4e, f).

The spectrum exhibits also a broad gamma peak approximately between 80 and 180 Hz, with the bicoherence peaks coupling gamma and theta (both 8 Hz peak and 16 Hz harmonic).

This behavior is qualitatively similar to the measurements discussed in Section 2. Compared to experimental data, the most obvious differences are the strength of the CFC effect (number and height of peaks in both spectrum and bicoherence), and the higher frequency of the “gamma” band, possibly due to the reliance of the DNS model on refractoriness as inhibitory feedback.

The results demonstrate that slow oscillatory forcing of a population of interconnected Hodgkin-Huxley neurons at low activity levels causes the emergence of new, high-frequency oscillations, that show statistically significant CFC patterns.

5 The mean field approximation model

To develop a mean field approximation compatible with the Hodgkin-Huxley population model used in the previous section, we first derive a lower dimension approximation of equations (Equations 1–3) and then apply a population average. The procedure follows the classical work of FitzHugh (1955), Abbott and Kepler (1990), and Kepler et al. (1992). The main distinction between their work and the present approach is the treatment of action potentials. Because we expect neurons to spend most of the time in the subthreshold regime, with intermittent and scarce action potentials (Figures 1–4b), we assume that the duration of action potentials is much smaller than subthreshold residence time scales and neglect action potentials.

5.1 Dimension reduction for subthreshold dynamics

The dimension reduction of equations (Equations 1–3) is based on the following observations. In the subthreshold regime E_L ≤ u≾−55 mV (see Figure 3):

(1) The small relaxation time of m (Figure 3e) indicates that it relaxes quickly to its equilibrium value; we neglect its explicit time dependency and replace it everywhere with m₀(u).

(2) The dependence of h on u is weak (Figure 3d); we a write h(u, t)≈h(u₀, t), with h₀(u₀) = h₀₀ and κ_h = κ_h(u₀) = κ_h0 constants u₀ is a reference value, say u₀ = E_L.

(3) In the subthreshold regime n≪1, and it appears at fourth power in the potassium current; we neglect n dynamics altogether.

Therefore, the subthreshold dynamics of equations (Equations 1–3) may be approximated as

$\begin{array}{llll}\frac{d}{dt}u(t)& =& j_L(u)+j_{\text{Na}}(u,t)+j_N(t),& (4)\end{array}$

$\begin{array}{l}{j}_L=g_L\left(u-E_L\right),j_{\text{Na}}=\frac{h(t)}{h_{00}}G(u),\text{ with}\\ G(u)=g_{\text{Na}}{h}_{00}{m}_0^3(u)\left[E_{\text{Na}}-u(t)\right],& (5)\end{array}$

$\begin{array}{llll}\frac{d}{dt}\frac{h(t)}{h_{00}}& =& -{\kappa }_{h0}\left[\frac{h(t)}{h_{00}}-1\right].& (6)\end{array}$

This approximation is not complete, however, because is does not account for the resetting effect an action potential.

To correct this, we note that the wide distribution of the threshold membrane potential (Figure 1b; see also Gerstner et al., 2014) indicates that a fixed-value threshold criterion is a severe simplification (e.g., Koch et al., 1995; Rinzel and Ermentrout, 1998). Adopting the more natural threshold criterion j_Na≈j_L codifies the variability of the threshold u value in Equation 4: the small post-spike value of h (Figures 3b, c) reduces the sensitivity of the sodium channel to membrane potential, and consequently, the ability of the neuron to fire.

To switch to a threshold criterion based on j_Na, it is convenient to change state variable from u to j_Na. Differentiating j_Na in Equation 5 to time, and using Equations 4, 6, after some straightforward algebra obtains

$\begin{array}{ll}\frac{d}{dt}{j}_{\text{Na}}=\frac{h}{h_{00}}{\phi }_N(j_{\text{Na}},j_N)-{\kappa }_j(j_{\text{Na}},h)j_{\text{Na}},& (7)\end{array}$

where φ_N represents rate of change rate j_Na induced by neurotransmitters in absence of refractoriness

$\begin{array}{l}{\phi }_N\left(j_{\text{Na}},j_N\right)=h{G}^{\text{'}}{j}_N,{\kappa }_j\left(j_{\text{Na}},{\varphi }_{Nt}\right)=-\left(h{G}^{\text{'}}\frac{j_L}{j_{\text{Na}}}+h{G}^{\text{'}}\right)\\ -\frac{1}{{\tau }_{h0}}\left(1-\frac{h}{h_{00}}\right),& (8)\end{array}$

and all dependencies on u become dependencies on j_Na and t. Replacing Equation 4 with 7 completes the change of variable from u to j_Na.

We can now approximate the threshold criterion as j_Na = j_T, where j_T is characteristic (constant) value, and represent the action potential resetting as a Dirac delta function (e.g., Roxin et al., 2011):

$\begin{array}{lll}\frac{d}{dt}\frac{j_{\text{Na}}}{j_T}& =\frac{h}{h_{00}}{\phi }_N-{\kappa }_j\left(\frac{j_{\text{Na}}}{j_T},\frac{h}{h_{00}}\right)\frac{j_{\text{Na}}}{j_T}-\frac{j_{\text{Na}}}{j_T}\delta \left(\frac{j_{\text{Na}}}{j_T}-1\right);& (9)\end{array}$

$\begin{array}{lll}\frac{d}{dt}\frac{h}{h_{00}}& =-{\kappa }_{h0}\left(\frac{h}{h_{00}}-1\right)-\frac{h}{h_{00}}\delta \left(\frac{j_{\text{Na}}}{j_T}-1\right),& (10)\end{array}$

where we assumed for simplicity that, post spike, j_Na = 0 and h = 0.

Equations 9, 10 represent a two-equation “leaky integrate and fire” model (LIF2). The model follows the same principle as other LIF models (e.g., Lapicque, 1907; Hill, 1936; Abbott, 1999; Dong et al., 2011), with the addition of Equation 10 accounting for excitability (h).

5.2 Mean field approximation

For a neural mass with no spatial dimensions, averaging equations (Equations 4–6) over the population and assuming that h and j_Na are uncorrelated obtains

$\begin{array}{llll}\frac{d}{dt}J& =& H\varphi -{\kappa }_{j}J-\langle \frac{j_{\text{Na}}}{j_T}\delta \left(\frac{j_{\text{Na}}}{j_T}-1\right)\rangle ,& (11)\end{array}$

$\begin{array}{llll}\frac{d}{dt}H& =& -{\kappa }_{h0}\left(H-1\right)+\langle \frac{h}{h_{00}}\delta \left(\frac{j_{\text{Na}}}{j_T}-1\right)\rangle ,& (12)\end{array}$

where 〈·〉 is the average operator, J = 〈j_Na/j_T〉, H = 〈h/h₀₀〉, and ϕ = 〈φ_N〉, with and κ_j≈κ_j0 constant. The δ terms are nonzero only for neurons firing in the interval [t, t+dt], i.e., they count the action potentials in the population, recorded in the unit of time, divided by the number of neurons in the population. This quantity is the firing rate N,

$\begin{array}{ll}N=\frac{1}{n}{\displaystyle \sum _{j=1}^{\nu }}\frac{n_s^j{\tau }_s}{dt},& (13)\end{array}$

where ν is the number of neurons firing in [t, t+dt], $n_s^j{\tau }_{sj}$ approximates the time spent firing by neuron j, τ_s is the mean duration of a spike, and $n_s^j$ is the number of spikes. Further writing ϕ = μ(N+Q), where Q is the external input, described as a firing rate, and μ is a “connectivity” parameter that describes the conversion of the presynaptic firing rates N+Q into postsynaptic j_N currents, Equation 11 become

$\begin{array}{ll}\frac{d}{dt}J=H\mu (N+Q)-{\kappa }_{j0}J-N,& (14)\end{array}$

$\begin{array}{ll}\frac{d}{dt}H=-{\kappa }_{h0}(H-1)-N{H}_T,& (15)\end{array}$

where H_T = h_T/h₀₀ is the threshold value of H. The system (Equations 14, 15) is closed by prescribing the relation between the state variable J and the firing rate N, also known as the activation function. Because τ_s is assumed infinitesimal, the time interval Δt → 0 can contain an arbitrary number of spikes. In general, Equation 13 implies that N≥0 and can be arbitrarily large at maximum activity level J = 1 when all neurons fire continuously. Here, we use the activation function derived by Qin and Sheremet (2022)

$\begin{array}{ll}N(J)=\alpha \frac{J}{1-J},& (16)\end{array}$

where α is the “susceptibility” parameter, related of the intensity of endogenous membrane potential fluctuations, which may also be understood as the susceptibility of a neuron to fire. The activation function (Equation 16) is approximately linear for J≪1 and is unbounded as J↗1.

Equations 14, 16 are the mean field approximation of the LIF2 model. The state of the neural field is described by two state variables: the mean current density J and the excitability H. Both are normalized 0 ≤ J, H ≤ 1; where H = 0 after an action potential and H = 1 in the resting state; and J = 0 after an action potential and J = 1 at the action potential threshold.

The equations have a straightforward interpretation in terms of energy balance. The variable J may be interpreted as the kinetic energy density (also, intuitively, as the “temperature”, by analogy with an ideal gas, where temperature and kinetic energy density are proportional). The firing rates N (internal) and Q (external, Equation 14) represent the energy flux into the system due to internal and external action potentials; the conversion into local kinetic energy is controlled by the excitability H and the synaptic connectivity parameter μ. The second term in Equation 14 is the natural dissipation rate of the kinetic energy density due leaky currents and refractoriness. The third describes energy loss due to post action potential refractoriness.

The neural mass is collectively characterized by four “material” parameters. Two characterize the mean neuron – the relaxation rates κ_h0 (Equation 15) and κ_j0 (Equation 15). The other two characterize the inter-neuron communication—the connectivity μ (Equation 14), and susceptibility α (Equation 16).

The mean field approximation serves the purpose of this demonstration by retaining the two elements we consider essential for CFC generation: negative feedback, represented by the excitability H; and nonlinearity, represented by the activation function and the first term in Equation 14.

6 Nonlinear coupling between slow and fast oscillations

Equations 14–16 represent an approximation of the DNS model described in Section 4, and can be used to gain insight into the CFC mechanism. For this, we use a standard asymptotic approach (e.g., Bender and Orszag, 1991) that takes advantage of the different scales of the coupled oscillations. For simplicity, all material parameters are approximated as constants.

In the theta-gamma coupling, theta is a low frequency, high power oscillation, while gamma a high frequency, low power one. This relationship may be formalized by introducing a small parameter ϵ≪1 and writing

$\begin{array}{llll}J(t,ϵt)& =& J_0(ϵt)+ϵJ_1(t,ϵt)+O({ϵ}^2),& (17)\end{array}$

$\begin{array}{llll}H& =& H_0(ϵt)+ϵH_1(t,ϵt)+O({ϵ}^2),& (18)\end{array}$

$\begin{array}{l}N=N_0(ϵt)+ϵN_1(t,ϵt),\text{ where }{N}_0=N(J_0),\text{ and}\\ N_1={\frac{\partial N}{\partial J}|}_0{J}_1,& (19)\end{array}$

where $\frac{\partial N}{\partial J}$ is a variational derivative. From Equation 16 one obtains

$\begin{array}{ll}\frac{\partial N}{\partial J}=\alpha {(1-J)}^{-2}\approx \alpha (1+2J+\dots ).& (20)\end{array}$

Theta (the slow oscillation, subscript 0), is represented as an order 1 component that evolves on the “slow” time ϵt. Gamma, (the fast oscillation, subscript 1), is an order ϵ component and assumed to evolve on two time scales: a “fast” time t, with a slow modulation of ϵt scale induced by theta. Substituting into Equations 14–16 and separating the powers of ϵ obtains, after some straightforward algebra, two sets of equations.

6.1 Leading order system

Over the fast t time scales the slow variables satisfy the equilibrium (no explicit time dependency) equations

$\begin{array}{ll}{H}_0(\mu N_0+Q)={\kappa }_{j0}{J}_0+N_0,& (21)\end{array}$

$\begin{array}{ll}{\kappa }_{h0}(1-H_0)-N_0{H}^T=0& (22)\end{array}$

Equation 21 represents the balance between incoming energy flux and energy losses due to the natural relaxation rate κ₀ and post spike resetting. Equation 22 describes the rate of growth of population refractoriness (1−H) as a consequence of action potentials.

From Equation 16, $J_0=\frac{N_0}{\alpha +N_0}$; solving Equation 22 in H₀ and substituting into Equation 21 obtains for the equilibrium firing rate the cubic equation

$\begin{array}{ll}(1-{\tau }_{h0}{H}^T{N}_0)(\mu N_0+Q)(\alpha +N_0)-{\kappa }_{j0}(\alpha +N_0)N_0-N_0=0,& (23)\end{array}$

where τ_h0 = 1/κ_h0. Keeping fixed the mean neuron material parameters κ_h0 and κ_j0 implies that the solutions of Equation 23 for a given external forcing Q₀ are uniquely defined by the network material parameters μ and α.

6.2 Second order system

At order ϵ, setting $\frac{\partial N}{\partial J}{|}_0{J}_1\approx \alpha J_1$ for low activity levels obtains

$\begin{array}{ll}\left[I\frac{d}{dt}-\Omega (J_0,H_0)\right]\left[\begin{array}{c}{J}_1\\ H_1\end{array}\right]=-\frac{d}{dt}\left[\begin{array}{c}{J}_0\\ H_0\end{array}\right];& (24)\end{array}$

where I is the the identity matrix, and Ω is a 2 × 2 matrix with entries

$\begin{array}{l}{\Omega }_{11}=H_0\mu \alpha -{\kappa }_{j0}-\alpha ,{\Omega }_{12}=\mu \alpha J_0+Q,{\Omega }_{21}=-\alpha H^T,\text{ and}\\ {\Omega }_{22}=-{\kappa }_{h0}.& (25)\end{array}$

Equation 24 describes the stability of the equilibrium states with respect to small perturbations (J₁, H₁). Seeking solutions for J₁and H₁∝e^λt obtains the eigenvalues

$\begin{array}{l}{\lambda }_{1,2}=\kappa \pm i\omega ,\text{ with }\kappa =\frac{1}{2}T,\omega =\frac{1}{2}\sqrt{\Delta },\text{ and}\\ \Delta =T^2-4D,& (26)\end{array}$

where T and D are the trace and the determinant of Ω,

$\begin{array}{lll}T& ={\Omega }_{11}+{\Omega }_{22}=\mu \alpha H_0-{\kappa }_{j0}-\alpha -{\kappa }_{h0},& (27)\end{array}$

$\begin{array}{l}D={\Omega }_{11}{\Omega }_{22}-{\Omega }_{12}{\Omega }_{21}=-{\kappa }_{h0}\left(\mu \alpha H_0-{\kappa }_{j0}-\alpha \right)\\ +\alpha H^T\left(\mu \alpha J_0+Q_0\right).& (28)\end{array}$

The stability of equilibrium states is also completely defined by the network parameters μ and α. Equilibrium states are unstable to small perturbations if growth rate κ>0 and stable if κ ≤ 0; if Δ < 0, the perturbations are oscillatory with eigenfrequency ω. Figures 5b, c shows the stability characteristics of the equilibrium states.

Figure 5

Figure 5. Distribution of characteristic parameters of the linearized equations (Equations 23, 24) in the (μ, α) parameter space for constant forcing Q₀. Each (μ, α) point in plane represents a unique networks configuration. (a) Current density J. (b) Excitability H. (c) Firing rate N. (d) Frequency ω. (e) Growth rate κ. (f) Schematic map of the stability of equilibrium states. The point A marks the weakly stable configuration used in the numerical simulations.

6.3 Linear equilibrium and stability

Figure 5 shows the characteristics of the linearized mean field Equations 23, 24 in the (μ, α) parameter space.

Equilibrium states correspond to the real roots of Equation 23. Triple equilibrium states exist for low susceptibility α and strong connectivity μ (typically with two stable points separated by an unstable one) where they occupy the cusp-like white region (this is a typical feature for this kind of systems, e.g., Ermentrout, 1998). Low susceptibility maintains stability at low activity by reducing the effect of endogenous activity on firing rates; strong connectivity ϵ maintains stability at higher activity levels by sustaining a high firing rate and recapturing a large amount of the energy released.

Low connectivity μ and low susceptibility α (lower left corner) correspond to an inefficient conversion of N₀ energy fluxes, and the system is driven primarily by external forcing Q. The opposite is true for high connectivity and susceptibility (upper right corner). As J₀ and N₀ increases, higher order terms in the expansion (Equations 17–19) become important, the relationship between (J₀, N₀) and external input Q₀ weakens, and the stability of the equilibrium point decreases.

6.4 The nonlinear CFC mechanism for slow/fast oscillations

Because the system of equations (Equation 24) is linear, it may be solved using elementary methods. Without going into details, if ω≠0, the eigenvalues (Equation 26) are distinct and the system may be uncoupled by diagonalizing the Ω matrix

$\begin{array}{ll}(I\frac{d}{dt}-\Lambda )P_1=-\frac{d}{dt}{P}_0,P_{0,1}=S\left[\begin{array}{c}{J}_{0,1}\\ H_{0,1}\end{array}\right]& (29)\end{array}$

where Λ is a diagonal matrix with entries λ_{1, 2} and S is the matrix of eigenvectors. The solutions of the problem (Equation 29) have the form

$\begin{array}{l}{P}_1^j=A^j(t)e^{-i{\theta }^j},\text{ with }{\theta }^j=-i{{\displaystyle \int }}^t\lambda (\zeta )d\zeta ,\text{ and}\\ A^j(t)=-{{\displaystyle \int }}^t\frac{d{P}_0^j(\zeta )}{d\zeta }{e}^{i{{\displaystyle \int }}^{\zeta }\lambda (\xi )d\xi }d\zeta .& (30)\end{array}$

where $P_{0,1}^j$ are the j-th, (j = 1, 2) components of the vectors P_{0, 1}, and λ is one of the two eigenvalues.

Note that in Equation 30 the amplitude and A^j, the phase θ^j, and frequency λ of the fast (gamma) oscillation P*j₁ all depend on the variables J₀ and H₀ that characterize the slow (theta) oscillations via the eigenvalues λ_{1, 2}. This solution describes explicitly the role of nonlinearity in coupling between all theta and gamma characteristics parameters, including all combination of amplitudes, phases, and frequencies.

The linearization achieved by expansion (Equations 17–19) preserves the nonlinear terms, but breaks them into factors of different scales and assigning their contribution accordingly. For example, the nonlinear term HJ in Equations 14–16 appears in Equation 24 as H₀J₁, mixing slow (theta, subscript 0) and fast variables. This reduces the evolution of the slow component (J₀, H₀) as sequence of equilibrium states (Equations 21, 22), and describes the fast oscillations (J₁, H₁) as a transient process of adjustment to the next equilibrium state (Equation 24).

Because throughout this evolution the system is approximately at equilibrium, each state may be represented in the (μ_A, α_A), Figure 5. A slowly varying forcing Q₀ shifts the topology of equilibrium regions with respect to a fixed point A = (μ_A, α_A). As Q₀ increases, the triple equilibrium cusp flattens, and shifts to a lower position relative to A (Figure 6a). The values of current density J, the eigenfrequency ω and the growth rate κ at point A also increase. Stability decreases globally in the (μ, α) space.

Figure 6

Figure 6. (a) Shifting of equilibrium stability as a function of forcing Q₀ (Equation 23). Compare with Figure 5. As the forcing increases, the configuration A becomes unstable. In the linearized model (Equation 24) the high frequency oscillations are unbounded. (b) In fact, the full nonlinear system undergoes a Hopf bifurcation from node to stable orbits. Theta-gamma CFC may be visualized as a slow Q₀ running through this range of values, periodically destabilizing and re-stabilizing the high frequency oscillations (gamma, limit cycles).

In summary, the linearized equations (Equations 23, 24) describe the nonlinear mechanism of CFC between slow and fast oscillations (theta-gamma) as a cyclic destabilization of the eigenfrequencies of the system.

6.5 Numerical simulations

The validity of the linearized equations (Equations 23, 24) is limited by the assumptions used by expansion (Equations 17–19). Whether it is realistic of not to assumes that theta and gamma are separated in both power and frequency by, say, ϵ~10⁻², such a scaling only describes short time behavior, because unstable eigenoscillations can grow over time and break the original ordering. To obtain long term solutions, Equations 14–16 must be integrated numerically.

6.5.1 Numerical setup

For numerical integration of the full system of equations (Equations 14–16), we represent forcing Q as a non-skewed, symmetric stochastic process with narrow-spectrum Gaussian pink noise ∝f^−1.5 superposed on a monochromatic oscillation of 8 Hz (Figure 7). The noise is added in order to provide a frequency smoothing for spectral estimators, but has otherwise no dynamical effect. By construction, the forcing time series does not contain any cross-frequency coupling, therefore its spectrum only exhibits the 8 Hz peak and the bicoherence is statistically zero (Figures 7b, c).

Figure 7

Figure 7. A segment of the time series Q(t) used as forcing in both linear and nonlinear numerical simulations. (a) A 1-second segment of the forcing time-series; (b) Spectral density of Q(t) power; (c) Q(t) bicoherence.

Here, we present numerical solutions for point A in Figures 5, 6, chosen because it is weakly stable at approximately middle of the Q oscillation range (Figure 6a; see also Figure 5e). In these conditions, eigenoscillations are stable throughout most of the lower forcing values, but unstable at high Q values.

To enable some comparison with measurements, we construct here a mock LFP by the simple superposition

$\begin{array}{ll}\text{LFP}~(Q+J).& (31)\end{array}$

In the CA1 region in the hippocampus, for example, Q might represent external synaptic excitation excitation though Schaffer collateral path, while J might represent the current source near stratum pyramidale. Both these signals contribute to LFP recordings in the rat hippocampus. However, this description does not account for other synaptic and somatic activities such as spiking and hyperpolarization activities, and other ionic processes (Buzsáki et al., 2012).

6.5.2 Results

Figures 8, 9 show that numerical solutions of the linearized (Equation 24) and the fully nonlinear mean field approximation (Equation 14–16) are qualitatively similar behavior. The gamma stabilization/destabilization cycle is obvious in the time series of the current density (Figures 8a, 9a). The effect is dampened in the mock LFP by the contribution of the Q (Figure 7a).

Figure 8

Figure 8. Numerical solution of the linear approximation (Equation 24) for neural network A (Figure 6), under the forcing Q defined in Figure 7. (a, b) time current density and mock LFP. (c) Spectral densities of current density J, forcing Q and Q+J, the mock LFP. (d, e) Bicoherence of the current density J, and Q+J, mock LFP.

Figure 9

Figure 9. Numerical solution of the fully nonlinear mean field approximation (Equation 14–16) for neural network A (Figure 6), under the forcing Q defined in Figure 7. (a, b) time current density and mock LFP. (c) Spectral densities of current density J, forcing Q and Q+J, the mock LFP. (d, e) Bicoherence of the current density J, and Q+J, mock LFP.

However, In the linearized representation, fast gamma oscillations are described as noise preferentially amplified at unstable eigenfrequencies, and can grow to arbitrarily large values. In the fully nonlinear model, the growth mechanism is initially the same, but as the gamma amplitudes increase, nonlinear terms become important and limit the growth. The solutions of the nonlinear system remain globally stable, and exhibit a bifurcation of the supercritical Hopf type (Figure 6b). This effect limits the power in the gamma band, as seen by comparing spectral densities (Figures 8c, 9c).

Theta harmonics are visible in both the spectrum and bicoherence. The bicoherence also exhibits peaks indicating significant phase correlations between theta and gamma (Figures 8c–e, 9c–e). The patterns of CFC between theta and gamma are in qualitative agreement with DNS results and experimental data. The also confirm that CFC between theta and gamma can arise strictly as a nonlinear response of the neural network to a Gaussian forcing.

7 Summary and discussion

The goal of this study is to demonstrate and analyze the nonlinear foundation of the CFC phenomenon.

The modeling framework is based here on a paradigm of “irregular”, low-intensity neural activity, defined as a state where neurons fire intermittently and spend most of the time in the subthreshold regime. These conditions allow neural population to organize into oscillatory patterns, with amplitudes measured by the LFP spectrum. These states are relevant for populations of Class II neurons, and have been shown to exhibit significant CFC bispectral patterns (e.g., Sheremet et al., 2016, 2019).

The nonlinear low-activity paradigm is supported by the analysis of LFP rat hippocampus recordings. The LFP spectrum and bicoherence (Jirsa and Müller, 2013; Sheremet et al., 2016, 2018b) shows distinctive CFC patterns that indicate coupling between theta and its harmonics, and between theta and gamma. Importantly, the growth of theta-gamma coupling with the power (amplitude) of the oscillations indicate that amplitude and phase types of CFC are related.

The numerical modeling framework discussed here has two components: a DNS model of interacting (Hodgkin and Huxley, 1952) neurons, and a simplified mean field approximation based on two-equation LIF model. The mean field equations are nonlinear, and retain the negative firing rate feedback, but have a much simpler structure than the DNS population model. In both models, the neural population is forced by a slow (theta) oscillation. In both models, Spectral/bicoherence of the mean membrane potential time series exhibit patterns similar to LFP recordings.

The DNS model is used to verify the nonlinear CFC mechanism directly. However, because the complexity of the DNS model makes further analysis difficult, to understand nonlinear CFC mechanism we turn to the mean field approximation model. An asymptotic expansion that takes advantage of the scale differences between theta and gamma can readily be used to linearize the system. This provides a few notable results:

(1) The linearized model describes the coupling between theta and gamma as a eigenfrequency (gamma) destabilization/restabilization cycle forced by theta: underlying mechanism of gamma forcing by theta is nonlinearity.

(2) The analytical solution of the linearized equations gives an explicit expression of the effect of nonlinearity on all gamma characteristics: amplitude, phase, and frequency. This demonstrates the role of nonlinearity as a unifying mechanism for all CFC types identified in literature. Indeed, the simplified modeling framework we use here allows for writing explicitly (Equations 26, 30) the relationship between all characteristics of the slow and fast (theta-gamma) oscillations: amplitude, phase, and frequency.

(3) The “gamma” frequencies produced by our simulations are different (higher) than experimentally observed. This is likely the result of relying on the short-time intrinsic Hodgkin-Huxley refractorinessU for negative feedback. Indeed, from Equations 24–30, the frequency λ (Equation 26) of the fast oscillation depends on the characteristics decay rates κ (Equations 27–28), originally defined in Equation 6 for the refractoriness, and Equation 8 for the sodium current. Note that in this mechanism, gamma oscillations are indeed controlled by the time scale of the negative feedback (e.g., in Equation 15, τ_h0 → 0 implies H → 1 and nonlinearity is disabled). However, further research into this dependency is beyond the scope of this study, and would be of doubtful usefulness, given the simplicity of the models used here (inhibition largely overrides these effects).

These results show that, under theta forcing, a population of Hodgkin and Huxley (1952) neurons can spontaneously generate high frequency (gamma) oscillations due to the intrinsic nonlinearity of the neural activity. The theta-gamma system of oscillations has CFC characteristics similar to observations. This suggests the nonlinearity of brain activity is a fundamental mechanism for the generation of theta-gamma CFC (see also Sheremet et al. 2019, 2020). This is consistent with general physics results and theories for the evolution of large systems as diverse as water waves, plasma, aggregation-fragmentation in particulate systems, chemical reactions, etc (e.g., L‘vov 1998; Zakharov et al. 1992; Connaughton et al. 2006; Proment et al. 2009). Our results suggest that nonlinearity also plays a role in brain activity CFC processes and cognitive activities.

However, we should stress that albeit fundamental, nonlinearity is only one of a number of CFC mechanisms identified in brain activity at various scales, that include, e.g., the celebrated ING/PING mechanisms based on inhibition, as well as interaction between different brain regions (e.g., Bartos et al., 2007; Tiesinga and Sejnowski, 2009; Buzsáki and Wang, 2012, and many others). Further research is needed to understand the importance of the nonlinear element in this complex context of brain activity.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

AS: Conceptualization, Formal analysis, Investigation, Methodology, Supervision, Visualization, Writing – original draft, Writing – review & editing. YQ: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Visualization, Writing – original draft.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. This work was supported by the McKnight BrainResearch Foundation, National Institute on Aging, grant number AG055544 and National Institute of Mental Health, grant number MH109548.

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.

Generative AI statement

The author(s) declare that no Gen AI was used in the creation of this manuscript.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

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

Footnotes

1. ^To use a phase model one has to find oscillators. Initially, the role of oscillators was played by Hodgkin-Huxley neurons (Hodgkin and Huxley, 1952) and derivative models (e.g., FitzHugh, 1955; Abbott and Kepler, 1990; Ermentrout and Kopell, 1986, and others). In time, neural “circuitry” has also been considered, such as central pattern generators (e.g., Marder and Bucher, 2001, probably not applicable to the cortex); interacting brain regions (e.g., Chehelcheraghi et al., 2017); and others.

References

Abbott, L. F. (1999). Lapicque's introduction of the integrate-and-fire model neuron (1907). Brain Res. Bull. 50, 303–304. doi: 10.1016/S0361-9230(99)00161-6

Crossref Full Text | Google Scholar

Abbott, L. F., and Kepler, T. B. (1990). “Model neurons: From Hodgkin-Huxley to Hopfield,” in Statistical Mechanics of Neural Networks, L. Garrido (Berlin: Springer Berlin Heidelberg), 5–18.

Google Scholar

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

PubMed Abstract | Crossref Full Text | Google Scholar

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

PubMed Abstract | Crossref Full Text | Google Scholar

Bender, C. M., and Orszag, S. A. (1991). Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory. Cham: Springer.

Google Scholar

Bergmann, T. O., and Born, J. (2018). Phase-amplitude coupling: a general mechanism for memory processing and synaptic plasticity? Neuron 97, 10–13. doi: 10.1016/j.neuron.2017.12.023

PubMed Abstract | Crossref Full Text | Google Scholar

Bragin, A., Engel, J., Wilson, C. L., Fried, I., and Buzsáki, G. (1999). High-frequency oscillations in human brain. Hippocampus 9, 137–142.3.0.

Google Scholar

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

PubMed Abstract | Crossref Full Text | Google Scholar

Buzsáki, G. (2006). Rhythms of the Brain. Oxford: Oxford University Press.

Google Scholar

Buzsáki, G., Anastassiou, C. A., and Koch, C. (2012). The origin of extracellular fields and currents-eeg, ecog, lfp and spikes. Nat. Rev. Neurosci. 13, 407–420. doi: 10.1038/nrn3241

PubMed Abstract | Crossref Full Text | Google Scholar

Buzsáki, G., and da Silva, F. L. (2012). High frequency oscillations in the intact brain. Prog. Neurobiol. 98, 241–249. doi: 10.1016/j.pneurobio.2012.02.004

PubMed Abstract | Crossref Full Text | Google Scholar

Buzsáki, G., Horváth, Z., Urioste, R., Hetke, J., and Wise, K. (1992). High-frequency network oscillation in the hippocampus. Science 256, 1025–1027. doi: 10.1126/science.1589772

PubMed Abstract | Crossref Full Text | Google Scholar

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

PubMed Abstract | Crossref Full Text | Google Scholar

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

PubMed Abstract | Crossref Full Text | Google Scholar

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

PubMed Abstract | Crossref Full Text | Google Scholar

Clemens, Z., Mölle, M., Erőss, L., Jakus, R., Rásonyi, G., Halász, P., et al. (2011). Fine-tuned coupling between human parahippocampal ripples and sleep spindles. Eur. J. Neurosci. 33, 511–520. doi: 10.1111/j.1460-9568.2010.07505.x

PubMed Abstract | Crossref Full Text | Google Scholar

Cohen, M. X. (2008). Assessing transient cross-frequency coupling in EEG data. J. Neurosci. Methods 168, 494–499. doi: 10.1016/j.jneumeth.2007.10.012

PubMed Abstract | Crossref Full Text | Google Scholar

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

PubMed Abstract | Crossref Full Text | Google Scholar

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

PubMed Abstract | Crossref Full Text | Google Scholar

Connaughton, C., Rajesh, R., and Zaboronski, O. (2006). Cluster-cluster aggregation as an analogue of a turbulent cascade: Kolmogorov phenomenology, scaling laws and the breakdown of self-similarity. Physica D 222, 97–115. doi: 10.1016/j.physd.2006.08.005

Crossref Full Text | Google Scholar

Coombes, S. L G. J, and Owen, M. R. (2003). Waves and bumps in neuronal networks with axo-dendritic synaptic interactions. Physica D 178, 219–241. doi: 10.1016/S0167-2789(03)00002-2

Crossref Full Text | Google Scholar

De Hemptinne, C., Ryapolova-Webb, E. S., Air, E. L., Garcia, P. A., Miller, K. J., Ojemann, J. G., et al. (2013). Exaggerated phase-amplitude coupling in the primary motor cortex in Parkinson disease. Proc. Nat. Acad. Sci. 110, 4780–4785. doi: 10.1073/pnas.1214546110

PubMed Abstract | Crossref Full Text | Google Scholar

Dong, Y., Mihalas, S., Russell, A., Etienne-Cummings, R., and Niebur, E. (2011). Estimating parameters of generalized integrate-and-fire neurons from the maximum likelihood of spike trains. Neural Comput. 23, 2833–2867. doi: 10.1162/NECO_a_00196

PubMed Abstract | Crossref Full Text | Google Scholar

English, D. F., Peyrache, A., Stark, E., Roux, L., Vallentin, D., Long, M. A., et al. (2014). Excitation and inhibition compete to control spiking during hippocampal ripples: intracellular study in behaving mice. J. Neurosci. 34, 16509–16517. doi: 10.1523/JNEUROSCI.2600-14.2014

PubMed Abstract | Crossref Full Text | Google Scholar

Ermentrout, G. B. (1981). n:m Phase-Locking of Weakly Coupled Oscillators. J. Math. Biol. 12, 327–342. doi: 10.1007/BF00276920

Crossref Full Text | Google Scholar

Ermentrout, G. B. (1998). Neural networks as spatio-temporal pattern-forming systems. Reports Prog. Phys. 61:353. doi: 10.1088/0034-4885/61/4/002

Crossref Full Text | Google Scholar

Ermentrout, G. B., and &Kopell, N. J. (1991). Multiple pulse interactions and averaging in systems of coupled neural oscillators. Math. Biol. 29, 195–217. doi: 10.1007/BF00160535

Crossref Full Text | Google Scholar

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

Crossref Full Text | Google Scholar

Ermentrout, G. B., and Terman, D. H. (2010). Mathematical Foundations of Neuroscience, volume 35 of Interdisciplinary Applied Mathematics. New York, NY: Springer.

Google Scholar

Fitzgerald, T. H., Valentin, A., Selway, R., and Richardson, M. P. (2013). Cross-frequency coupling within and between the human thalamus and neocortex. Front. Hum. Neurosci. 7:84. doi: 10.3389/fnhum.2013.00084

PubMed Abstract | Crossref Full Text | Google Scholar

FitzHugh, R. (1955). Mathematical models of threshold phenomena in the nerve membrane. Bullet. Mathemat. Biopyhs. 17, 257–278. doi: 10.1007/BF02477753

Crossref Full Text | Google Scholar

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

Google Scholar

Gwon, D., and Ahn, M. (2021). Alpha and high gamma phase amplitude coupling during motor imagery and weighted cross-frequency coupling to extract discriminative cross-frequency patterns. Neuroimage 240:118403. doi: 10.1016/j.neuroimage.2021.118403

PubMed Abstract | Crossref Full Text | Google Scholar

Hasegawa, H. (2000). Responses of a Hodgkin-Huxley neuron to various types of spike-train inputs. Phys. Rev. E 61, 718–726. doi: 10.1103/PhysRevE.61.718

Crossref Full Text | Google Scholar

Hill, A. V. (1936). Excitation and accommodation in nerve. Proc. Royal Soc. London 119, 305–355. doi: 10.1098/rspb.1936.0012

Crossref Full Text | Google Scholar

Hodgkin, A. L. (1948). The local electric changes associated with repetitive action in a non-medullated axon. J. Physiol. 107, 165–181. doi: 10.1113/jphysiol.1948.sp004260

PubMed Abstract | Crossref Full Text | Google Scholar

Hodgkin, A. L. (1958). The croonian lecture - ionic movements and electrical activity in giant nerve fibres. Proc. Royal Soc. London 148:1–37. doi: 10.1098/rspb.1958.0001

PubMed Abstract | Crossref Full Text | Google Scholar

Hodgkin, A. L., and Huxley, A. F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. 117, 500–544. doi: 10.1113/jphysiol.1952.sp004764

PubMed Abstract | Crossref Full Text | Google Scholar

Hoppensteadt, F. C., and Izhikevich, E. M. (1997). Weakly connected neural networks, volume 126 of Applied Mathematical Sciences. Cham: Springer.

Google Scholar

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

PubMed Abstract | Crossref Full Text | Google Scholar

Isler, J. R., Grieve, P. G., Czernochowski, D., Stark, R. I., and Friedman, D. (2008). Cross-frequency phase coupling of brain rhythms during the orienting response. Brain Res. 1232, 163–172. doi: 10.1016/j.brainres.2008.07.030

PubMed Abstract | Crossref Full Text | Google Scholar

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

PubMed Abstract | Crossref Full Text | Google Scholar

Jirsa, V., and Müller, V. (2013). Cross-frequency coupling in real and virtual brain networks. Front. Comput. Neurosci. 7:78. doi: 10.3389/fncom.2013.00078

PubMed Abstract | Crossref Full Text | Google Scholar

Jirsa, V. K., and Haken, H. (1996). Field theory of electromagnetic brain activity. Phys. Rev. Lett. 77, 960–963. doi: 10.1103/PhysRevLett.77.960

PubMed Abstract | Crossref Full Text | Google Scholar

Jirsa, V. K., and Haken, H. (1997). A derivation of a macroscopic field theory of the brain from the quasi-microscopic neural dynamics. Physica D 99, 503–526. doi: 10.1016/S0167-2789(96)00166-2

Crossref Full Text | Google Scholar

Kadmon, J., and Sompolinsky, H. (2015). Transition to chaos in random neuronal networks. Phys. Rev. X 5:041030. doi: 10.1103/PhysRevX.5.041030

Crossref Full Text | Google Scholar

Kepler, T. B., Abbott, L. F., and Marder, E. (1992). Reduction of conductance-based neuron models. Biol. Cybern. 66, 381–387. doi: 10.1007/BF00197717

PubMed Abstract | Crossref Full Text | Google Scholar

Koch, C., Bernander, O., and Douglas, R. J. (1995). Do neurons have a voltage or a current threshold for action potential initiation? J. Comput. Neurosci. 2, 63–82. doi: 10.1007/BF00962708

PubMed Abstract | Crossref Full Text | Google Scholar

Kovach, K., Oya, H., and Kawasaki, H. (2018). The bispectrum and its relationship to phase-amplitude coupling. Neuroimage 173, 518–539. doi: 10.1016/j.neuroimage.2018.02.033

PubMed Abstract | Crossref Full Text | Google Scholar

Kuramoto, Y. (1975). Self-entrainment of a population of coupled non-linear oscillators. In International Symposium on Mathematical Problems in Theoretical Physics, pages 420-422. Lecture Springer Notes in Phys. 39:13365. doi: 10.1007/BFb0013365

Crossref Full Text | Google Scholar

Kuramoto, Y. (1984). “Chemical oscillations, waves, and turbulence,” in Springer Series in Synergetics (Cham: Springer-Verlag).

Google Scholar

L‘vov, V. S. (1998). Nonlinear physics: Universality of turbulence. Nature 396, 519–521. doi: 10.1038/24998

Crossref Full Text | Google Scholar

Ladenbauer, J., Ladenbauer, J., Külzow, N., de Boor, R., Avramova, E., Grittner, U., et al. (2017). Promoting sleep oscillations and their functional coupling by transcranial stimulation enhances memory consolidation in mild cognitive impairment. J. Neurosci. 37, 7111–7124. doi: 10.1523/JNEUROSCI.0260-17.2017

PubMed Abstract | Crossref Full Text | Google Scholar

Lapicque, L. (1907). Recherches quantitatives sur l'excitation electrique des nerfs traitée comme une polarization. J. Physiol. Pathol. Gen. 9, 620–635.

Google Scholar

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

PubMed Abstract | Crossref Full Text | Google Scholar

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

PubMed Abstract | Crossref Full Text | Google Scholar

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

PubMed Abstract | Crossref Full Text | Google Scholar

Mainen, Z. F., Joerges, J., Huguenard, J. R., and Sejnowski, T. J. (1995). A model of spike initiation in neocortical pyramidal neurons. Neuron 15, 1427–1439. doi: 10.1016/0896-6273(95)90020-9

PubMed Abstract | Crossref Full Text | Google Scholar

Marder, E., and Bucher, D. (2001). Central pattern generators and the control of rhythmic movements. Curr. Biol. 11, R986–R996. doi: 10.1016/S0960-9822(01)00581-4

Crossref Full Text | Google Scholar

Marshall, L., Helgadóttir, H., Mölle, M., and Born, J. (2006). Boosting slow oscillations during sleep potentiates memory. Nature 444, 610–613. doi: 10.1038/nature05278

PubMed Abstract | Crossref Full Text | Google Scholar

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

PubMed Abstract | Crossref Full Text | Google Scholar

Muller, L., Chavane, F., Reynolds, J., and Sejnowski, T. J. (2018). Cortical travelling waves: mechanisms and computational principles. Nat. Rev. Neurosci. 19, 255–268. doi: 10.1038/nrn.2018.20

PubMed Abstract | Crossref Full Text | Google Scholar

Munia, T. T. K., and Aviyente, S. (2019). Time-frequency based phase-amplitude coupling measure for neuronal oscillations. Sci. Rep. 9:12441. doi: 10.1038/s41598-019-48870-2

PubMed Abstract | Crossref Full Text | Google Scholar

Nakao, H. (2016). Phase reduction approach to synchronisation of nonlinear oscillators. Contemp. Phys. 57, 188–214. doi: 10.1080/00107514.2015.1094987

Crossref Full Text | Google Scholar

Palkar, G., Wu, J., and Ermentrout, B. (2023). The inhibitory control of traveling waves in cortical networks. PLoS Comput. Biol. 19, 1–32. doi: 10.1371/journal.pcbi.1010697

PubMed Abstract | Crossref Full Text | Google Scholar

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

PubMed Abstract | Crossref Full Text | Google Scholar

Proment, D., Nazarenko, S., and Onorato, M. (2009). Quantum turbulence cascades in the Gross-Pitaevskii model. Phys. Rev. A 80, 1–4. doi: 10.1103/PhysRevA.80.051603

Crossref Full Text | Google Scholar

Qin, Y., Menara, T., Bassett, D. S., and Pasqualetti, F. (2021). Phase-amplitude coupling in neuronal oscillator networks. Phys. Rev. Res. 3:023218. doi: 10.1103/PhysRevResearch.3.023218

Crossref Full Text | Google Scholar

Qin, Y., and Sheremet, A. (2022). Mesoscopic collective activity in excitatory neural fields, part I: governing equations. arXiv [preprint] arXiv:2206.07966. doi: 10.48550/arXiv.2206.07966

Crossref Full Text | Google Scholar

Rasch, B., and Born, J. (2013). About sleep's role in memory. Physiol. Rev. 93, 681–766. doi: 10.1152/physrev.00032.2012

PubMed Abstract | Crossref Full Text | Google Scholar

Rinzel, J., and Ermentrout, B. (1998). “Analysis of neural excitability and oscillations,” in Methods in Neuronal Modelling: from Synapses to Networks, C. Koch, and I. Segev (Cambridge, MA: MIT Press), 251–291.

Google Scholar

Rosenblum, M., Tass, P., Kurths, J., Volkmann, J., Schnitzler, A., and Freund, H.-J. (2000). Detection of phase locking from noisy data: application to magnetoencephalography. In Chaos In Brain? World Scientific. 4, 34–51. doi: 10.1142/9789812793782_0004

Crossref Full Text | Google Scholar

Roxin, A., Brunel, N., Hansel, D., Mongillo, G., and van Vreeswijk, C. (2011). On the distribution of firing rates in networks of cortical neurons. J. Neurosci. 31, 16217–16226. doi: 10.1523/JNEUROSCI.1677-11.2011

PubMed Abstract | Crossref Full Text | Google Scholar

Sheremet, A., Burke, S. N., and Maurer, A. P. (2016). Movement enhances the nonlinearity of hippocampal theta. J. Neurosci. 36, 4218–4230. doi: 10.1523/JNEUROSCI.3564-15.2016

PubMed Abstract | Crossref Full Text | Google Scholar

Sheremet, A., Kennedy, J. P., Qin, Y., Zhou, Y., Lovett, S. D., Burke, S. N., et al. (2018a). Theta-gamma cascades and running speed. J. Neurophysiol. 121, 444–458. doi: 10.1101/422246

Crossref Full Text | Google Scholar

Sheremet, A., Qin, Y., Kennedy, J. P., Zhou, Y., and Maurer, A. P. (2019). Wave turbulence and energy cascade in the hippocampus. Front. Syst. Neurosci. 12:62. doi: 10.3389/fnsys.2018.00062

PubMed Abstract | Crossref Full Text | Google Scholar

Sheremet, A., Zhou, Y., Kennedy, J., Qin, Y., Burke, S., and Maurer, A. (2018b). Theta-gamma coupling: a nonlinear dynamical model. BioRXiv. doi: 10.1101/304238

Crossref Full Text | Google Scholar

Sheremet, A., Zhou, Y., Qin, Y., Kennedy, J. P., Lovett, S. D., and Maurer, A. P. (2020). An investigation into the nonlinear coupling between CA1 layers and the dentate gyrus. Behav. Neurosci. 134, 491–515. doi: 10.1037/bne0000366

PubMed Abstract | Crossref Full Text | Google Scholar

Skaggs, W. E., McNaughton, B. L., Permenter, M., Archibeque, M., Vogt, J., Amaral, D. G., et al. (2007). EEG sharp waves and sparse ensemble unit activity in the macaque hippocampus. J. Neurophysiol. 98, 898–910. doi: 10.1152/jn.00401.2007

PubMed Abstract | Crossref Full Text | Google Scholar

Soltesz, I., and Deschênes, M. (1993). Low- and high-frequency membrane potential oscillations during theta activity in CA1 and CA3 pyramidal neurons of the rat hippocampus under ketamine-xylazine anesthesia. J. Neurophysiol. 70, 97–116. doi: 10.1152/jn.1993.70.1.97

PubMed Abstract | Crossref Full Text | Google Scholar

Staresina, B. P., Bergmann, T. O., Bonnefond, M., Van Der Meij, R., Jensen, O., Deuker, L., et al. (2015). Hierarchical nesting of slow oscillations, spindles and ripples in the human hippocampus during sleep. Nat. Neurosci. 18:1679. doi: 10.1038/nn.4119

PubMed Abstract | Crossref Full Text | Google Scholar

Stiefel, K. M., and Ermentrout, G. B. (2016). Neurons as oscillators. J. Neurophysiol. 116, 2950–2960. doi: 10.1152/jn.00525.2015

PubMed Abstract | Crossref Full Text | Google Scholar

Stimberg, M., Brette, R., and Goodman, D. F. M. (2019). Brian 2, an intuitive and efficient neural simulator. Elife 8:e47314. doi: 10.7554/eLife.47314

PubMed Abstract | Crossref Full Text | Google Scholar

Strogatz, S. H. (2000). From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators. Physica D 143, 1–20. doi: 10.1016/S0167-2789(00)00094-4

Crossref Full Text | Google Scholar

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

PubMed Abstract | Crossref Full Text | Google Scholar

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

PubMed Abstract | Crossref Full Text | Google Scholar

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

PubMed Abstract | Crossref Full Text | Google Scholar

Zakharov, V., L'Vov, V., and Falkovich, G. (1992). Kolmogorov Spectra of Turbulence 1. Wave Turbulence. Berlin: Springer.

Google Scholar