For periodic signals, Fourier coefficients provide the right mathematical tool for extracting the main frequencies. To capture the pseudo-periodicity of gamma rhythms, neuroscientists have used the following computational tool (Henrie and Shapley, 2005; Chariker et al., 2018).

The idea is to fix a time interval of suitable length T and to compute Fourier coefficients on [0, T] as if the signal was periodic with period T, i.e., let

and define the power concentrated at frequency k to be |f^(k)|2. To capture gamma-band frequencies, the time interval T is usually chosen to be between 100 and 500 ms: too short of an interval will fail to capture the relevant frequencies, and too long of an interval is ineffective since the signal is not truly periodic.

For a fixed initial time t₀, we perform the computation above repeated with many time shifts, i.e., we sample the signal on [t, t+T] for t = t₀, t₀+dt, t₀+2dt, ⋯  for small dt, and the squares of the computed Fourier coefficients, |f^(k)|2, are averaged over all of the samples to obtain what is called the power spectral density (PSD) and often referred as the Welch method in computational software, refer to Welch (1967).

We consider here the system

which in polar coordinates is

In this system, the unit circle has period 2π, and it attracts all the trajectories distinct from the origin. To put the rhythm in the gamma range, we slow time by a factor of 3 to obtain a frequency of ~50 Hz. Adding white noise, we obtain the following system of stochastic differential equations (SDE):

where (B¹, B²) denotes a two dimensional Wiener process.

Typical solutions of Equation (4)- with μ_x = μ = μ_y at various levels of noise are shown in Figure 1. The traces are plots of the y-coordinate; the wriggly curves are sample paths of the SDE over many cycles. PSDs at corresponding noise levels are also shown.

At μ = 0.2, the peaks of the y-trace are somewhat irregularly spaced thanks to the noise, but the rhythm is too regular to resemble gamma rhythms produced by the real cortex. Increasing μ to 0.3 increases the amount of variability, but the rhythm is still too regular; in particular, it does not degenerate as can be seen by the empty spot near the origin in the phase plane trajectory. At μ = 0.4, the noise has significantly broadened the PSD (which is desirable), but local properties of Brownian paths also begin to manifest in the y-traces in the form of short rises and falls occurring at rapid successions. Such high-frequency oscillations on top of the main gamma rhythm are not typical of the behavior of membrane potentials in gamma activity. At μ = 0.5, the range of signal frequency becomes a little too broad, threatening to obstruct the main rhythm.

Our conclusion from the study above is that the addition of noise to purely periodic dynamics produces variability that goes some distance toward simulating the irregular character of gamma rhythms, but the use of a single parameter, namely μ= noise level, is too rigid: it is not irregular enough at low noise-levels and introduces undesirable features when noise-level is tuned too high.

Continuing to work with limit cycles subjected to white-noise forcing, we give two examples below to show the beneficial effects of an additional parameter.

Our first example adds shear to the deterministic dynamics. We introduce a new parameter α to Equation (4) to give

The deterministic of this equation, reads in polar coordinates as,

The idea is that the farther a trajectory is from the center, the higher its angular velocity. Figure 2A shows the case of μ = 0.2 with α = 3. It is evident that at this noise level, shear produces significantly more varied behaviors, a fact confirmed by a much broader spectrum.

In a second example, we use anisotropic noise, i.e., the equation is as in (4) with μ_x = 0.5 and μ_y = 0.2. The idea is that the larger μ_x would add variability, while the smaller μ_y would not result in the unwanted Brownian structures in y-traces. These expectations are confirmed in Figure 2B.

Remark 1. We conclude from the study above that because gamma signals are multi-faceted, to properly simulate them one needs to be able to control—independently—their frequencies, amplitudes, and degrees of degeneracy (i.e., the way the rhythm degrades from time to time and re-emerges). There has to be a mechanism for the system to switch from one regime to another at seemingly random times. Limit cycles with noise controlled by one or two parameters can reproduce certain aspects of the signal but do not possess sufficient flexibility. We propose in the sections to follow a system of ODEs with several parameters designed to control directly the properties we think are important and to produce the irregularity of gamma characteristics via stochastic motion in parameter space.

Below u and v represent the absolute values of the E and I-conductances of a typical neuron in a local population. We propose that their dynamics be described by

where a₁, a₂, b, and c are fixed parameters with a1=-0.01, a2=0.1, b=11.9, c=6.6×10-4. These values are chosen to reproduce the nullclines in Figures 3, 4, the dynamical significance of which are explained below. The key parameters are ϵ, γ, and K. They will be discussed at length below; for now, think of them as taking values in

We will adopt the following notation:

The nullclines of Equation (7) are then given by

for the first equation, and

for the second equation. Note that the polynomial f(u) = −K(u−a₁)(u−a₂) reaches its maximum for u = 0.5(a₁+a₂) with a value of 0.25K(a2-a1)2.

Our choice of parameters allows only one intersection between the quadratic function f(u) and the linear function g(u) in the positive quadrant. This is ensured by the condition:

which gives K < 176.8.

There are four stationary points in the positive quadrant, the region of interest. They are

where u^* is the positive solution of

and v^* = bu^*+c. Refer to Figure 1, which gives a sense of the global dynamics and basic structures of the system.

Theorem 1 establishes rigorously the region of interest for this system.

Theorem 1. The positive (u, v)-quadrant is invariant under the dynamics defined by Equation (7), and there exists a bounded absorbing set to which all solutions enter.

Proof. The fact that the positive quadrant is positively invariant follows from the fact that u = 0, v_t = γv(−v+c) and ϵu_t = −u(K(u−a₁)(u−a₂)), v = 0, are solutions lying on the u and v-axes. For the existence of an absorbing set, we compute

By using Young inequality, for u, v > 0,

we can split the right hand side of (8) into a polynomial in u of degree 4, and a polynomial in v of degree 3, with both negative leading coefficients. Since the solution lies in the positive quadrant, by polynomial comparison, we obtain:

where M > 0 and N > 0 can be chosen independently of initial conditions. Integrating (9) leads to the existence of an absorbing set attracting all trajectories.

Remark 2. Equation (7) is inspired by the FitzHugh-Nagumo (FHN) and Leslie-Gower (LG) models. Starting from the FHN system, a well-studied system with a slow-fast structure, we replaced the cubic nullcline with a parabola in the first equation and added factors of u and v respectively in front of the first and second equations to constrain the solutions to the upper-right quadrant. Such modifications were necessary because in the classical FHN system, one variable describes voltage and the other is the so-called recovery variable, but the dynamics of voltage excursions and recovery follow different time courses than E- and I-conductances. Also, both variables in FHN take negative values, which the magnitudes of conductances do not. The variables in the LG system, on the other hand, have been shown (refer to Ambrosio et al., 2018) to have a dynamical character closer to that of the REI mechanism in Chariker et al. (2018). Relaxing the fast constraint (i.e., allowing ϵ to be larger) induces correlated dynamics in u and v that are strikingly similar to those observed in experimental and numerical simulations of conductance dynamics; compare, e.g., Figure 1 in Okun and Lampl (2009) and Figures 2B,E in Chariker et al. (2018).

Recall that we have three parameters: ϵ, γ, and K. We first explain the role of ϵ, fixing for now γ = 1, and studying the dynamics as ϵ is varied for each value of K. The dynamics of the system from ϵ very small to very large for K = 60 (a fairly typical value of K) are summarized in Figure 2. For ϵ ≪ 1, Equation (7) describes a slow-fast system with a limit cycle as can be seen in the top two panels of Figure 2. As ϵ is increased, this limit cycle turns into a sink somewhere between ϵ = 0.1 and ϵ = 1. For large ϵ, e.g., at ϵ = 10, one can show that the critical attractive manifold is the line Δ:v = g(u). All the trajectories reach this line fast and then follow it slowly toward the fixed point (u^*, v^*) as can be seen in the bottom right panel.

Below we will give the analysis for ϵ very small, as well as the Hopf bifurcation that takes the limit cycle to the sink.

Continuing to keep γ = 1, we first examine the dynamics of Equation (7) as functions of K and ϵ. Simulation results are shown in Figure 4. Notice first that these results are consistent with those in Figure 3 with regard to increasing ϵ for fixed K. What is new here is the effect of varying K for each ϵ. Figure 4 shows clearly that larger K corresponds to larger excursions by u and v. This means

Finally, we examine the effect of varying γ. From the equations, it is clear that trajectories of Equation (7) will trace out the same curves as long as γϵ remains constant; and that varying γ keeping γϵ fixed corresponds to changing the speed with which one moves along these curves. For example, at K = 60, for values of ϵ = 0.1 and γ = 1, numerical simulation gives a limit cycle with period ~44 ms (equivalently frequency around 22 Hz). For ϵ = 0.01, γ = 10, the period becomes 4.4 ms (frequency around 225 Hz).

Proposition 2. For each fixed K, the curves traced out by the trajectories of Equation (7) depend only on ϵγ. Fixing K and ϵγ, and varying γ, velocities are proportional to γ; in particular, the frequency of the limit cycle is proportional to γ⁻¹.

The meaning and main general effects of the variation of parameters ϵ, K, and γ in Equation (7) can be summarized as follows:

As we will show momentarily, these are the parameters we need to vary to produce the irregularities seen in gamma rhythms.

The first set of experimental results we used to challenge our model was that reported in Xing et al. (2012). In this article, the authors studied gamma rhythms in awake vs. anesthetized monkeys. LFP from V1 (the primary visual cortex) in response to high contrast sinusoidal grating patches were recorded and the resulting data was analyzed. After the initial power increase (which we do not model), peak gamma frequency was found to be about of 60 Hz in the awake and 40 Hz in the anesthetized monkeys studied (Figure 1 of Xing et al., 2012). Time frequency analysis of single trials confirmed that oscillations in the awake animal were faster (Figure 2). Also, for anesthetized monkeys, the signal was weaker, having an amplitude about 60% that of the awake (Figure 1). In addition to the PSD, which describes spectral power averaged over time, the authors of Xing et al. (2012) (Figure 2) studied temporal structures of gamma-band activity. They found intermittent bursts of activity lasting for small fractions of a second suggesting some short term history dependence.

As proof of concept, we challenged our model to produce two signals with the characteristics of brain rhythms of awake and anesthetized monkeys as reported in the experimental article above. Our goal is not a perfect match with data but to demonstrate how these quantities can be manipulated through parameter selection in a model like the one presented in Sections 3 and 4. Time traces, PSDs, and spectrograms (Fourier power computed on shifted time intervals and plotted as a function of time, refer to Section 2.1) for two signals intended to simulate these two very different cortical states are presented in Figure 8; the exact parameters used are given in the legend. The amplitudes and peak frequencies of the PSD plots are in agreement with the data and the spectrogram shows small bursts of activity.

It is well-known to visual neuroscientists that in primate contrast response, peak frequency and spectral power increase with contrast (Henrie and Shapley, 2005). For definiteness, we challenged the model to reproduce characteristics of the signals in Jia et al. (2013). Figure 2C of Jia et al. (2013) shows that in response to a large grating (10° in diameter), peak frequency increased from about 30 to 44 Hz, and stimulus-induced gamma power was enhanced more than 3-fold as the contrast was increased from about 6 to 50%. Selecting suitable parameters from our model, we were able to build two signals having these spectral characteristics. They are shown in Figure 9A.

It is shown by Brunet et al. (2014) that the repeated presentation of a visual grating stimulus to monkeys resulted in a steady increase of visually induced gamma-band activity in V1 and V4 (and in the synchronization of the two rhythms in these two areas). The authors proposed this as a plausible way to maintain effective stimulus signaling in the face of dwindling firing rates, presumably due to adaptation. Figure 1 in Brunet et al. (2014), which shows LFP traces and power spectra from a recording session with an awake monkey, shows peak frequencies increasing slightly but hovering mostly around 60 Hz, consistent with the values in Example 1. A striking feature here is that the PSDs become increasingly sharply peaked, with significant increases in gamma power at these frequencies with stimulus repetition. In Figure 9B, we present the PSD of 3 signals with these properties generated by our model.

The examples above demonstrate that the two dimensional ODE system with stochastically varying coefficients presented in Sections 3 and 4 of this article is sufficiently flexible that through parameter selection, one can reproduce, on demand, a variety of characteristics observed in gamma-band rhythms in the real cortex.

An early and well-known model is PING (Whittington et al., 2000; Börgers and Kopell, 2003); similar models include (Ermentrout and Kopell, 1998; Tiesinga et al., 2001) among others. These models were the first to use non-linear dynamics to explain gamma rhythms. The original PING model produces highly regular population spikes, capturing successfully the oscillatory behavior of gamma rhythms but not their irregular character. There were several follow-up studies e.g., (Borgers et al., 2005; Börgers, 2017) in which network models were used to produce more nuanced spike patterns.

Another body of study that received much attention is (Brunel and Hakim, 1999; Brunel, 2000). In these studies, the authors started with networks of sparsely coupled integrate-and-fire neurons, and let system size tend to infinity while keeping the number of connections an infinitely small fraction of system size. Arguing that in such a limit distinct neurons are likely to have disjoint sets of presynaptic cells, the authors of Brunel and Hakim (1999); Brunel (2000) modeled neuronal dynamics by an equation consisting of a deterministic part describing meanfield activity plus a Gaussian noise that is independent of the neuron to neuron, and gamma rhythms were modeled as regimes following a supercritical Hopf bifurcation. These are the first reduced models of gamma rhythms that we know of. Another much cited article is by Brunel and Wang (2003). Here, the authors assumed that gamma rhythms consisted of purely periodic motion plus a noise term, and focused on the dependence of the period on various factors.

Experimental studies from the last 15 years brought to light some intriguing features of gamma rhythms produced by the real cortex, stressing their broad-band nature (Henrie and Shapley, 2005). They show that the amplitudes of the oscillations in local field potentials can be quite large (Henrie and Shapley, 2005), far from regimes that emerge following Hopf bifurcations. Another feature of interest revealed by experimental results is that the power and frequencies of the oscillations wander (Xing et al., 2012), with the same patterns often persisting for tens, sometimes up to 200 ms, indicative of some form of short-term memory.

More recently, a number of detailed biological network models have appeared showing that gamma rhythms with the properties above occur naturally as a consequence of Excitatory-Inhibitory interaction in local neuronal processes, calling attention to concepts such as multiple firing events (MFE) (Rangan and Young, 2013) and recurrent excitation-inhibition (REI) (Chariker et al., 2018). The analysis in these recent articles provided an understanding, at least qualitatively, of how both the oscillatory and irregular characters of gamma rhythms come about.

This article is concerned not with physiological processes associated with gamma rhythms but with the mathematical properties of the signal itself. As these rhythms are naturally produced and have very intriguing signatures, we sought to understand how they can be produced using reduced neuronal models defined by low dimensional dynamical systems. As noted above, reduced models of gamma rhythms have been studied before, but how well they capture the broad-band, irregular nature of gamma rhythms had not been evaluated up until now.

Thus, we began in Section 2 with a study of the main group of reduced models in the literature, namely those described by periodic dynamics perturbed by white noise (Brunel, 2000; Brunel and Wang, 2003). We found that these models do reproduce some aspects of the irregular side of gamma behavior but the multi-faceted nature of these behaviors required more multi-dimensional control. To that end, we proposed, in Sections 3 and 4, a model inspired by the FitzHugh-Nagumo system. We do not claim that the model we proposed is the only viable model, far from it, but the following properties of this model are of note: As parameters are varied, the dynamical regimes described range from stable equilibria to Hopf bifurcations to robust limit cycles. Other parameters offer direct control of the amplitudes and frequencies of the periodic regimes. We found also that varying parameter randomly but continuously, such as random walks in parameter space, capture more realistically the wandering nature of gamma characteristics.

In the final section, to demonstrate versatility we challenged the model to reproduce several sets of experimental data, and the model performed satisfactorily. We remark that the modeling approach employed here can be used to study other rhythms provided that we adjust augment the model to accommodate the characteristics of the rhythm in question. A challenging example is signals with superposition of rhythms. In this case, an idea would be to consider two coupled models of ours, the strength of coupling depending on the interaction of these rhythms.