## Abstract

Gamma oscillations of the local field potential are organized by collective dynamics of numerous neurons and have many functional roles in cognition and/or attention. To mathematically and physiologically analyse relationships between individual inhibitory neurons and macroscopic oscillations, we derive a modification of the theta model, which possesses voltage-dependent dynamics with appropriate synaptic interactions. Bifurcation analysis of the corresponding Fokker–Planck equation (FPE) enables us to consider how synaptic interactions organize collective oscillations. We also develop the adjoint method (infinitesimal phase resetting curve) for simultaneous equations consisting of ordinary differential equations representing synaptic dynamics and a partial differential equation for determining the probability distribution of the membrane potential. This method provides a macroscopic phase response function (PRF), which gives insights into how it is modulated by external perturbation or internal changes of parameters. We investigate the effects of synaptic time constants and shunting inhibition on these gamma oscillations. The sensitivity of rising and decaying time constants is analysed in the oscillatory parameter regions; we find that these sensitivities are not largely dependent on rate of synaptic coupling but, rather, on current and noise intensity. Analyses of shunting inhibition reveal that it can affect both promotion and elimination of gamma oscillations. When the macroscopic oscillation is far from the bifurcation, shunting promotes the gamma oscillations and the PRF becomes flatter as the reversal potential of the synapse increases, indicating the insensitivity of gamma oscillations to perturbations. By contrast, when the macroscopic oscillation is near the bifurcation, shunting eliminates gamma oscillations and a stable firing state appears. More interestingly, under appropriate balance of parameters, two branches of bifurcation are found in our analysis of the FPE. In this case, shunting inhibition can effect both promotion and elimination of the gamma oscillation depending only on the reversal potential.

## 1. Introduction

It is known that various biological systems use collective rhythms that are composed of noisy individual oscillations with their complex interactions [1–4]. Among them, local circuits in the cerebral cortex and hippocampus consist of thousands of interacting excitatory and inhibitory neurons that together often generate macroscopic oscillations called the local field potential (LFP). Many recent experiments have shown that gamma oscillations (30–80 Hz) of the LFP may play functional roles in cognition and/or attention [1,5–7]. Furthermore, inhibitory neurons are known to be a key factor for gamma oscillations [5,8–12], which emerge in the LFP even when individual neurons do not fire at every gamma cycle [13,14]. This phenomenon is sometimes called sparse firing. Because of the difficulty of mathematical treatments of this type of activity, much remains unknown about the relationship between the physiological properties of synaptic interactions and macroscopic gamma oscillations. The relationship between external inputs (or changes of internal state) and macroscopic response also remains unclear.

In the case of a limit cycle oscillator, the phase response function (PRF), which is also called the infinitesimal phase resetting curve [15], is a powerful tool for investigating how such rhythms are modulated by external interactions [16,17]. With respect to the PRF for population activity, theoretical studies have used the Fokker–Planck equation (FPE) [18–20] to derive macroscopic PRFs and have shed light on the theoretical relationships between the microscopic oscillations and the macroscopic ones. Furthermore, a recent experimental study determined the macroscopic PRF in a local neuronal population in the CA3 region of the hippocampus [21]. Therefore, it is important to develop theoretical frameworks including the derivation of the macroscopic PRFs to achieve an in-depth understanding of the mathematical and physiological relationships between individual neurons or synapses and macroscopic oscillations.

In this study, we focused on a neuron model, the so-called theta model, which has the simplest form as well as equivalent dynamics to the conductance-based Class I neuron models because it is mathematically the normal form of the saddle-node infinite period [15]. We then derived a version of the theta model that incorporates conductance-based synapses, but remains tractable for numerical and mathematical analysis. The advantage of the theta model is that the resulting FPE lies on a circle (periodic boundary conditions) so that there are no discontinuities (as in leaky integrate-and-fire (LIF) models). Furthermore, the continuous nature of the system and its periodicity in phase space allow us to define a related adjoint operator and thus compute macroscopic PRFs. In addition, we can also use numerical continuation to study bifurcations in the model. The model also retains the physiological properties of conductance-based synapses so that properties like shunting inhibition can be analysed.

Next, to obtain macroscopic PRFs for the neuronal population model, the adjoint method was developed for simultaneous equations composed of ordinary differential equations (ODEs) for synaptic dynamics and a partial differential equation (PDE) for determining the probability distribution of the membrane dynamics. After confirming the validity of the proposed analytical frameworks, we focused on the case of an inhibitory population and used the macroscopic PRFs and bifurcation diagrams to investigate how the synaptic rising time constant, synaptic decaying time constant and shunting inhibition affected the macroscopic gamma oscillations. We finally studied locking to periodic stimuli.

## 2. Material and methods

### 2.1. Derivation of the modified theta model

We consider a homogeneous network of *N* = 1000 inhibitory neurons as described by the quadratic integrate-and-fire (QIF) model. The QIF model is a representative model of Class I neurons and is widely used for computational studies [22–24]. Here, *V*^{(i)}(*t*) denotes the potential of the *i*th neuron at time *t*. Then, the dynamics of *V*^{(i)}(*t*) are represented as
2.1where *ξ* represents stochastic fluctuations with and is the magnitude of the fluctuation, *C* = 1 μF cm^{−2} is the membrane capacitance and *I* (μA cm^{−2}) is a constant current. *g*_{L} = 0.1 mS cm^{−2} is the leak conductance and *g*_{syn} is the synaptic conductances. In absence of synapses and external current, the firing threshold is *V*_{T} = −55 mV. That is, when *V*(*t*) crosses the threshold, then the voltage grows faster than exponentially and blows up to infinity in a finite amount of time. At this point, *V*(*t*) is then reset to *−∞* where it will then tend towards *V*_{R}. The resting potential *V*_{R} is −62 mV. The reversal potential of inhibitory synaptic currents is *V*_{syn} = −70 mV.

The dynamics of the GABAergic (GABA_{A}) synaptic conductances are second-order exponentials
2.2where *t*^{(k)} is the firing time of the pre-synaptic cell, *τ*_{r} = 0.5 ms is a rise time and *τ*_{d} = 5 ms is a decay time [25]. To match the peak conductance of GABA on interneurons to 6.2 nS, which is based on the earlier model and experimental studies [25–27], we set 2.9 × 10^{−}^{4} cm^{2} as the area of a neuron [28] and thus obtained the value Under the condition that the synaptic connections between neurons are random with a probability of *p*_{syn}, it is known that the synaptic influence on each neuron can be homogenized as an approximation [25]. Then, the above equation reads as
2.3Here, we introduce *θ* that satisfies the relation
2.4

Then, equation (2.1) reads as 2.5We call this the modified theta model because it is a physiologically precise version of the theta model, which is well known as a reduced model of Class I neurons [15]. It should be noted that the conversion to the modified phase model is valid even under strong noise or synaptic interactions because it is not derived from a phase reduction, which is valid only in the case of weak interactions [11,29,30], but rather by the transformation of variables.

#### 2.1.1. Fokker–Planck equation and macroscopic phase reduction by the adjoint method

We developed an appropriate adjoint method for the population of the neurons, each of which is described by the modified theta model (equations (2.3) and (2.5)). First, we set the current *I* = 2, noise intensity *σ* = 2 and probability of connections *p*_{syn} = 0.2 as nominal values. Physiologically, *I*, *σ* and *p*_{syn} depend on the state of neurons and networks; therefore, they should be treated more flexibly compared with other parameters. The FPE [31,32] of the neurons and synaptic equation are described as
2.6and
2.7where *c*_{1} = 2/(*V*_{T} − *V*_{R}), *c*_{2} = (2*V*_{syn} − *V*_{R} − *V*_{T})/(*V*_{T} − *V*_{R}), *c*_{3} = −1*/ τ_{r}τ_{d}, c_{4} = −(τ_{r} + τ_{d})/τ_{r}τ_{d}, n_{syn} = p_{syn} × N = 0.2 × 1000 = 200 and A(t) is a firing probability of neurons. Because A(t) is obtained by the flux at θ = π [18,33], it satisfies A(t) = g_{L}P(π,t*)/

*C*.

By introducing *G* = (*g*_{1}, *g*_{2})^{T} = (*g*_{syn}, d*g*_{syn}/d*t*)^{T}, the synaptic dynamics can be re-written as
2.8

The adjoint method is a useful analytical method for deriving the PRF from a limit cycle oscillator [15]. Recently, it has been extended to the FPE by Kawamura *et al*. [19,20]. In the case of our population model (equations (2.6) and (2.7)), the adjoint equation and the dual product for the adjoint method should be derived as a combination of the conventional adjoint method for ODE [15] and the adjoint method for PDE, which was recently developed [19,20]. First, we introduce *P*_{0}(*θ*, *t*) and *G*_{0}(*t*) = (*g*_{0}(*t*), d*g*_{0}(*t*)/d*t*)^{T} as the limit cycle solution for *P* and *G*. They are defined in the region of 0 ≤ *t* < *T*_{macro}, where *T*_{macro} is the period of the limit cycle oscillation.

In this article, *t* = 0 is set so that *P*_{0}(*π*, *t* = 0) is the maximum firing probability. *P*_{0} and *G*_{0} satisfy the relations
2.9and
2.10We then decompose *P*(*θ*, *t*) and *G* as *P*(*θ*, *t*) = *P*_{0}(*θ*, *t*) + *Q*(*θ*, *t*) and *G* = *G*_{0} + *H* with *H* = (*h*_{1}, *h*_{2})^{T}, respectively. The linearized equations for *Q*(*θ*, *t*) and *H* are
2.11
2.12

These equations can be rewritten in matrix representation as
2.13where
2.14
2.15
2.16
2.17The phase space of equation (2.13), *C*_{0}, is the direct sum of two subspaces of *C*_{1} and *C*_{2} as
2.18
2.19
2.20Then, their dual spaces are defined as
2.21
2.22
2.23

It can be analytically confirmed that is a solution of equation (2.13).

From equation (2.13), the linear operator for this phase space is given as 2.24

By defining the dual product as
2.25the adjoint operator *L** is determined by the relation
2.26and the adjoint equation is derived as
2.27and
2.28

We need the zero eigensolutions of the adjoint equation, and with the normalizing condition 2.29

In actual calculations, we integrate equations (2.27) and (2.28) backwards in time from arbitrary initial conditions to obtain the zero eigensolution (i.e. and ) of the adjoint equation, because functional components other than the zero eigensolution have positive eigenvalues; therefore, they eventually vanish (in reverse time).^{1}

#### 2.1.2. Macroscopic phase response function

Here, we consider the macroscopic PRF. When the external perturbation is added to the equation of synapses, the macroscopic PRF is simply provided by On the other hand, when the external perturbation is applied to all of the neurons, the macroscopic PRF is obtained as
2.30where *u*(*θ*) = *c*_{1}(1 + cos*θ*)/*C* is the sensitivity of a single neuron to the phase *θ* [20]. This equation provides unique insights into how the external perturbations and/or changes of internal states affect macroscopic properties of gamma oscillations and will be discussed in following sections.

## 3. Results

In order to study the dynamics of large networks of neurons with synaptic coupling, we have extended the so-called theta model (a continuous version of the QIF) to incorporate conductance-based synapses. We did this in order to explore how aspects of synaptic coupling, for example shunting inhibition, affect the ability of inhibitory networks to produce population rhythms. The advantages of the theta model over other simple spiking models are that the model is continuous and certain numerical computations are much easier. For large numbers of neurons with random coupling, we can reduce the system to a single PDE for the firing rate and synaptic activity (see Material and methods). We can thus turn the focus onto the dynamics of a deterministic system for which there are many available tools.

### 3.1. Population dynamics of the modified theta models

The numerical computation results for the population of the modified theta models (equations (2.3) and (2.5)) with nominal values of parameters are shown in figure 1*a*,*c*, where the macroscopic gamma oscillation can be observed. In figure 1*c*, we plot the total synaptic conductance as a macroscopic version of the population dynamics. Figure 1*b* shows an example of the membrane potential of a single neuron, which is obtained by reversal transform of equation (2.4). We can see that in this parameter set, there is sparse firing in which the individual neurons do not fire at every gamma cycle. We also confirmed that the numerical simulation of the corresponding FPE (equations (2.6) and (2.7)) accurately matches the simulated population dynamics (figure 1*d*). Because our model does not incorporate the effect of resetting and a refractory period, an equivalent representation of this model with resetting and a refractory period is discussed in the electronic supplementary material, text S1 and figure S1.

### 3.2. Bifurcation diagrams of gamma oscillation

As we are looking only at an inhibitory network (with no ‘exotic’ currents like the calcium T-current), the network needs tonic drive and feedback inhibition in order to produce population rhythms. Thus, with no synaptic coupling or with reduced drive, the network will produce asynchronous constant in time dynamics. However, as the drive increases and with sufficient negative feedback, a population rhythm emerges through a loss of stability of the constant asynchronous state via a Hopf bifurcation (HB). Thus, we wish to study the stability of the constant asynchronous state as the drive and negative feedback change. The parameter regions for macroscopic gamma oscillations can be identified by bifurcation analyses of the discretized FPEs [31,34].^{1} We chose the applied current, *I*, and the connection probability, *p*_{syn}, as two parameters of interest. We note that *p*_{syn} essentially sets the magnitude of the synaptic connections via the quantity *c*_{5} in the FPE (equations (2.6) and (2.7)). By varying, say, *I* with *p*_{syn} fixed, we detect HBs. We then compute a two-parameter curve of HBs as a function of the two parameters *I* and *p*_{syn}. Figure 2*a*,*d* shows the two-parameter curve for noise amplitudes, *σ* = 1, 2. Below the curve (small inputs and sparse connectivity), the constant asynchronous state is stable and above it, there will be oscillations in *g*_{syn} and thus any macroscopic measure of the behaviour. We also see that the numerical computations of the modified theta models (equations (2.3) and (2.5)) follow the bifurcation analyses and the numerical solutions of the FPE. The small fluctuations in both the stable regions (figure 2*b*,*e*) and the gamma oscillatory region (figure 2*c*,*f*) are due to the finite-size effect. We note that the curve of HBs is raised up in the higher noise case which means that larger inputs and more densely connected networks are needed for gamma oscillations to emerge and overcome the stability of the asynchronous state.

### 3.3. Macroscopic phase response function by the adjoint method

The adjoint method provides a rigorous technique for the determination of the PRF and has been successfully applied to individual neural oscillators. The PRF of an oscillating system measures the shift in timing of the oscillation owing to a stimulus delivered at a specified phase in the ongoing oscillation. The solution to the adjoint equation is proportional to the PRF when the perturbations are small and, thus, the PRF is useful in determining, for example, how changes in the parameters affect the frequency of an oscillator or whether or not it will lock onto a periodic stimulus. For the coupled PDE/ODE system in our population model, obtaining the adjoint solution requires solving a certain linear PDE along with some normalization requirements. We solve the linear PDE by backwards integration (see Material and methods) and to check whether we have found the correct solution, we need to look at an inner product that involves both the solution to the adjoint and the population oscillation. We confirmed the validity of the adjoint method developed in this study by the dual product of the zero eigensolution and its dual and comparison with the phase shift occurring when a perturbation is applied. Figure 3*a* shows a two-dimensional solution of the adjoint method, and figure 3*b* shows one-dimensional solutions of To check the numerical computation, we checked to see whether equations (2.25) and (2.29) were satisfied for each state on the limit cycle orbit. Figure 3*c* shows that the sum of the two terms corresponding to and is, indeed, constant for all time.

Equation (2.30) provides a formula for the PRF of the network when an identical current perturbation is applied to every cell. Thus, we can compare this bulk PRF to the phase shifts induced by a direct perturbation of the network with a short current injection. Therefore, another confirmation of the adjoint method was performed by comparison with PRFs by direct perturbation of the FPE (equations (2.6) and (2.7)) and the results are shown in figure 4*a*. We applied a square wave (amplitude 50, duration 0.01 ms) current pulse to the whole population at different points in the cycle and used this to compute the phase shift and the approximate PRF. Both PRFs determined by the adjoint method (with discretization of the PDE into 200 and 5000 bins) coincided with the results obtained by the direct perturbation of the FPE, which again indicated the appropriateness of the adjoint method. Similarly, the validity of was confirmed by the perturbation (amplitude 2, duration 0.01 ms) to the left hand side of equation (2.7) (figure 4*b*).

One of the main purposes for computing the PRF is to study the effects of coupling and periodic forcing of oscillatory systems. Coupling between inhibitory networks could be studied using this method, but, as most coupling between populations is mediated by excitatory outputs, we will not pursue the coupled problem further. However, periodic driving can also be analysed using the method and the experimental driving of inhibitory populations is possible [12]. Thus, to further test the utility of the adjoint method as a means of studying populations, we applied a weak sinusoidal current to the full population. We used current *I*(*t*) in equation (2.6) as
3.1where *I*_{0} = 2 and *I*_{app} and *ω*_{app} are the amplitude and frequency of the sinusoidal oscillation, respectively. The region of synchronization can be analytically estimated by the phase coupling function obtained by convolving the periodic current with the PRF [16,17]. The coupling function is defined as
3.2where *Z*_{macro}(*Θ*) is the macroscopic PRF equation (2.30). *Γ*(*Φ*) is shown in figure 4*c*. The maximum and the minimum value of *Γ* indicate tolerance of frequency mismatches between *ω*_{app} and the natural frequency of gamma oscillation without oscillatory inputs, *ω*_{nat}, for synchronization to a unit amplitude of oscillation (i.e. *I*_{app} = 1). That is, if the frequency difference is greater or less than the magnitude of *Γ*, then 1 : 1 locking is not possible.

Figure 4*d* shows the result of both the theory based on equation (3.2) and direct simulation of the FPE. We plot the difference between the network frequency and the forcing frequency as a function of the difference between the natural (with *I*_{app} = 0) and the forcing frequency. We can see that the synchronized region by numerical computations matches the analytical prediction, which is inside of the black thick lines.

### 3.4. Modulation of gamma oscillation by varying synaptic time constants

We used the macroscopic PRF obtained by adjoint method to predict the effect of changing the synaptic time constant on the frequency of gamma oscillation.

We denote *R*_{r}(*t*) as the sensitivity of the small change in *τ*_{r} to d^{2}*g*/d*t*^{2}. This sensitivity can be derived as the partial differentiation of d^{2}*g*_{0}/d*t*^{2} by *τ*_{r} and turns out to be
3.3

Therefore, the frequency changes of the macroscopic oscillation are evaluated as 3.4

Similarly, the influence of *τ*_{d} on the period of the macroscopic oscillations can be obtained by
3.5with
3.6Figure 5 shows the frequency of the macroscopic oscillation with different values of *τ*_{r} and *τ*_{d} obtained by numerical computation of the population of the modified theta model as well as those analytically predicted by the phase reduction. The slopes of the straight lines were obtained by equations (3.4) and (3.5) with the nominal values of parameters (*τ*_{r} = 0.5 and *τ*_{d} = 5). In both cases in which *τ*_{r} and *τ*_{d} were changed, the numerical computations matched the analytical predictions provided by the phase reduction across almost the entire region. The discrepancy at *τ*_{r} = 0.1 could be due to the fact that *τ*_{r} is getting close to zero where the model is singular. The plots nevertheless show that the population PRF can be used to estimate the frequency dependence on the characteristics of the synaptic inhibition.

To capture the influence of changing *τ*_{r} and *τ*_{d} on the macroscopic frequency in a wide parameter range of gamma oscillation, we evaluated Δ*ω*_{r} and Δ*ω*_{d} as
3.7by the adjoint method when changing *I* and *p*_{syn} with different values of *σ* (figure 6*a*,*b*,*d*,*e*). We confirmed that the signs of Δ*ω*_{r} and Δ*ω*_{d} were always negative over the entire regions shown in figure 6. In addition, Δ*ω*_{r}/Δ*ω*_{d} (figure 6*c*,*f*) varies between 0.85 and 1.87. It was larger than 1 over the almost the entire regions. However, it was smaller than 1 with *I* ≃ 1 and *σ* = 1. In general, Δ*ω*_{r}/Δ*ω*_{d} increased with increasing *I* and *σ*, whereas it was less affected by *p*_{syn}.

### 3.5. Effect of shunting inhibition

Experimental evidence has shown that the GABA synapses sometimes increase membrane potential for *V*_{R} ≤ *V*_{syn} [35]. This phenomenon is called shunting inhibition, and it was recently argued that it promotes macroscopic gamma oscillation by improving its robustness [35]. Here, we test this hypothesis using the FPE.

First, we numerically solved the FPE and evaluated the maximum firing probability, the average firing probability and *ω* as in figure 7*a*,*b* in the range of −70 ≤ *V*_{syn} ≤−55. Increasing the maximum value of the firing probability indicates that the amplitude of the gamma oscillation increases, which means that the robustness of the gamma oscillation is increased by shunting inhibition. The average firing probability and *ω* also increase as *V*_{syn} increases. These properties of shunting inhibition are the same as those reported previously [35].

Next, we evaluated the PRFs with shunting inhibition as shown in figure 7*c*. As *V*_{syn} increases, the amplitudes of PRF become smaller while maintaining the waveform, indicating that the phase of the oscillation is stable and robust against external interactions. (The larger is the PRF, the more sensitive is the oscillator.) Therefore, the results of PRF also support the finding that shunting inhibition improves the robustness of the gamma oscillations. The maximum value of PRF almost remains constant in the first part of the shunting inhibition in the range of *V*_{syn} <−62.5 mV and decreases in the range of *V*_{syn} >−62.5 mV (figure 7*d*), which differs from the other indexes (maximum, average of firing probability and *ω*) that are monotonic over the entire range of −70 ≤ *V*_{syn} ≤−55.

We computed the bifurcation diagram for the FPE using *V*_{syn} as a parameter and changing the other parameters to investigate whether the above effects of shunting inhibition hold in general. Surprisingly, when we set *p*_{syn} = 0.05, *I* = 2, the gamma oscillation disappeared above *V*_{syn} ≃ − 69.0 mV and a stable fixed firing state appeared (figure 8*a*–*d*). This result means that increasing *V*_{syn} can destroy the macroscopic gamma oscillation, which is in contrast to the study by Vida *et al*. [35]. More interestingly, when we increased the current to *I* = 3 with the same synaptic probability, there were two branches of HB for the macroscopic oscillations: below *V*_{syn} ≃ − 63.0 mV and above *V*_{syn} ≃ − 56.5 mV. Between these bifurcations, there was a region that had a stable fixed firing state (figure 8*e*–*h*), although the numerical computation fluctuated near the gamma frequency because of the vicinity of the bifurcation point and the finite-size effect. Therefore, the effect of shunting inhibitions is largely altered by the balance of the inhibitory synapses *p*_{syn} and current *I*. Finally, the stable region and oscillatory region are shown in the two-dimensional bifurcation diagram for *I* and *V*_{syn} (figure 8*i*).

## 4. Discussion

In this study, we first derived a modified theta model that possesses voltage-dependent dynamics and appropriate forms and strengths of the synaptic interactions. These properties are not incorporated into the conventional theta model which is the normal form for the saddle-node infinite cycle bifurcation [15]. Unlike related integrate-and-fire models with conductance-based synapses, the modified theta model is continuous with no resetting. By letting the number of neurons tend to infinity, we derived a hybrid PDE/ODE for the coupled neurons and their synaptic gates. The PDE for the population dynamics has simple periodic boundary conditions, and because it is continuous, requires no special methods for solving it. Furthermore, solutions to the discretized PDE can be numerically continued and bifurcations are easily detected using AUTO or other packages. Thus, we were able to find the parameter regions for macroscopic gamma oscillations and show that these arise via a supercritical HB. We found that the macroscopic oscillation tends to emerge with larger currents (*I*), higher rates of synaptic couplings (*p*_{syn}) and the weaker noise (*σ*). These results match the experimental findings that the pharmacological activations of interneurons by agonist of metabotropic glutamate receptors or kainate receptors induce large amplitude gamma oscillations and the loss of their interactions by GABA_{A} receptor antagonist bicuculline eliminates such inhibition-based gamma oscillations [8–10]. Thus, the analyses of the proposed model aid our understanding of the mechanism by which the inhibition-based gamma oscillation arises.

Populations of neurons are subject to external perturbations, so it is important to understand how perturbations affect the population rhythm. For limit cycle oscillations, the PRF describes how perturbations shift the timing of the rhythm. In the limit of short, weak perturbations, the shift in timing is described by the solution to the adjoint equation for the linearization about the limit cycle. Here, we derived the adjoint equation for the linearized Fokker–Planck and synaptic equations along with the limit cycle orbit. Accurate PRFs for the macroscopic oscillation are provided by solutions to the adjoint equation. As the adjoint method for solving a combined problem of ODE and PDE has not been described, we successfully developed a method by introducing the appropriate dual product. The validity of the derived PRFs was confirmed by the dual product of zero eigensolution and its dual (figure 3*c*), by direct perturbation of the FPE (figure 4*a*,*b*), by the region of synchronization to oscillatory inputs (figure 4*d*) and by the frequency shift caused by changing the synaptic time constants (*τ*_{r} and *τ*_{d}) to the population of the modified theta model (figure 5*a*,*b*). Although the PRF derived by the adjoint method is generally used to evaluate the effect of external perturbations [16,17], in our study, the macroscopic PRF was successfully applied to evaluate not only the effects of external perturbations but also how intrinsic parameter changes alter the macroscopic oscillatory dynamics. Recent experiments have used optogenetic tools to selectively drive subpopulations of neurons with periodic stimuli [12,21,36]; the adjoint method allows us to investigate the effects of such perturbations on the macroscopic gamma oscillations analytically.

Many studies have used the LIF model [25,35,37] or other neuronal models [26,38] to investigate the properties of population rhythms. Their analyses have often been confined to simulations of the population of individual neurons. This is partly because the FPEs for their models require complex boundary conditions, which makes mathematical and numerical treatments somewhat limited. In particular, it is hard to perform numerical bifurcation analysis on these PDEs. Although Brunel *et al*. derived the spike rate response from the FPE to evaluate the properties of the population of the LIF model [25,37], the resulting formulae are complicated and also valid only in the linear regime. By contrast, the population activity of the standard theta model is rather easily numerically analysed because the corresponding FPE can be treated with periodic boundary conditions and has no discontinuities [18,31,39]. The physiological origin of the gamma oscillations is not easy to discuss in the conventional theta model because the synaptic activity comes through additive currents rather than conductance-based synapses. By contrast, the modified theta model and its analytical framework developed in this study provide mathematically rigorous analyses for determining how collective dynamics depend on physiologically relevant parameters. Some theoretical work has been done on this with regard to inhibitory networks. Brunel *et al*. [25] showed that increases of the rise (*τ*_{r}) and decay (*τ*_{d}) times of GABA synapses decreased the frequency of gamma oscillations and that the sensitivity of *τ*_{d} was much lower than that of *τ*_{r}. Maex *et al*. also demonstrated a significant frequency decrease by increasing *τ*_{d} [40]. Using the adjoint method, we derived expressions for the sensitivity of the frequency to the synaptic time constants (figure 6). Our results, like those of previous authors, show that the frequency decreases with the increase in both the synaptic decay and the synaptic rise time. We also showed that the ratio of the sensitivities to the time constants, Δ*ω*_{r}/Δ*ω*_{d}, is smaller than 1; i.e. the changes in *τ*_{d} affect the oscillatory frequency more than *τ*_{r}, under the condition of weak current (*I* ≃1) and noise (*σ* = 1). This result differs from that of Brunel *et al*. [25], who focused on much higher frequency (more than 100 Hz) states of gamma oscillation where the slower scale of the decay time may not be as relevant.

In a pair of comprehensive papers, White *et al*. [41] and Chow *et al*. [42] systematically explored the effects of synaptic time constants on the firing period of individual neurons in inhibitory networks. In their analyses, the sensitivity of the firing period to the synaptic time constant is largely dependent on *I* and it is larger when the applied current is smaller, which is different from our finding that the macroscopic frequency is more sensitive to the time constants when *I* is large (figure 6*a*,*b*,*d*,*e*). Their result makes sense for individual neurons since at low drives, the neuron is close to the bifurcation from rest and, so, will be more sensitive to perturbations. Our situation is different since there are two ‘competing’ factors. To better understand why our results are different from the Chow *et al*. and White *et al*. work, the results of the adjoint method were compared for a higher current *I* = 4 and a lower one *I* = 1.5 with fixed parameters of *p*_{syn} = 0.2 and *σ* = 2. Figure 9*a* shows the PRFs and figure 9*b* shows the firing probability of the periodic orbit with *I* = 4 and *I* = 1.5. The macroscopic PRF is smaller with larger *I*, which would be related to the fact that the PRF of an individual neuron is smaller with larger *I* [43] (i.e. the individual neuron is farther to its bifurcation point). On the other hand, the large amplitude oscillation in firing probability is obtained with *I* = 4 (figure 9*b*), because it is far from the supercritical HB point. From equations (3.4) and (3.5), the sensitivity of frequency, *∂ω*/*∂τ*_{r} and *∂ω*/*∂τ*_{d}, is obtained by averaging of the time course of and respectively. Thus, we show these time courses in figure 9*c*,*d*. Regardless of smaller PRF of *I* = 4, they have larger negative values, which leads to high sensitivity. These larger values are caused because the large amplitude oscillation in firing probability influences these time courses through *R*_{r} in equation (3.3) and *R*_{d} in equation (3.6). Thus, the high sensitivity is obtained even though the macroscopic PRF is larger under the condition of smaller *I* (figure 9*c*,d). We can also see that this feature is apparent in the case of *τ*_{r} in figure 9*c*, which leads to large Δ*ω*_{r}/Δ*ω*_{d}.

From these sensitivity analyses, we see that the analyses of firing probability allow for a more quantitative evaluation of population activity, that would be difficult with only numerical computations of the finite-size neuron network [41,42]. The macroscopic PRF with its link to the microscopic one (equation (2.30)) enables us to understand the complicated relationships between microscopic structure and macroscopic oscillatory dynamics, and it cannot be achieved by the phase reduction of the simplified firing rate models of neuronal population [21,44].

Other relatively minor (hence small) types of interactions can also be taken into account as perturbations in this framework. In the electronic supplementary material, we evaluated the influence of a small number of gap junctions on the frequency of oscillation and obtained good agreement between theoretical prediction and numerical computations (see the electronic supplementary material, text S2 and figure S2 for details).

Shunting inhibition has been experimentally shown to promote macroscopic gamma oscillations [35]. Using weak coupling analysis and in the strongly oscillatory regime (all neurons fire on every cycle), the Gutkin group has shown that shunting inhibition can destabilize the perfectly synchronized state [11,30]. Our results show that the effects of shunting are complicated and depend on other parameters in the model. For example with lower bias currents, increasing the reversal potential of the GABA synapse destroys the gamma oscillation while at higher bias currents, there is only a small interval of reversal potentials where the gamma oscillations are destroyed. Thus, far from the bifurcation point (high drive), shunting inhibition can enhance gamma but near the bifurcation to gamma (low bias), the reduction in inhibitory drive (increased GABA reversal) pushes the population towards an asynchronous steady state. These non-trivial properties are observed not only by the modified theta model, but also by the conductance-based neuron model (see the electronic supplementary material, text S3 and figure S3 for details). Furthermore, it is also known that the reversal potential is higher than the resting membrane potential in the immature brain [45]; the systematic analyses of our model would be effective to understand the dynamics of the immature brain.

In the sparse firing state, it is difficult to systematically explore parameter space because the simulations must necessarily be long and involve many cells (in contrast to systems where every neuron fires in each cycle). For this reason, population models are a very useful tool as they result in PDEs that are deterministic. The proposed model enables us to evaluate the effect of synaptic interactions more rigorously than the conventional theta model. We demonstrated that the precise evaluation of synaptic interactions is important because the time constants and reversal potential greatly altered macroscopic properties (stability, frequency and PRF) of gamma oscillations as they are varied even within biologically plausible ranges. Thus, the proposed model and methods could be a key technique for elucidation of the micro–macro relationship of the gamma oscillation and its functional roles. Although we restricted ourselves only to networks of inhibitory neurons so as to limit the number of parameters studied, the same methods could be used to analyse networks of excitatory and inhibitory neurons. Such extensions would be very useful in the study of phase-locking of gamma across different areas as we could then use the derived adjoint functions (PRFs) to couple macroscopic population rhythms. Furthermore, we believe that the adjoint method for simultaneous equations composed of ODEs for interactions and a PDE for the probability distribution, that we have developed in this study, will provide a useful starting point for an in-depth understanding of various collective oscillations in biology that are composed of individual noisy oscillators with complex interactions including circadian oscillation in suprachiasmatic nucleus [2], quorum-sensing bacteria [3] and somite segmentations in vertebrate embryos [4].

## Funding statement

This study was supported by Grants-in-Aid for Scientific Research (23680024 to K.K. and 23240065 to Y.J.). K.K. was also supported by SCOPE (112103010). G.B.E. was supported by National Science Foundation grant no. DMS1219754.

## Footnotes

↵1 The FPE is numerically integrated by the Crank–Nicolson method. In the case of solving the adjoint equation, the Lax–Wendroff scheme is further adopted in order to improve the precision of the numerical integration. We basically divide the phase (−

*π*,*π*] into 200 bins except for figure 3 and blue lines in figure 4*a*,*b*, in which 5000 bins are adopted for accurate computation. In figure 4*a*,*b*, the results of 200 and 5000 bins are in good agreement, indicating that 200 bins are enough to calculate macroscopic phase response functions. For the bifurcation analyses by XPPAUT, centred difference for 100 bins is adopted and numerically integrated by CVODE [31]. This scheme preserves the total probability. Because of this conservation, there is always a zero eigenvalue. To eliminate this eigenvalue and thus fix the mass, we reduce the dimension to 99 equations and set, say,*P*_{1}d*x*= 1 − d*x*∑^{N}_{j=2}*P*_{j}. Here, d*x =*2*π*/100. With the elimination of the zero eigenvalue, we can use AUTO to compute steady states and their stability.

- Received January 18, 2014.
- Accepted February 27, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.