Distinct effects of heterogeneity and noise on gamma oscillation in a model of neuronal network with different reversal potential

Gamma oscillation is crucial in brain functions such as attentional selection, and is inextricably linked to both heterogeneity and noise (or so-called stochastic fluctuation) in neuronal networks. However, under coexistence of these factors, it has not been clarified how the synaptic reversal potential modulates the entraining of gamma oscillation. Here we show distinct effects of heterogeneity and noise in a population of modified theta neurons randomly coupled via GABAergic synapses. By introducing the Fokker-Planck equation and circular cumulants, we derive a set of two-cumulant macroscopic equations. In bifurcation analyses, we find a stabilizing effect of heterogeneity and a nontrivial effect of noise that results in promoting, diminishing, and shifting the oscillatory region, and is largely dependent on the reversal potential of GABAergic synapses. These findings are verified by numerical simulations of a finite-size neuronal network. Our results reveal that slight changes in reversal potential and magnitude of stochastic fluctuations can lead to immediate control of gamma oscillation, which would results in complex spatio-temporal dynamics for attentional selection and recognition.


Model
We start with a neuronal population described by the MT model with heterogeneity and noise. The MT model is a physiologically precise version of the theta model 23 . The phase of the i-th MT neuron satisfies the following differential equation: where C = 1(µF/cm 2 ) is the membrane capacitance, g L = 0.1(µS/cm 2 ) is the leak conductance, V T = −55(mV) is the firing threshold, V R = −62(mV) is the resting potential, and V syn is the reversal potential of GABAergic synaptic currents. The entire term in square bracket serves as coupling function. By changing the value of V syn , the sign of the term in square bracket will also change. If the whole term in square bracket is positive, then the effect of synapse is to increase the membrane potential of the i-th neuron, otherwise decrease the membrane potential. We note that the phase θ i is transformed from a version of quadratic integrate-and-fire neuron model (QIF) and the membrane potential of i-th neuron can be evaluated by V i = V T +V R 2 + V T −V R 2 tan θ i 2 (see "Methods" for details). ξ i (t) represents noise with �ξ i (t)� = 0 and �ξ i (t), ξ i (t ′ )� = δ ij δ(t − t ′ ) , and σ is the magnitude of noise. I(µA/cm 2 ) represents the input current. To employ two-cumulant truncation, following previous studies 32- 34 , we adopt a Cauchy-Lorentz distribution r(I) = 1 π � (I−η) 2 +� 2 , as the distribution of input currents, where η and are the center and width of the distribution, respectively. is the scale of heterogeneity, with larger for larger heterogeneity in the neuronal population. The i-th neuron fires when θ i exceeds π and modulates the membrane potential of the connected neuron by GABAergic synapses. g i syn represents the dynamics of GABAergic synaptic conductance. With mean-field approximation of the random and sparse connectivity, the dynamics of conductance obeys the following equation: where τ d = 5(ms) is the decay time constant, ḡ peak = 0.0214 (mS/cm 2 ) is peak conductance, N is the number of neurons in the neuronal population. In the numerical simulation of finite neurons, we set N = 3000 , which is considered to be an appropriate size for a typical layer within a single column 35 . P syn is the probability of random synaptic connections between neurons. A(t) is the firing rate of the neuronal population. Note that Eq. (2) is derived by mean-field approximation of the initial starting point (see "Methods" for details). Equations (1) and (2) constitute the microscopic model for numerical simulation as well as the starting point for deriving the macroscopic model.
For the derivation of a macroscopic model, firstly, we derive the Fokker-Planck equation (FPE) to describe the state of an infinite size neuronal network, and expand the probability density function (PDF) of FPE in a Fourier series. Next, we introduce the circular cumulant referred to in the novel dimension reduction method proposed and obtain the first two cumulants with the smallness assumption 32 . Since the σ term in Eq. (1) is multiplicative, which is additive in the reference paper, some modifications are required. (For a step-by-step derivation, see "Methods".) Then, the two-cumulant macroscopic model is derived as, , i is an imaginary unit, and * denotes a complex conjugate. Z and κ are the first and second order cumulants, respectively. Note that Eq. (3a) is in a form of the addition of four terms, that the first three terms are noise-free terms and the last term represents the effect of noise. If we set κ = 0 , the first three terms are exactly same as the dimension reduction result of the Ott-Antonsen Ansatz 36 . The firing rate A(t) of the macroscopic model can also be derived, Equations (2), (3) and (4) constitute the entire macroscopic model with two cumulants Z and κ . To begin with the analysis of this macroscopic model, we first describe the main analysis and validation idea of this paper. Our focus is whether the firing rate A(t) is stable or not. Although it is hard to analyze the probability density function in the Fokker-Planck equation Eq. (9) itself, we can analyze the equilibrium point and its stability of the two-cumulant macroscopic model derived in our study [Eqs. (2), (3) and (4)]. When the linearized equation around the equilibrium point has only negative eigenvalues in the real part, the firing rate A(t) is also stable in time. It undergoes a Hopf bifurcation when a pair of eigenvalues cross the imaginary axis due to certain parameter changes. This bifurcation finally results in the oscillation of the firing rate A(t). Since we consider five parameters in all in the macroscopic model ( η , V syn , P syn , , σ ), in order to scrutinize the effect of and σ , we plot 2-dimension bifurcation diagrams using two of the three variables ( η , V syn , P syn ) and study the changing curves due to the changing of and σ . In the validation part, we compare the bifurcation diagrams of the macroscopic model with the results of numerical simulations of finite neurons ( N = 3000 ) based on the microscopic model [Eqs. (1) and (2)]. All bifurcation diagrams of two-cumulant macroscopic model are plotted with XPPAUT 37 . All www.nature.com/scientificreports/ numerical simulations of the microscopic model are integrated by the Euler method, with a time step t = 0.01 (ms), using MATLAB (2019b, http:// www. mathw orks. com/ produ cts/ matlab/).

Results
In the primary research, we investigated the Hopf bifurcation with respect to the single variable η (Supplementary Material Fig. S1). Since η determines the center of distribution of input current I, increasing η could increase the average input current to neurons, thus activate the system. With different parameter settings, the system could achieve one of two states in the long run: stationary state or oscillatory state. We also showed the transition between two states by gradually increasing η with time in Supplementary Material Fig. S2, which reflects that the behavior of neuronal network turned from stationary state to oscillatory state in a certain value of the parameter.
Distinct roles of heterogeneity and noise in ( η , V syn ) plane. First, we analyze different roles of heterogeneity and noise in the ( η , V syn ) plane, as shown in Fig. 1A and 1F. The plane is divided into two regions by each curve: the region of the stationary state and that of the limit-cycle oscillation corresponding to the state with high gamma oscillation power. The curve itself serves as the boundary of two different states. Figure 1A shows the role of heterogeneity. When we keep the magnitude of noise σ = 0.1 fixed, with increasing heterogeneity , the curve moves rightward, leading to the enlargement of the stationary region. With any specific V syn , larger heterogeneity always tends to enlarge the stationary region of η . This enlargement effect is weakened with the increase of V syn , given that in the rightward movement of the curves in the lower half-plane is more obvious than in the upper half-plane. The range of and σ are set in 0.03-0.06 and 0-0.25 separately. In terms of , the changing of bifurcation diagrams is monotonic, so we only choose several values to show the tendency. In terms of σ , the range of σ is limited by two-cumulant truncation, because if σ is too large, cumulant higher than order two cannot be ignored. Next, we employ numerical simulation of finite neurons by [Eqs. (1) and (2)] to reveal the dynamics of single neurons and neuronal population. Positions marked with a plus sign and a multiplication sign are both on the right side of the red curve ( = 0.04 ), which is the oscillatory side, and on the left side of the green curve ( = 0.05 ), which is the stationary side. Figure 1B contains raster plots showing the behavior of all 3000 neurons, obtained by numerical simulation, under [ η = 1.59 , V syn = −56.5 , plus sign]. In order to eliminate the influence of the initial condition, the raster plots and time-courses are segments of simulation starting from 1000 (ms). In Fig. 1B, the raster plot with red dots corresponds to the behavior under = 0.04 at the position marked with a plus sign, and the one with green dots is the behavior under = 0.05 at the same position in the parameter plane. Note that the conditions of the raster plot in Fig. 1B match those in the same color in Fig. 1A. Such color matching is also applied to Fig. 1C-E and Fig. 1G-J. The raster plots show that when the parameter setting is on the oscillatory side of the bifurcation curve, neurons tend to fire synchronously. However, when the parameter setting is on the stationary side of the bifurcation curve, neurons fire asynchronously, which seems like random firing on the raster plot. Figure 1C is the time-courses of g syn by numerical simulation of finite neurons, under [ η = 1.59 , V syn = −56.5 , plus sign]. Figure 1C, D, it is clearly shown that the time-courses of the same position with different result in different states, very low amplitude oscillation (green curves) or gamma oscillation (red curves), which show excellent correspondence with the regional division in the bifurcation diagram in Fig. 1A. The raster plot Fig. 1D and the time-course Fig. 1E also agree with the regional division in Fig. 1A. The fluctuation close to stationary state (green curves) is due to the finite size effect because we assume the neuronal network has an infinite size in the derivation of the macroscopic two-cumulant model, while setting 3000 neurons in microscopic numerical simulation. The finite size effect could be regarded as extra noise of the mean field of order ∼ N −1/2 , where N is the size of finite ensemble 38 . In the stationary region close to the bifurcation curve, since there are always two conjugate eigenvalues close to imaginary axis, the extra drive by finite size effect yields small-amplitude resonance. In the oscillatory region, extra fluctuation by finite size effect makes the oscillation on each cycle slightly different with each other.
The effect of noise in the ( η , V syn ) plane is shown in Fig. 1F. We keep heterogeneity = 0.04 fixed, with the increase in the magnitude of noise σ . Unlike the simple moving effect of heterogeneity shown in Fig. 1A, the change in curves due to σ is more complicated and largely dependent on V syn . For small V syn in the lower half-plane, a larger magnitude of noise σ tends to shrink the stationary region of η , while for large V syn , a larger magnitudes of noise σ tends to enlarge the stationary region of η . For even larger V syn , an increasing σ seems to have no effect on the stationary region. Figure 1G contains raster plots of [ η = 1.6 , V syn = −57.3 , square], which clearly shows synchronous or asynchronous activity under different magnitude of noise σ . Comparing Fig. 1B with 1G, we can observe that two directions to stabilize oscillation (increasing or σ ) show similar raster plots, at least in the vicinity of bifurcation. Figure 1H is the time-courses of numerical simulation of finite neurons, under [ η = 1.6 , V syn = −57.3 , square] . The square is on the right side of the red curve ( σ = 0.1 ) and on the left side of the green curve, while the situation marked by the circle is exactly opposite. Figure 1G,H show that small noise gives rise to the gamma oscillation (red), and large noise gives rise to the stationary state (green). However, in Fig. 1I,J, small noise stabilizes the oscillation (red), while large noise arouses high amplitude oscillation (green). The raster plot and time-courses in Fig. 1I,J also agree with the result of the bifurcation diagram in Fig. 1F.
Besides, we investigated the effect of firing threshold ( V T ) in Supplementary Material Fig. S3. The bifurcation analysis shows that increasing firing threshold monotonically stabilizes the neuronal network. This result is in agreement with intuition since a higher firing threshold means harder to fire for single neurons. We also use numerical simulation of 3000 neurons to show that no matter the reversal potential higher or lower the firing threshold, the above result is correct.
In order to further validate the results of the bifurcation diagram obtained from the two-cumulant macroscopic model [Eqs. (2), (3) and (4) 1) and (2) (1) and (2)]. To obtain the total frequency power, we perform a fast Fourier transform to the time-course of [g syn ] and integrate the power spectral density over the frequency domain of gamma oscillation . To normalize the power difference between stationary points and oscillatory points and avoid negative value, we adjust the integration of power spectral density P γ by introducing a transformation P ′ γ = ln(10 4 P γ + 1) . The value used to plot heatmaps in Figs. 2 and 3 is the integration of P ′ γ over the frequency band of gamma oscillation . The power spectrum is showed in Supplementary Material Fig. S4.  Fig. 1F. Although there are some randomly scattered dots near the boundary of the two regions due to the effect of noise in the numerical simulation, these heatmaps and curves achieve an excellent agreement, indicating the correctness of the bifurcation diagrams shown in Fig. 1A,F. Besides, we found that the center frequency of oscillation monotonically increases with V syn in a wide frequency range, which includes both lower gamma and higher gamma range. (Supplementary Material Fig. S4).
Distinct roles of heterogeneity and noise in ( P syn , V syn ) plane. We next investigate the role of heterogeneity and noise in the ( P syn , V syn ) plane, and consider the rate of random synaptic connection P syn . Figure 3A represents the effect of heterogeneity . When we keep the magnitude of noise σ = 0.1 fixed, with the increase of heterogeneity , the regions inside the curves shrink, which enlarges the stationary state region, similar to Fig. 1A. With any given V syn , the oscillatory region of P syn monotonically shrinks. Figure 3B shows the role of noise in ( P syn , V syn ) plane. When we keep the heterogeneity = 0.02 fixed, shown as Fig. 3B, with increasing magnitude of noise σ , the curves twist and regions rotate. By further analyzing the nontrivial rotation, it is shown that the oscillatory regions of P syn are dependent on V syn . When increasing V syn , the oscillatory regions of P syn move rightward, while keeping the area of oscillatory region.
The validation of bifurcation analyses in the ( P syn , V syn ) plane is shown in Figs. 3C-E and F-H, which are plotted by the same method as Fig. 2. Figure 3C-E represent the effect of increasing heterogeneity with values: 0.015, 0.02, 0.025. Figure 3F-H represent the effect of increasing magnitude of noise σ with values: 0, 0.01, 0.02. Note that Fig. 3C,G are the same figures under [ = 0.02 , σ = 0.1 ], serving as the home position. Considering heterogeneity alone, the agreement of heatmaps and curves in Fig. 3A,C-E is excellent, showing the gradually shrinking oscillatory region with increasing . In terms of magnitude of noise σ , shown as Fig. 3B and F-H, one observes a clear rotation phenomenon of the oscillatory region in the heatmaps just as predicted by the bifurcation curves. This means that in some regions, gamma oscillation emerges with increasing noise, while in some other regions it stabilizes. Although there is some considerable difference between bifurcation diagrams and heatmaps due to finite size effect of numerical simulation, especially in Fig. 3H, the qualitative rotation effect can be clearly observed.
Besides ( η , V syn ) plane and ( P syn , V syn ) plane, we also investigated ( η , P syn ) in Supplementary Material www.nature.com/scientificreports/ models. Moreover, we confirmed that, similar to other bifurcation analysis ( Figs. 1 and 3), the effect of V syn is complex, which cannot be simply classified as stabilizing or facilitating oscillations.

Discussion
In this work, we formulated a set of macroscopic low-dimensional differential equations from an ensemble of spiking neuron models that enable us to analyze the entraining of gamma oscillation. The strength of our model is that it possesses voltage-dependent dynamics and analytically bridges micro-macro dynamics for gamma oscillations considering both noise and heterogeneity. Moreover, unlike numerical simulation of a large number of neurons by microscopic model, since our two-cumulant model is analytical, it explicitly shows the effect of each parameter and don't suffer from numerical issues such as computational cost. Though relatively complicated in the form, to some extent, our two-cumulant macroscopic model can mechanistically explain the reason behind phenomenon. In Supplementary Material Section 6, we analyzed how suppresses the oscillation by changing the eigenvalue of the Jacobian matrix of the system and found that it is similar to the effect of heterogeneity in Kuramoto model 1,2 . Instead, the rotation in Fig. 3B is really nontrivial and is found using our macroscopic equations in this study.
In research into synchronization of general coupled oscillators, both heterogeneity of the oscillator ensemble and induced noise have shown significant effects on synchronization. According to some previous research, heterogeneity suppresses synchronization 39 , while noise promotes the onset of synchronization in some cases 40,41 . In this work on neuronal populations, similar to previous research, heterogeneity diminishes macroscopic oscillation (Figs. 1A and 3A). The role of noise is largely dependent on synaptic reversal potential V syn , which is closely related to the coupling function in general coupled oscillators. With a different setting of V syn , noise promotes, www.nature.com/scientificreports/ suppresses or changes the region of synchronization (Figs. 1F and 3B). Therefore, our model demonstrates that even small changes in the coupling function can alter the effect of noise in collective dynamics. As the MT model incorporates voltage dependent membrane dynamics and synaptic interactions in a physiologically plausible manner, we obtain macroscopic Eq. (3) with somewhat harder computation than previous research. Ratas et al. 33 analyzed macroscopic dynamics of quadratic integrate-and-fire neurons. When considering both heterogeneity and noise together, they simply put forward the argument that effects of heterogeneity and noise are topologically similar in a bifurcation diagram. However, in our work, we found that the effect of heterogeneity and noise on gamma oscillation is different. Heterogeneity of a neuronal network stabilizes the oscillation. The reversal potential only monotonically affects the extent of this stabilization, but no qualitative effect. Nevertheless, the effect of noise is qualitatively dependent on reversal potential. Comparing with the work from Ratas et al., there are several differences that need to be mentioned. The first is we adopt a more detailed single neuron model with reversal potential V syn , which has been reported to contribute to the synchronous discharges for epilepsy in an in vitro experiment 13 . V syn determines whether the synapse increase or decrease membrane potential in our model.The second is they only consider excitatory neurons, while we varied the value of reversal potential within the experimentally observed range for GABAergic neurons. It basically decreases, but can increase in some cases, the membrane potential of the connected neuron. Thirdly, their single neuronal model is close to the firing threshold, which is not our case.
Gamma oscillation in a neuronal population plays important roles in brain functions such as attentional selection 3,42 , signal discrimination 43,44 and learning 45 . Therefore, controlling gamma oscillation is functionally required. Especially in the case of attentional selection, it is required for the higher visual cortex to shift the gamma synchronized target immediately to the attended area of the lower V1 population 42 . The reversal potential V syn , that is physiologically altered by intracellular calcium concentration through calcium-dependent Cl − transporters 28,29 , largely altered the dynamics of neuron.In our model, since the membrane potential is typically changing between the resting potential V R and the firing threshold V T , if V syn < V R , the synapse would always decrease the membrane potential. On the other hand, if V syn > V T , the synapse would always increase the membrane potential. This qualitative effect of V syn on synapse also shapes the bifurcation curve. Our results further reveal that V syn leads to immediate control of the gamma oscillation and synchronization. We also note that stochastic fluctuation in the membrane potential is increased by acetylcholine 46,47 , thus σ is also considered to vary more rapidly than other properties of the neuronal network. From these points, our results bridge gamma oscillation and ion channel, and imply that the dynamics of gamma oscillation can be well controlled by slight changes (in mV or sub-mV order) of V syn and σ , which have a potential influence on attentional selection and other cognitive functions. Because gamma oscillation exhibits complex spatio-temporal dynamics 3,4,6 and both reduced and excess oscillation are found in diseases such as autism 14 and schizophrenia 15 , our analyses shed light on another influence of the reversal potential and stochastic fluctuations on brain functions through gamma oscillation.

Methods
Transformation from quadratic integrate-and-fire model to modified theta model. The quadratic integrate-and-fire (QIF) model is a typical model for class I neurons near the firing threshold 48 . The dynamics of the membrane potential of the i-th neuron satisfy the following equation: where V T = −55 is the firing threshold, and V R = −62(mV) is the resting potential. The sign of V i − V syn determines whether the synapse increases or decreases the membrane potential: when negative, it increases the membrane potential, otherwise decreasing the membrane potential. The original form of g i syn (t) in Eqs. (1) and (5) obeys the following first-order equation: where ḡ peak = 0.0214(mS/cm 2 ) is obtained by dividing peak conductance of GABA on interneurons (6.2 nS) 52 by membrane surface area of neuron ( 2.9 × 10 −4 cm 2 ) 53 . δ(·) is the Dirac delta function and t (k) is the time of firing of pre-synaptic neurons connected to i-th neuron. With mean-field approximation, the sum of the Dirac delta function can be replaced, which transforms to Eq. (2).
In order to avoid the membrane potential value V i (t) jumping from +∞ to −∞ when firing, a phase variable θ i (t) can be introduced in Eq. (5), where −π < θ ≤ π . This variable transformation turns the QIF model into a modified theta model 23 . Compared with the infinite value of V i (t) , the phase value θ i (t) simply crosses the value of θ i (t) = π at these firing moments. Eq. (5) can be transformed into Eq. (1).
Derivation of two-cumulant model. Here, we elaborate on details of the derivation from the microscopic model [Eqs. (1) and (2)] of the two-cumulant macroscopic model [Eqs. (2), (3) and (4]. The dynamics of phase oscillator θ in Eq. (1) can be separated into a deterministic part u and a stochastic part w: www.nature.com/scientificreports/ where •dW is the (Stratonovich) white noise differential 49 . In the thermodynamic limit of an infinite oscillator ensemble ( N → ∞ ), its state can be described by ρ(θ, I, t) which satisfies the Fokker-Planck equation (FPE). Note that the multiplicative form of the stochastic parameter (with the Stratonovich interpretation) in our model Eq. (8c) is different from the additive noise term in the reference paper 32 . The corresponding FPE in our case is obtained as 50 It is helpful to separate the AC and DC parts of u at this stage as where c 3 = c 2 1 /(4C 2 ) , j ≥ 1 and α −1 = 0 . We note that the amplitude α j = π −π ρ(θ, I, t)e ijθ dθ is the Kuramoto-Daido order parameters 51 at a given I. Next, we integrate α j (I, t) over the distribution r(I) as Then, on the assumption that α j (I, t) is analytic in the upper half-plane of the complex variable I 36 , the integral can be computed by the residue on the upper half-plane as As a result, we substitute Ż j for α j in Eq. (12) and obtain an infinite series of ordinary differential equations for the order parameter Z j (t), where Z j series start from j = 1 , and set Z 0 = 1 , Z −1 = 0.
In order to obtain a set of low-dimensional reduced equations of the system Eq. (15), we follow a novel method 32 based on circular cumulants. Order parameters Z j = �e ijθ � can be regarded as the jth moment of the random variable e iθ , which are determined by a moment-generating function Thus, the related order parameters and its time derivatives can be written as, (11) ρ(θ, I, t) = r(I) 2π www.nature.com/scientificreports/ Substituting F(k, t) for Z j (t) in Eq. (15) using Eq. (17), and comparing the exponent of k on both sides of the equation, the partial differential equation for F(k, t) follows The circular cumulants introduced in 32 are determined by a cumulant-generating function defined as From Eqs. (16) and (19), one can derive the relationship between order parameters Z j (t) and circular cumulants χ j (t) . For the first two cumulants, Applying ∂ t to ψ in Eq. (19), we obtain the relationship of ∂F ∂t and ∂ψ ∂t as following Exerting Eq. (18) into Eq. (21), the partial differential equation for ψ(k, t) can be derived as where A, B, C, D contain the second, third, fourth, fifth partial derivatives of F(k, t) to t: The complete form of ψ is shown in Eq. (26). On the assumption that the smallness of the third cumulant is O(σ 4 ) , which vanishes in an approximation 32 , we only take the first two cumulants into consideration. By applying ψ(k, t) = χ 1 (t)k + χ 2 (t)k 2 to Eq. (22), we finally obtain the first two cumulants in the macroscopic model Eq. (3), where χ 1 = Z denotes the first cumulant and χ 2 = κ denotes the second cumulant. In order to achieve an explicit form of firing rate A(t) in Eq. (2), which is determined as one requires an explicit form of ρ(θ, I, t) first, coming from Eq. (11). With the approximation of only two cumulants, the moment-generating function F(k, t) in Eq. (16) can be simplified as F = exp[kZ + κ(k 2 /2)] . Following the assumption 32 that F ≈ [1 + κ(k 2 /2)] exp[kZ] under smallness of κ , the moments Z j can be derived by performing ∂ t to F as Z j = Z j + [j(j − 1)/2]κZ j−2 . Apply Z j to Eq. (11), and the summation of Fourier series is (20) Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.