Analytical insights on theta-gamma coupled neural oscillators.

In this paper, we study the dynamics of a quadratic integrate-and-fire neuron, spiking in the gamma (30-100 Hz) range, coupled to a delta/theta frequency (1-8 Hz) neural oscillator. Using analytical and semianalytical methods, we were able to derive characteristic spiking times for the system in two distinct regimes (depending on parameter values): one regime where the gamma neuron is intrinsically oscillating in the absence of theta input, and a second one in which gamma spiking is directly gated by theta input, i.e., windows of gamma activity alternate with silence periods depending on the underlying theta phase. In the former case, we transform the equations such that the system becomes analogous to the Mathieu differential equation. By solving this equation, we can compute numerically the time to the first gamma spike, and then use singular perturbation theory to find successive spike times. On the other hand, in the excitable condition, we make direct use of singular perturbation theory to obtain an approximation of the time to first gamma spike, and then extend the result to calculate ensuing gamma spikes in a recursive fashion. We thereby give explicit formulas for the onset and offset of gamma spike burst during a theta cycle, and provide an estimation of the total number of spikes per theta cycle both for excitable and oscillator regimes.

1 Introduction

Oscillations of neural activity are ubiquitous in the brain in many frequency bands [1], and it has been often argued that they play a functional role in cortical processing [2-4]. Physiological experiments and computational models have shown that ongoing brain oscillations are involved in sensory-motor functions [5], synaptic plasticity [6], memory formation and maintenance [7], among many other cognitive tasks. Indeed, it has been reported [2] that intrinsic brain rhythms can bias input selection, temporally link neurons into assemblies, and facilitate mechanisms that cooperatively support temporal representation and long-term consolidation of information. Notably gamma oscillations (>30 Hz) are prominent in neocortex during attention [8], sensory processing [9,10], or motor control tasks [11], together with slower rhythms in the theta (3–8 Hz) or delta (1–3 Hz) range that have also been linked to various aspects of cognitive processes like working memory or the transmission of sensory and motor signals.

Many recent contributions point to nontrivial interactions among different frequency bands [12-14], such as phase-amplitude [15,16] or phase-phase coupling [17,18] that can facilitate the simultaneous integration of multiple layers of information [19]. The hippocampus is a privileged site for observing such interactions [11,20], since theta and gamma waves are particularly strong and reliable in that region [21]. Another particular case is represented by perception of speech signal performed by auditory cortex. In fact, to capture the many different relevant features of speech (i.e., syllables, vowels, consonants, etc.), the brain must be able to parse the speech signal over these many time-scales at the same time. A number of recent works introduced the hypothesis that a network of nested theta (3–8 Hz) and gamma (30–100 Hz) rhythms could accomplish this task [22-24], given their matching in frequency with syllabic and phonemic time-scale, respectively. Since there is no external onset signaling the presence of an incoming syllabic content, the phase of the gamma rhythm needs to be reset by some intrinsic mechanism, e.g., by theta input [23]. It becomes therefore important to know the time to first spike, which would be a measure of the speed of gamma phase resetting, as well as the time to last spike and the spiking frequency during excitable period.

There is a large literature on mathematical analysis of single frequency oscillators in networks of cortical circuits [25-31], and much work has been done in computational modeling of neural oscillations [2,32,33]. There is also a significant number of mathematical studies on cross-frequency interactions, however, most of that analysis is limited to the cases of weak coupling [34-37]. Strong coupling case has been analyzed either with pulsatile coupling [25,38,39] or with semianalytical and computational techniques [40-42]. Importantly, the question of how strong continuous coupling between slow and fast oscillations influences frequency and time of fast spikes has not been treated analytically, at least to the best of our knowledge. Yet experimental data suggest that phase-amplitude coupling in the brain is continuous (i.e., low-frequency phase is conveyed through local field potential, a continuous signal) and strong [15,39,41], so this will be the regime we aim to study in the present work.

In this article, we provide analytical insights on the precise spiking times of a simplified Pyramidal Interneuron Network Gamma (PING) [41] during theta modulation. Two separate cases are studied: In the first setting, which we will refer to as oscillatory regime, the gamma network behaves as an intrinsic oscillator whose spike frequency is modulated by the theta phase; in the second, named excitable regime, gamma spikes are only evoked when input coming from the theta oscillator is strong enough. In the latter case, the system is in an “excitable” regime, where theta pushes gamma back and forth across a Saddle-Node on Invariant Circle (SNIC) bifurcation. The analysis can be generalized beyond theta-gamma nested oscillations; indeed it describes any coupling between low and a high frequency rhythms [43], provided that the latter is produced through feedback inhibition to the excitatory cell. To compute the time to the first gamma spike, we used different approaches for the two regimes: In the oscillatory case, we reduce the system in order to describe its dynamics with the Mathieu equation [44], and in the excitable case we apply an extension of geometric singular perturbation theory [45-47]. We then use a combination of the two to get successive spike times and an estimation of the total number of spikes per theta cycle.

The paper is organized as follows:

1. In Sect. 1, we introduce the system to be studied.

2. In Sect. 2, we consider the system in the oscillatory regime and compute time to first gamma spike using Mathieu functions. We found that spike time is mainly determined by the magnitude of theta-gamma coupling (λ) and of theta frequency.

3. In Sect. 3, we turn our attention to the excitable regime where theta phase determines the magnitude of input, thereby causing the gamma circuit to spike.

4. Finally, we show that our approach gives results in agreement with direct numerical simulations of the system of interest.

In our analysis, we use tools from geometric singular perturbation theory. This approach normally fails in proximity of nonhyperbolic points, as it would be the case for the system considered in the present paper, but the blow-up method extension provided in [48] allows us to compute approximations of the passage time to the first spike in the excitable case, and it is used both in the oscillator and excitable cases to estimate the duration of inhibition and the passage time of subsequent spikes. The latter estimates are based on the idea that inhibition puts the system in a state of quasi equilibrium; consequently, they work well if inhibition is strong and excitation not too high.

2 Theta-Gamma Coupled Oscillator

We consider a minimal formulation of a theta modulated gamma spiking network. One single Excitatory Gamma (EG) neuron (), modeled as a θ-neuron [49] receives an excitatory input coming from an oscillator (Θ) whose natural frequency lies in the theta band. The canonical θ-neuron model is described by a phase variable lying on a one-dimensional circle in the range , a spike is produced when . The EG neuron participates in a PING rhythm, although in our case the inhibitory gamma neuron is instantaneously enslaved to the excitatory cell, meaning that every excitatory spike would immediately prompt a simultaneous inhibitory spike [32]. This allows us to suppress the explicit dynamics of the inhibitory gamma neuron and focus on inhibitory synaptic dynamics only. Our system can be described by the following equations:

where is the instantaneous phase of the slow rhythm variable (delta/theta frequency band, i.e., 1–8 Hz), which provides the sinusoidal modulatory input to the EG cell; is the variable representing the activation of the inhibitory synapse; represents constant driving input to excitatory gamma neuron; λ is the strength of theta-gamma coupling; is the inhibitory synaptic strength; ω has been chosen so that frequency falls into the theta range; is a scaling parameter that scales inversely with the time constant of synaptic inhibition; is a second, slower, scaling parameter that has been chosen such that , an assumption that is motivated by biophysical considerations and, in addition, keeps the three time scales (theta rhythm, synaptic inhibition, and excitatory membrane potential) well separate.

We will consider two cases: the oscillator case, defined by , and the excitable case, defined by . The characterizing feature of the oscillator setting is that – subsystem in (1) is an intrinsic oscillator at every stage of a Θ-cycle, i.e., the total current input to EG neuron is always positive. In the excitable case, on the other hand, part of theta oscillation period is such that subsystem of (1) has an attracting quasisteady state, i.e., the total input to the EG neuron is negative or positive depending on Θ-oscillator phase. If , the net input to EG neuron is always negative and the gamma circuit is always silent.

3 Time to First Spike, Oscillator Case

Let us consider the case in which constant driving term in system (1) is positive and such that, in absence of theta modulation, the EG neuron would fire periodically with a spiking frequency in the gamma range (30–150 Hz). We assume that the dynamics of the theta oscillator is at least one order of magnitude slower than synaptic decay, so that almost no residual inhibition is present at the beginning of a new theta cycle. In order to obtain an equation in the form of Mathieu equation, we first perform a change of variables in system (1) , going from θ-variable to membrane voltage . As it is known [40], θ-neuron is formally equivalent to the Quadratic Integrate and Fire (QIF) neuron:

The two neural models are formally equivalent if we define the reset conditions as

where the time of spike. We have omitted the synaptic input dynamics since we assume that inhibition is directly enslaved to spikes coming from EG neuron, hence the inhibitory synapse stays inactive up to the first EG spike. We restate the system in (2) as a single equation, assuming by convention that at (or at the beginning of a new theta cycle):

For , this equation has an exact solution in terms of Mathieu functions. This can be found by imposing one more change of variable:

where the prime mark denotes the time derivative (a similar transformation is used in [50] where the cosinusoidal forcing term was replaced by exponential decay, leading to a different solution of the corresponding differential equation). Hence, we write (3) as a second-order differential equation:

If parameters a, q, z are rescaled as following:

then Eq. (5) has the form of a Mathieu equation:

To interpret Eq. (7), we need temporal rescaling from t to z, and as a consequence the period of cosinusoidal term, which in Eqs. (1) and (2) was , becomes . The solutions to Eq. (7) are linear combinations of the even and odd Mathieu functions [44], and , respectively. The solution

obeys the desired initial conditions, where the dot indicates the derivative with respect to z. Because of the change of variable in (4), the spiking times in the absence of inhibition correspond to the zeros of the solution of (7) (Fig. 1). Hence, by scaling back to the original variables and looking at the first zero of we obtain the time to the very first spike . We numerically compute the time to the first spike as a function of parameters a and q, i.e., and λ. The subsequent spikes, on the other hand, depend on inhibition, and thus cannot be described by (3) alone. We looked for solutions of (2) with initial condition and . Figure 2 shows the time to first spike as a function of λ with fixed at three values, , , and . Note that for the dependence on λ is strong, but for larger values of the sensitivity of with respect to λ is smaller, since becomes the dominant input term. In the next section, we will consider the case when inhibition is strong and fully controls the gamma spikes.

Fig. 1

Oscillatory regime—dynamics of theta-modulated EG neuron in absence of inhibition. Red line: membrane voltage of EG neuron in presence of theta modulation (blue line) without inhibitory synaptic input (). Solution to Mathieu equation is plotted in yellow

Fig. 2

Oscillatory regime—time to first spike. Time to first spike as a function of coupling constant λ, for , and . Black line: simulation; red dotted line: solution to Mathieu equation. ,

4 Time to First Spike, Excitable Case

The excitable case implies that , and . Under these assumptions, the gamma spikes are only possible when Θ lies in a proper subinterval of , which corresponds to the values of Θ for which . This ensures that the dynamics of (1) cross the SNIC bifurcation for a certain value of . We carry out the computation with the initial conditions

where is defined by the condition , i.e. is a stationary fixed point in the absence of Θ positive input. The gamma neuron relaxes to when theta modulation is turned off. Note that is not required for our solution to be applicable, since, for any initial value , quickly converges to . At the end of every theta cycle goes back to zero, since its decay constant is one order of magnitude bigger than , and the EG cell has stopped firing once inhibition pushed it below the SNIC bifurcation. We start by computing an estimate of the time to the first gamma spike. System (1) involves two time scales, one that controls the intrinsic dynamics of the EG neuron and the other comes from Θ modulation. In the excitable case, rather than using the approach based on Mathieu functions, we use geometric singular perturbation theory. This approach leads to explicit estimates of the onset and duration of the gamma burst, it gives some geometric insights and can be applied in a more general setting. In order to compute the time at which the fast and the slow dynamics intersect, we need the value of Θ corresponding to , i.e., where the SNIC bifurcation takes place. Simple algebra shows that this occurs when

To ensure that (9) has solutions, we verify that the RHS of (9) is in the interval . For the upper bound, we have

The lower bound is obtained as follows:

given that . Let now be a solution of (9) satisfying and let us consider system (1) with initial conditions (8): it is clear that any trajectory of the system can roughly be divided into two separate chunks (see Fig. 3). Starting from point , the system immediately enters the slow motion part of the trajectory, which is adjacent to the nullcline . The slowest region of motion lies in the vicinity of the singular point , where both and its derivative with respect to are zero. Once the trajectory has gone beyond the singular point, φ turns positive again and grows quadratically in magnitude. This way quickly reaches the value , since it is well known that any unbounded solution of the theta neuron for positive net input explodes in finite time. At the same time Θ increase is of order . Now let us start by computing the time spent along the fiber which is close to the nullcline, and then direct our attention to the motion in the neighborhood of the singular point. The time needed for Θ to grow from −π to equals

Fig. 3

Excitatory regime—phase plane. Phase portrait of system (1) for , , , . The blue line represents the trajectory of the system when starting from initial conditions : it passes along the nullcline (in red) and then quickly escapes to once past the singular point . The dotted purple line shows that, for any starting point where and , the trajectory converges to the blue line

When Θ reaches , is (recall that this is also ) close to the threshold value of . In order to estimate the time that EG neuron needs to produce the first spike, i.e., to reach , we need to examine the behavior when close to point . We first translate the variable Θ, introducing . Using Taylor expansion around and , we can write

We transform (1) to the coordinates , taking into the account the expansion (11) and ignoring which remains zero until the first gamma spike (we omit the tilde for the simplicity of notation). The resulting system is

with .

We then rescale the variables of (12) as follows:

In terms of the rescaled variables, with the tilde omitted, system (12) becomes

This means that the nullcline of system (13) for , defined by the parabola , is a good approximation for the nullcline of system (1) in the neighborhood of the singular point , or equivalently, in coordinates , the point . Thus system (13) has the same form as system (2.5) in [48] and can be restated as a Riccati equation:

By performing a change of variable, it can be shown that (14) is equivalent to a second-order Bessel equation. The following function is a general solution of (14):

where are Bessel functions of the first kind of order ν. The only solution approaching the left branch of the nullcline parabola for is the one obtained by choosing , thus we pick this value of c. The inverse of function , namely , defines the trajectory of x as a function of y (Fig. 3). Unfortunately, due to its highly nonlinear form, it is impossible to compute directly.

We use these results together with Theorem 2.1 and Remark 2.11 in [48] in order to derive the following estimate (the result dates back to much earlier, see for example [51]).

Proposition 1Letandsatisfy. Also fix. Consider a family of solutions of (13) with initial conditionsand. Letbe the intersection point of this trajectory with the line. Then, for sufficiently smallδ,

wherestands for the smallest positive zero of the denominator in (15):

From now on, we will use the numerical approximation

Note that the solution with initial conditions (8), transformed to the coordinates , satisfies the assumptions of Proposition 1. Therefore, estimate (17) holds.

Now let be the time the of the first gamma spike, i.e., when . From (16), it is easy to see that, after scaling back to the original variables (, Θ, ), can be written as

with

The term in (17) and the following one, of order in , happen to be zero in the theta neuron model (as well as in the QIF model) when there is no excitatory feedback from the EG cell to the theta band oscillator (see the Appendix). The next nonzero term in (17) is then of order , which represents the error with respect to the time at which the true trajectory of the system reaches . The value of δ does not have to be small, on the contrary our approximation works better when δ is such that the trajectory of the system is close to the asymptote , as it is the case for the EG cell spike threshold .

Predictions of Proposition 1 are illustrated in Figs. 3 and 4.

Fig. 4

Excitable regime—time to first spike. Time to first spike as a function of coupling constant λ, for , and . Black full line: simulation; red dotted line: analytic solution. ,

5 Subsequent Gamma Spikes, Oscillator Case

In the oscillator case, we assume that inhibition is strong enough to push the system below the SNIC bifurcation, regardless of the value of Θ, i.e.,

If the opposite is true, the system does not encounter the bifurcation, since is always greater than zero, and our analysis cannot be applied to subsequent spikes. We wish to derive an estimate on the number of EG spikes occurring along one Θ period, which is given by the time needed for Θ to grow from −π to π, thus equal to . Let be the subsequent gamma spikes and , the corresponding values of Θ. Let be the relative time after at which the total driving input to the EG neuron reaches zero from negative values:

From now on, we use the fact that , and relabel . Hence, we can write

We expect to be of order from (20), is then large compared to . We can thus replace (20) by a simpler formula:

Now

We denote the time interval between two successive gamma spikes by and use the following estimate:

where

Estimate (24) is obtained analogously as (17). We can write the modulated instantaneous Interspike Interval (ISI), i.e., the instantaneous period, as

We derive the intrinsic period of PING in absence of any modulation:

After some algebra and performing a second-order Taylor expansion around , i.e., when theta excitation is minimal, the ISI becomes

The smallest ISI is obtained by expanding around , i.e., when theta excitation is maximal:

from (27) and (28) we can estimate respectively the lowest and highest gamma frequencies attained during theta modulation. We then derive an expression for the interval of time ΔT between the first to the last gamma spike within half a Θ-period:

where is the number of gamma spikes as Θ varies between −π and 0. Writing

and

we obtain

As provides an estimate for the number of gamma spikes as Θ grows from −π to 0, the total number of gamma spikes M is the largest integer such that

This formula works well, especially when inhibition is sufficiently strong. Figure 5 shows two cases where the formula gives the exact prediction of the number of gamma spikes and a good approximation of spike times. It is worth to mention that in the oscillatory case there is no phase reset at the end of a theta cycle, meaning that the initial conditions are never the same at the beginning of a theta oscillation. As a consequence, the result in (31) does not hold as a rigorous solution but as an average estimate, and the exact number of spikes can still vary over different trials. In Fig. 6, we show the direct comparison between the predictions of the formula and the simulation, as a function of λ. Note that for λ very small the estimate of the formula is too big. There we would need to include more terms in the ε expansion to get a more accurate prediction. When λ is large, for fixed , the positive input is such that inhibition is not sufficient to periodically time the spikes. As a consequence, the estimate of formula (31) becomes too small.

Fig. 5

Oscillatory regime—dynamics. Two cases in which formula (31) gives a correct prediction of the number of gamma spikes and a fair estimate of spike times. Upper panel: instantaneous firing frequency of gamma cell obtained from simulation (full red line) and from Eq. (25) (black dotted line). Lower panel: the simulation of EG cell membrane potential is shown in red while black dotted lines represent firing times predicted by our analysis; the blue curve shows theta modulation (). Left: , , , , . Right: , , , ,

Fig. 6

Oscillatory regime—number of spikes. Number of spikes in the simulation (blue) and the prediction of the formula as a function of λ varying from 0 to , . In the left panel, , . In the right panel, ,

6 Subsequent Gamma Spikes, Excitable Case

6.1 Second Gamma Spike

Let , be defined by

Note that

but, similarly to the oscillator case, we expect to be of order from (32). This allows us to neglect the last two terms on the RHS of (33). Hence, by (32) and (33), we have

Further, we write

and substitute into (34) getting

It then follows from (35) that keeping only the leading terms:

6.2 Subsequent Gamma Spikes

We can write a variant of estimate (31) for the excitable case:

where now the extremes in the integrals are chosen to be the times of the first and last gamma spikes (i.e., the times when the EG neuron crosses the SNIC bifurcation respectively from below and above), assuming that these would be approximately symmetric with respect to the Θ cycle.

This formula works adequately for large inhibition and relatively small (negative) . Otherwise, due to the intricate interplay between the growth of Θ and the decay of s almost to 0 (witnessed in the computation of ), it is not sufficient to have just the lowest terms of the ε expansion of . Figure 7 shows two cases where the formula gives the exact prediction of the number of gamma spikes. In Fig. 8, we show the direct comparison between the predictions of the formula and the simulation, as a function of λ. For λ large inhibition is too weak to time the spikes and the estimate of the formula becomes too small.

Fig. 7

Excitable regime—dynamics. Three cases for which formula (31) gives a correct prediction of the number of gamma spikes. Plot colors as in Fig. 4. and , i.e. the theta phases where the first and last Hopf bifurcation approximately take place, are shown in cyan. , and . Left panel: , . Right panel: ,

Fig. 8

Excitable regime—number of spikes. Number of spikes in the simulation (blue) and the prediction of the formula as a function of λ. Left panel: , , , . Right panel: , , ,

7 Conclusions and Future Directions

In this paper, we investigated how a continuous, strong, low frequency (1–10 Hz) modulation determines the spiking properties of a simplified PING oscillator. This work has been particularly motivated by recent investigation on the role of theta-gamma interactions in processing speech signals [52]. Syllabic input are in fact known to possess a quasiperiodic structure matching theta frequency [24]. Within this framework, theta-modulated gamma spikes need to be aligned to the onset and the offset of linguistically relevant chunks [23]. It is then crucial to understand the timing of gamma spikes and the way they are influenced by theta input, since theta is supposed to detect the presence of long timescale syllabic content. It remains to be unveiled whether the scaling we analytically determined here is produced in more realistic models for speech processing [52] currently under development. Indeed, this result could also be used for other purposes: investigating how theta fluctuations modulate gamma firing in the hippocampus; determining the impact of alpha oscillations on higher frequencies (including gamma), which are thought to carry bottom-up information in visual perception. Indeed timing of first spike is assumed to be particularly relevant in visual cortex, since it is has been shown that it would facilitate the neural encoding of stimuli [53].

To explore the dynamics of the system, we split the problem into two parameter regimes: In the first, the frequency of gamma spikes is only modulated by theta phase, while in the second the gamma cell would only fire if forced by theta input. In the former regime, by restating the problem in form of a Mathieu differential equation and looking at the first zero of the Mathieu function solving the initial value problem, we were able to find the time to first gamma spike. In the latter, we separate the dynamics into three time scales, one characterizing EG neuron dynamics in absence of any external input, one for theta dynamics, and one for synaptic inhibition, and we approximate the time to first spike by using an extension of the geometric singular perturbation theory based on the application of the blow-up method [46,48].

Computations align with the intuition (arising from the fact that is a type I neuron) that time to first spike decreases in both cases with coupling strength λ and constant driving current . Interestingly, in the excitable case, we found that time to first spike depends approximately on , which implies that it saturates rapidly as λ grows. As a second notable result, scales as , being the speed of theta cycle. Building on these results, we subsequently computed the time to successive spikes in both regimes, where inhibitory synaptic decay time becomes an important factor. For both regimes were able to compute approximate spike times and predict the exact number of spikes per theta cycle (and instantaneous frequency of firing as a direct consequence) in a range of parameter values that leads to firing within the gamma frequency band.

In the present work, we analyzed a simple system in which coupling was limited to a feedforward theta-gamma connection. It would be a natural next step to extend the analysis to bidirectional coupling by including a feedback from gamma spikes to the Θ oscillator. A second assumption we made in constructing our system stated that the gamma circuit internal delay between excitatory and inhibitory spikes was negligible, meaning that both cells would fire at exactly the same time. To make the model more biologically appealing, one could relax this hypothesis by introducing a synaptic delay after an excitatory spike and study the correspondent system (i.e., a full PING). For relatively short delays, we would expect the results obtained in this paper to hold at least qualitatively. Throughout this paper, we considered gamma to be a simplified PING generator, on the other hand it still remains an open question whether the same characteristics of theta-gamma modulation we explored here would still be found in a different gamma generator, e.g., an Interneuron Network Gamma (ING) network [54] that can still be implement with Type I neurons as in the case of this work.

Appendix

We show that the term of order in expansion (17) is zero in our model. The subsequent term, of order is also zero, but we do not include the result here since computations are long and heavy. The interested reader could derive this result from [51].

In the absence of any synaptic input from the inhibitory neuron, we restate system (1) to reduce the complexity of notation:

Then, following [51], we write the expansion

which approximates the time to the first spike in the excitable case.

Coefficient in (41) might be written as

where S stands for the coordinates of the singular point and subscripts indicate the derivatives, i.e., is the first derivative of φ with respect to x, taken at . It is easy to verify that any derivative of φ with respect to x of order n, for n odd, is equal to zero at S. Furthermore, since is constant in system (40), is clearly zero. Hence, .

Competing Interests

The authors declare that they have no competing interests.

Authors’ Contributions

LF, MK, AH, BG designed research; LF, MK performed research; LF, MK, AH, BG wrote the manuscript.