An E-I Firing Rate Model for Gamma Oscillations

Firing rate models for describing the mean-field activities of neuronal ensembles can be used effectively to study network function and dynamics, including synchronization and rhythmicity of excitatory-inhibitory populations. However, traditional Wilson-Cowan-like models, being amenable to mathematical analysis, are found unable to capture some dynamics such as Interneuronal Network Gamma oscillations (ING) although use of an explicit delay can help. We resolve this issue by introducing a mean-voltage variable (v) that considers the subthreshold integration of inputs and works as an effective delay in the negative feedback loop between firing rate (r) and synaptic gating of inhibition (s). Here we describe a three-variable r-s-v firing rate model for I-I networks, which is biophysically interpretable and capable of generating ING-like oscillations. Linear stability analysis, numerical branch-tracking and simulations show that the rate model captures many of the common features of spiking network models for ING. We also propose an alternative r-u-s model for ING which considers an implicit delay by a pre-synaptic variable u. Furthermore, we extend the framework to E-I networks. With our six-variable r-s-v model, we describe for the first time in firing rate models the transition from Pyramidal-Interneuronal Network Gamma (PING) to ING by increasing the external drive to the inhibitory population without adjusting synaptic weights. Having PING and ING available in a single network, without invoking synaptic blockers, is plausible and natural for explaining the emergence and transition of two different types of gamma oscillations. Author Summary Gamma rhythm (30-90 Hz) is a neuronal oscillation broadly observed across various animal species. It is associated with cognitive functions and gamma oscillation dysfunction has been found in neurological diseases. Traditional studies of gamma oscillations are primarily conducted in detailed but computational demanding spiking network models. Motivated by the growing size of accessible datasets and large-scale modeling, we develop mean-field representations of gamma oscillations that are computationally lightweight and mathematically analyzable. We resolve the incapability of generating oscillation in inhibitory network with traditional Wilson-Cowan firing rate model by introducing implicit delay with a mean membrane potential variable and generalize the model for excitatory-inhibitory network to capture richer dynamics. Our work provides a potential building block for multi-area large-scale modeling to help understanding the functional roles of circuits at different levels.

In the rhythmic PING state (Fig 1A), the I cells receive subthreshold external input and therefore 183 can only be recruited by the E cells. Their firing shuts down the E cells as well as their own 184 firing, hence their firing episode is brief. The recovery of E cells has two phases: a gradual phase 185 (during the first 15 ms after an I-cell episode) of slowly increasing E-cell firing from near silence 186 as synaptic inhibition wears off followed by a rapid phase (5 ms preceding the next I-cell 187 episode) where E cell firing accelerates because of recurrent excitation. This autocatalytic rapid 188 rise in E-cell activity recruits I cells and gives rise to another I-cell episode. The reciprocal 189 interaction between the E-I cells leads to the steady alternating PING-pattern of network activity: 190 increasing activity of E cells that recruits nearly synchronous firing of I cells. The two-phase 191 aspect of E-cell recovery can be understood by examining the dynamics of emergence and 192 termination of PING during a duration of I-cell drive that suppresses I cells (S1 Fig).   (Fig 1C), the cell-averaged time 215 courses are clearly rhythmic. Most prominent are the highly synchronized I-cell firing episodes 216 that generate a strong synchronized inhibitory synaptic output, gsyn,II and gsyn,EI with sharp onset 217 and exponential decay τsyn,I. These conductance transients lead to currents Isyn,II and Isyn,EI (dashed 218 blue, from I to E) that shut down I-cell and E-cell activity. Recovery of the E cells is set by the 219 slower decay time of gsyn,EI and by the gradual re-establishing of asynchronous E-cell firing and 220 a single inhibitory population always has one stable fixed point and therefore cannot generate an 266 oscillation (see Methods 2 for the proof; see also [34] immediately feeds back to suppress r. One resolution is to introduce an explicit delay for 277 activating s as in the r-s delay model as shown in [32]. However, the explicit delay presents 278 challenges for analytic treatment, so we propose involvement of an implicit delay. In the spiking 279 network, explicit delay is not required for generating ING (although a delay is incorporated in 280 several other models for ING, and found to render the behavior more robust, see e.g., [42]). The 281 subthreshold build-up of single neuron membrane potential contributes to the delay activation of 282 the synaptic gating and therefore is a candidate delay mechanism for a firing rate model. 283 Motivated by the spiking network simulation in Fig 1F, we introduced an additional dynamic 284 variable into the traditional Wilson-Cowan model that describes the temporal integration of mean 285 membrane potential v over cells. The model equations are given as follows. 286 where the input-output function f(v) is a sigmoid function 288 289 q(r) describes the dependence of the synaptic activation rate on the presynaptic firing rate by 290 In the r-s-v model, the additional v-variable can represent the subthreshold integration of inputs 292 and acts as an effective delay in the negative feedback loop between firing rate r and synaptic 293 gating of inhibition s. 294

295
The three-variable r-s-v model for an I population captures many dynamic properties of ING in 296 our simulations of an I-I network of , as developed in their systematic 297 study of conditions for coherent gamma oscillation in an I-I spiking network. We find striking 298 similarities in comparing the time courses of averaged firing rate, membrane potential and 299 synaptic gating variable from the spiking network simulation with the analogous variables r, v, s 300 respectively in firing rate model. In the isolated I-I network (Fig 2A top), the time courses for I 301 cells are similar with those in the spiking E-I network during ING (Fig 1F). The temporal 302 sequencing of v, then r and then s is successfully reproduced in the firing rate model (Fig 2B  303 top). The amplitude differences in the transients between the rate and spiking models are not 304 concerning since there is no uniform rule for normalization in these two models and the 305 amplitude of instantaneous firing rate of the spiking network depends sensitively on the window 306 size. Moreover, the level of synchrony of the spiking network depends on factors such as 307 heterogeneity in input drive and connectivity, while the firing rate model, by its nature, assumes 308 the heterogeneity without a quantification. However, one difference is worth noticing between 309 the models. In the spiking network model, the firing rate is pulsatile and drops to near zero 310 between spike clusters, matching with the fact that the I-cells are silent. This is not the case in the 311 r-s-v firing rate model. There, the firing rate is driven by f(v). After a firing event and before v 312 gets small enough such that f(v) ≈ 0, allowing r to decay to zero and flatten out, the inhibitory 313 activation s decays enough so that the subthreshold v starts increasing, able to activate the next 314 cycle. We also compare the time course of r-s-v model with the r-s delay model ( Fig 2C) and 315 observe the similarity except that r flattens a bit between I-firing events. With long explicit delay 316 (not shown), the r-s delay model can achieve the pulsatile firing rate as in the spiking network. We next compared the ING oscillation's dependence on the input II in these three models. Here, 348 we focus on the temporal maximum, minimum and mean of the variable r (Fig 2 middle)  when II is small (< 0.68), the network shows sparse asynchronous firing (Fig 2A middle). With 352 moderate II, the cells are above threshold and the recurrent inhibition restricts the firing to a brief 353 window when II + Isyn,II is at a minimum -the network is in ING oscillation. As II increases, the 354 intrinsic frequency of each cell also increases, leading to faster network rhythms (Fig 2A  355 bottom). When the input drive II is large enough (> 2.12), it dominates over the synaptic currents.  The ING oscillation frequency increases with input drive as for the spiking model (Fig 2B  362 bottom). The r-s delay model has similar characteristics as the r-s-v model in terms of the firing 363 rate (Fig 2C middle), but the network frequency has a nonmonotonic dependence on II (Fig 2C  364 bottom). To demonstrate the robustness of the r-s-v firing rate model, we explored the dependence of the 369 ING oscillation on the time constants (Fig 3) and other synaptic parameters such as the maximal 370 conductance and the reversal potential (Fig 4). Hopf bifurcation for oscillation offset at high τr is omitted for better resolution in region of 384 interest. The dependence on τr has similar features as for τs; τr can neither be too large or small 385 for oscillation. Notice that the oscillation amplitude shrinks rapidly as τr increases, indicating 386 weak modulation of firing rate (weak synchrony) unless τr is relatively small. The network 387 frequency decreases significantly as τr increases and plateaus beyond a few ms. The numerical observations can be validated by linear stability analysis of the three-variable 417 system (see Methods 3). We first showed that an exchange of the fixed point's stability can 418 change only by Hopf bifurcation. Therefore, to ensure the system has a stable limit cycle solution 419 for an ING oscillation, Hopf bifurcation must exist. Then by the trace condition for a Hopf 420 bifurcation, we derived a necessary condition for oscillation given by 421 The equality is attained when the effective time constant of the three variables are the same, i.e. 425 This confirms that none of the time constants can be zero or infinity, suggesting the existence of 427 two Hopf bifurcations with respect to each time constant. 428

429
We further examined how other synaptic dynamics related factors influence the existence 430 properties of ING in the rate model, in particular the effect of the synaptic activation time 431 constant τs /γI, the maximal inhibitory synaptic conductance gII and inhibitory reversal potential 432 + 2 on ING in the firing rate model (Fig 4). The II -range for oscillation shrinks a bit then stabilizes as γI increases. For a fixed II, the 442 frequency varies little with γI; the color bands are horizontal, as is the plot in the middle panel. the allowable range for γI increases but still has little impact on frequency (Fig 4A bottom). 461 462 Secondly, the maximal inhibitory synaptic conductance gII also needs to be moderate (Fig 4B  463 top). If recurrent inhibition is too weak, there is no oscillation, interpretable as failure to 464 synchronize in the I population, while too strong recurrent inhibition suppresses the activity of 465 the I population, i.e. very low firing rate. These effects match with the observations in spiking 466 network model [27]. The oscillation frequency has a non-monotonic dependence of gII, where the 467 minimal frequency is attained with a moderate synaptic conductance (Fig 4B middle). The loss 468 of oscillation at high gII can be compensated by increased II (Fig 4B bottom). As gII and II 469 increase together, the amplitudes of the input drive and the synaptic current both grow but 470 remain comparable, leading to an increase in network frequency. 471 472 Thirdly, hyperpolarizing inhibition rather than shunting inhibition benefits an ING oscillation. In 473 the rate model, the mean membrane voltage variable is normalized such that the rest potential 474 vrest is zero and the excitatory reversal potential satisfies , FFF = 1. Assuming, in unscaled variables, 475 Vrest = −65 mV and , FFF = 0 mV, the control inhibitory reversal potential is given by + 2 = −0.1 476 corresponds to + 2 = −71.5 mV. In the one-parameter bifurcation diagram of the inhibitory 477 reversal potential + 2 (Fig 4C top), the Hopf bifurcation with large + 2 value shows + 2 cannot be 478 too high above the rest potential, vrest = 0, for oscillation. The Hopf bifurcation with small + 2 479 value is far below a biological realistic value and not of interest. In the biologically realistic 480 regime (say + 2 > −0.5), the oscillation frequency increases monotonically with the inhibitory 481 reversal potential (Fig 4C middle). With fixed + 2 (Fig 4C bottom, vertical slices), the II range for 482 oscillation shrinks and disappears as + 2 gets to around 0.2. As II increases, the Hopf bifurcation The synaptic activation does not begin rising until nearly the peak-time of firing rate due to the 508 delay by u (compare the red vertical dashed lines in panels B and E). The firing rate stays near 0 509 for 5-10 ms before the next cycle begins. Indeed, the example time course of the r-u-s model (Fig 5A) shows activation of variables in the 520 order of r, u and s. One major difference between the r-u-s model and the r-s-v model is the 521 silent phase between cycles. In the r-u-s model, the synaptic activation s does not respond to the 522 rise of firing rate r directly, creating a time window for −ws + II to remain far below threshold 523 and therefore r stays near 0. From the time courses, s increases quickly as r crosses the threshold 524 in the r-s-v model ( Fig 5B) and s peaks before r falls to one-half its maximum value; while in the 525 r-u-s model, s begins to increase later, when r nearly reaches its peak, and s doesn't peak until r 526 is near to its lowest level. Moreover, in the r-u-s model, the oscillation amplitude is higher and 527 consequently the frequency is lower. This is because the firing rate r has longer time to increase 528 before the synaptic gating variable s is activated and begins to turn it off. The r-u-s model 529 depends on the input drive II similarly to the r-s-v model, but with a larger amplitude for firing 530 rate ( Fig 5C) and with a larger dynamic range (II ∈[1.05, 4.03]) that spans lower gamma 531 frequencies (f ∈[60.2, 92.4] Hz). A major parameter-dependent difference between the models is 532 on the steepness of the input-output function as we describe below. 533 534 We find that for generation of ING oscillations the input/output functions and synaptic activation 535 function of our rate models cannot be shallow. We describe the conditions about steepness, first, 536 with numerical calculation and then obtain insight from an analytic treatment. Using the 537 steepness (i.e., 1/(4κ), see Methods 3.3 for details) as a control parameter, we see a Hopf 538 bifurcation that sets the maximal κf ( Fig 6A) or κF ( Fig 6C) for oscillation, indicating that the 539 input-output function should be steep enough (κf or κF small enough). Notice that the allowable 540 range of steepness for the r-u-s model (near 0.15) is almost 6 times that for the r-s-v model (near 541 0.025), indicating a significantly relaxed steepness constraint for the r-u-s model. In both 542 models, the oscillation amplitude (peak firing frequency during an ING cycle) decreases as the 543 input-output function's steepness decreases (Fig 6A bottom for r-s-v model and Fig 6C bottom  544 for r-u-s model). The reason is that a steep f or F leads to a faster change (large dr/dt) in the 545 firing rate r and thus a larger oscillation amplitude. As a consequence, more time is needed for s 546 to decrease enough so that the synaptic inhibition is small enough and the next cycle of spiking 547 can occur. In terms of the change in the input drive II, the frequency of the r-u-s model is more 548 robust than that of the r-s-v model (Fig 6D versus Fig 6B). This is again due to the steepness 549 difference in the input-output function. In the r-s-v model, II drives the membrane voltage v. As 550 II increases, v passes through θf from below and r increases quickly. While in the r-u-s model, II 551 is part of the input argument to the less-steep function F, so the same change in II will lead to a 552 smaller change in the magnitude of r and thus a smaller change in frequency. If q(r) is linear, as often assumed in such formulations q(r) = r, then we see from s > 0 that 587 to be steep. The corresponding analysis applies similarly to r-u-s model (not shown). In the r-s 590 delay model, the steepness constraint is considerably relaxed. It can be shown both numerically 591 (S5 Fig) and analytically (S6 text) that as the explicit delay increases, less steepness in input-592 output function is required for oscillation.
where the E and I populations have different thresholds and degrees of steepness for the 602 functions f and q (see Methods 1.4). The six-variable system can be split into two subsystems, 603 where the E and I subsystems are described by the last and first three equations respectively. 604 Each subsystem acts on the other only by the conductance-driven synaptic coupling, i.e., sE 605 drives vI and sI drives vE. 606 607 Without fine tuning the parameters, we can reproduce with the firing rate model (Fig 7) the 608 transition from PING to ING by increasing the input drive II as seen in the spiking network (Fig  609   1). The example time courses for the E (top) and I (bottom) populations for PING ( Fig 7A) and 610 ING (Fig 7B) share similar features with the spiking network (Fig 1E and 1F). faster than that of PING. As we discussed for the spiking network, the oscillation period for 618 PING involves two phases -the decay of synaptic inhibition and the gradual activation of the E 619 population. That said, we find one notable difference between the rate and spiking models: the 620 activation of the E population in the firing rate model takes longer than in the spiking network  Hopf bifurcations corresponding to their own onset and offset (as in Fig 7C and 7D). For 692 somewhat larger IE, the PING and ING regimes are merged; at a critical value of IE, the middle 693 two Hopf bifurcations collide and disappear (Fig 7E and 7F)  We are motivated to develop rate models as complementary, or alternative depending on project 707 goals, to detailed spiking network models. We seek frameworks that retain biophysical 708 plausibility, mechanistic interpretability, mathematical tractability and that enable geometrical 709 and dynamical systems viewpoints. Rate models (also referred to as neural mass models) can be 710 computationally efficient and analyzable, providing opportunities to gain insights into local 711 network and multi-area brain dynamics [41,43,44]. In this study we develop and analyze rate 712 We also proposed a novel method to generate mean-field time courses and therefore effective 739 bifurcation diagram for spiking network models, which can be applied as a general way for 740 comparisons between spiking network models and firing rate models. shunting inhibition can not only support gamma oscillation but also improve the oscillation's 788 robustness. We tested this by adding a u-process to our r-s-v model so that both voltage-789 dependent synaptic current and presynaptic delay were represented. We found that, for a fixed 790 input drive, the maximal inhibitory reversal potential increased slightly and then decreased with 791 increasing τu (results not shown). We conclude that hyperpolarizing inhibition is still more 792 ING in a network of cells with Type II excitability. We are unable to make a comparison just 797 now for a rate model since we have yet to implement Type II cells in our rate models, but this is 798 an interesting future direction. 799 800

801
In three rate models, we found that the ING oscillation emerges and terminates via supercritical 802 Hopf bifurcations as the drive current II passes through critical values. Hence, the firing rate 803 oscillation amplitude increases gradually from zero with II at the onset of oscillations and 804 decreases to zero at the offset. We interpret this low amplitude mean-field oscillation as weakly- In an inhibitory spiking network with heterogeneous driving current, we found (but not shown 810 here) that I cells that receive strong drive typically fire more than those with lower drive. That is 811 to say, cells with small input drive skip cycles. We also observed in a network with 812 heterogeneous connectivity that cells with high in-degree skip more cycles than cells with low 813 in-degree tends to fire in each cycle (results not shown). spiking network and found the mechanism that generates higher frequency wins. We observed 840 qualitatively similar behavior in our six-variable E-I firing rate model (Fig 7) and in the spiking 841 E-I network of Hodgkin-Huxley-like units (Fig 1): PING transitions into ING as the drive to I 842 population increases. With strong drive to I cells, the ING frequency in the isolated I-I system is 843 high. The E cells fail to follow high frequency because either (1) The inhibition is too strong, so 844 the E cells are nearly silent, or (2) The network frequency is too fast, so the E cells spike 845 throughout the interval between two I-episodes. In contrast, when the drive to I cells is moderate, 846 the E cells are not overly suppressed and the time between two I-episodes is long enough, 847 allowing for both synaptic inhibition to decay and for E cells to recover. Transitions from PING 848 to faster ING in an E-I spiking network were also reported in [35] for increasing drive to I cells 849 and in [60] by increasing the synaptic conductance between I cells. 850

851
In our E-I rate model, the relative difference of PING and ING frequencies depend on the drive 852 to both E and I populations. As illustrated in Fig 7A and 7D, PING is significantly lower in 853 frequency than ING and also compared to PING in the spiking models. While the rate model's 854 PING frequency is insensitive to the drive in the I population, it depends strongly on the drive to 855 the E population. For example, PING frequency doubles by an increase of 0.2 in IE ( Fig 7E) and 856 with little change in the ING frequency. For this study we did not fine-tune parameter values 857 (such as time constants) of the E-I firing rate models to achieve strict quantitative agreements 858 with data or spiking models, although this could be done in future work. Steep gain function in rate model 899 We have found that some gain function in the rate models needs to have a steep slope. This 900 feature needs closer investigation. This constraint may be linked to the lack of a regenerative 901 process. A formal derivation of these rate models from, and comparison with, spiking models 902 would likely shed some light here. In rate models, the gain function is typically related to the 903 cellular frequency-current relation which can have a steep rise near the onset, smoothed a bit by 904 the assumed effects of noise and heterogeneity. However, we note that the operating range of 905 ING in our spiking models is not necessarily in this steeply portion of the frequency-current 906 relation. Interestingly, the three-variable rate model for a network of heterogeneous QIF units 907 We encourage efforts for further analysis and derivations that may reveal how intrinsic 934 properties together with synaptic dynamics, gsyn(t), and interactions between units can that lead 935 to rate models that provide meaningful and easily interpretable descriptions. We can hope for 936 rate models that are applicable across the range from weakly modulated sparse firing to more 937 synchronized episodes and across a range of frequencies, that can describe the situation of E-I as 938 well as allowing for potential dominance of I-I. Can we bridge from stochastic dynamics, say, 939 based on the integrate-and-fire family as starting points that lead to low-order MF models and 940 perhaps parameter-tunable based on electrophysiological data? Maybe the WC-framework can 941 be tweaked and extended to overcome some of the current limitations. What form of gain 942 function, formally derivable, is more appropriate than an assumed high-gain sigmoidal-like 943 input-output relation? By developing procedures or guidelines to fit rate-model parameters to 944 simulations of spiking models, one could, hopefully, illustrate how to approximately fit the firing 945 rate model to experimentally recorded data such as cell excitabilities, synaptic input kinetics and 946 connectivity. 947 948 Looking forward 949 We hope that our results for these idealized extensions of the WC-framework provide options for 950 further applications and extensions of rate models to describe qualitative aspects of E-I local 951 network dynamics. The role of inhibition in balancing and stabilizing excitation is being well-952 explored with rate models (e.g., [64]). Destabilization of the inhibition-stabilized state and 953 emergence of network oscillations is a natural application (e.g., [65]). In consideration of co-954 existent slower dynamical mechanisms, such destabilization could describe with rate models, 955 say, gamma oscillations nested within a theta rhythm (e.g., [66,67] The input-output function takes the same form with the r-u-s model except for different 993 parameter values of θF and κF. 994 The control parameters that generate the time course in Fig 2C are listed in Table 3. 996 997  1000 We introduced an alternative way to include an implicit delay in Fig 5. The model's equations 1001 are 1002 The input-output function is 1004 As for the r-s-v model, the scaling factor of synaptic activation rate is 1006 The control parameters that generate the time course in Fig 5B are listed in Table 2. The values 1008 were chosen to be the same as those for the r-s-v model when the corresponding parameters are 1009 shared in the two models. All parameters were fixed except for the inverse steepness κF and the 1010 external drive II in the bifurcation diagrams.
The E and I populations have different parameters for the input-output function, where the 1020 function for I population is steeper than E population.

rE and rI -steady state relationships (SSRs) for the six-variable r-s-v model
1033 To obtain the rE and rI -SSRs, we first set each equation in the system into steady state. Then we 1034 write vE as a function of rE using drE/dt = 0 and sE as a function of rE using dsE/dt = 0. The 1035 relations are given by 1036 1 = 1 '( ( 1 ) 1037 1 = γ 1 1 ( 1 ) + $1 γ 1 1 ( 1 ) + 1 1038 Similarly, we write sI, vI as functions of rI. 1039 Graphically, the λ-values can be interpreted as the roots of a vertically shifted cubic curve, since 1096 the first term is a cubic term with 3 negative real roots and the second term is positive. By  1097 plotting the left-hand side as a function of λ and checking the number of intersections with the 1098 zero-line, it can be seen that the shifted-up cubic curve has 3 negative roots if the shift-up term is 1099 small; it has 1 negative root and a pair of complex roots if the shift-up term is large. Therefore, 1100 the only way for the fixed point of the system to exchange stability is via Hopf bifurcation. 1101 1102

Dependence of the existence of Hopf bifurcation on parameters 1135
In order to have a stable periodic solution, the existence of a Hopf bifurcation is necessary. 1136 Therefore, the minimum of the right-hand side of Eq 2 given by Cauchy-Schwarz inequality 1137 needs to be smaller than the left-hand side. That is, the three effective time constants being close 1138 to each other is a condition that favors the existence of oscillation. Suppose one of the time 1139 constants τr, τs, τv is very small (goes to 0) or very large (goes to ∞), then the right-hand side of 1140 Eq 2 will go to ∞ and there can be no solution for Hopf bifurcation. where q(V) = 1 for 0.5 ms after V crosses 0 mV from below and q(V) = 0 otherwise. The rising 1193 rates are αE = 14/3 ms −1 , αI = 53/11 ms −1 , and the decay rates are βE = 1/2 ms −1 , βI = 2/11 ms −1 . 1194