Role of NMDA conductance in average firing rate shifts caused by external periodic forcing

A widely accepted view of computations in the brain relies on population coding, where the neural ensemble firing rate is modulated in a stable manner to transmit information and perform various cognitive tasks. At the same time, oscillatory neural activity is specifically modulated in frequency, coherence and power during cognitive performance. How the firing rate and oscillations interact remains a salient question. In this paper, we develop a theory for the interactions between oscillatory signals and the firing rate of neural populations based on activity of non-linear voltage-dependent NMDA synapses. Notably, we show under which conditions oscillatory inputs can control the mean firing rate without loss of stability. Using mathematical analysis and simulations of mean-field models, we demonstrate that presence of NMDA synapses on both the excitatory and the inhibitory neurons is critical for sinusoidal oscillations to significantly and stably increase the firing rate. We characterize the oscillation-induced mean firing rate shift as a function of the fast and slow synaptic weights and demonstrate the parameter region, in which the effect under investigation is mostly pronounced. Results of our work may help identify the properties of neural circuits that allow for constructive control of the firing rate codes by large-scale neural oscillations.


I. INTRODUCTION
Oscillatory activity is abundant across neural networks of the brain, and its profile is sensitive to various experimental and natural conditions (Wang, 2010). However, the relation between neural oscillations and the most basic functions of the brain circuits such as encoding and redistributing of information is poorly understood. One of the most prominent principles of information encoding in the brain is population rate coding (Rolls, 2011). According to this principle, information is stored as distribution of firing rates across a population. Firing rate in this case is defined as the average number of spikes emitted in a relatively prolonged time interval. In this context, a question of how might the firing rate modulations, necessary for information transmission and task execution, may interact with the oscillatory profile of the neuronal activity, comes to the forefront.
In many cases, oscillatory profile is closely related to population firing rate . Also, there is evidence that oscillatory activity is not merely a byproduct of rate coding, as it could be modulated independently of firing rates (Fries et al., 2001). In principle, oscillatory activity could be related to rate coding in several ways. First, it could serve as an independent dimension of neural code. Second, it could serve as a "metadata" that controls transmission of rate-coded information in the brain (Akam, Kullmann, 2014;Fries, 2005). Third, it could provide a modulatory scaffold for rate coding which stabilizes the code or makes the coding impossible by destabilizing the necessary firing rate dynamics (Dipoppa, Gutkin, 2013;Schmidt et al., 2018).
In order to algorithmically implement the aforementioned functions, there should be mechanisms of transformation between population rate code and oscillatory profile (Akam, Kullmann, 2010;Akam, Kullmann, 2014). The conversion from firing rate to oscillatory properties is quite well understood at the level of simple spiking networks (Brunel, Wang, 2003); also certain biophysical mechanisms were described for microcircuits (Roopun et al., 2008), as well as for single neurons. However, the opposite conversion, i.e. the effect of oscillatory properties on the mean firing rate, is much less studied.
To change the mean firing rate of a system by oscillatory entrainment, the system should be nonlinear. One candidate property is non-linearity of the f-I curves (gain functions) of neurons in the subthreshold state. The positive half-wave of input oscillations brings a neuron close to the threshold, thus increasing the firing probability, while the negative half-wave does not strongly affect this probability.
Thus, in average, oscillatory input modulation could increase the mean firing rate (Rolls, 2011;Salinas, Sejnowski, 2000); this effect has been explored mathematically in the limit of weak oscillations (Voronenko, Lindner, 2017). However, given that in very noisy state (typical for cortical networks) the gain functions are close to linear, the aforementioned mechanism should produce quite subtle effect.
Another possibility to link oscillations and firing rate modulation is to consider slow non-linear mechanism, such as NMDA conductance or voltage-dependent ion channels; this idea was proposed in the study (Akam, Kullmann, 2010), but has not been implemented in that study. In the present paper, we theoretically investigate the role that slow voltage-dependent NMDA currents could play in shifting the time-averaged state of an excitatory-inhibitory system in the presence of external zero-mean sinusoidal forcing. In order to separate the effect of NMDA non-linearity from the aforementioned effect of gain function non-linearity, we consider a low-dimensional population model with the gain functions linearized about the unforced equilibrium. Slopes of the gain functions, as well as the time constants that define the dynamics of population firing rates and voltages, are numerically derived from simulations of the uncoupled networks of leaky integrate-and-fire neurons. In our study, we mainly concentrate on excitatory effect produced by the external forcing. We demonstrate that in the presence of NMDA-receptors on the excitatory neurons only, the positive shift of the firing rate could only be very slight (1-2 Hz) without loss of stability, while the positive inhibitory firing rate shift could be more pronounced. However, when NMDA-receptor are also present on the inhibitory neurons, the system could be made more excitable without loss of stability, and the large positive shifts of both excitatory and inhibitory firing rates could be achieved.
The paper is organized as follows.
(1) We describe our model and its linearized version. (2) We provide analytical expressions for the relation between the parameters of the external forcing and the shift of the time-averaged state of the system that the forcing produces. (3) We describe the procedure of parameter selection that provides a realistic state of the system. (4) We perform geometrical analysis of the phase plane for the system with NMDA-receptors on the excitatory population only. We find the timeaveraged steady-states for the unforced and for the forced system. We also discuss the stability-related limitations for the steady-state shift that follows from our geometrical consideration. (5) We confirm predictions of our geometrical analysis by explicitly simulating the corresponding population model. (6) We perform the geometrical analysis and the confirming simulation for the system with NMDA-receptors on both populations (excitatory and inhibitory). We discuss how adding NMDA-receptors to the inhibitory population helps to overcome the limitations on the steady-state shift. (7) Finally, we compare the two aforementioned cases for a range of values of a system parameter that controls excitability of the fast subsystem.

A. The low-dimensional neural circuit system model
In order to describe the neural population dynamics of an excitatory-inhibitory neural circuit we use a low dimensional model. The two different neural populations are recurrently intercoupled with instantaneous excitatory AMPA and inhibitory GABAA, as well as slow non-linear excitatory NMDA connections. The network also receives external excitatory inputs represented by sum of white-noise signal with non-zero mean and sinusoidal signal with zero mean. We describe our system by the following six variables: (1) the excitatory and inhibitory population firing rates: , where the index a could be replaced by e for the excitatory population or by i for the inhibitory population (this notation will be used along the paper). Here NMDA ae J is the strength of the NMDA-related coupling between the excitatory population and the population a ; the function NMDA g describes the dependence of NMDA-current on the membrane voltage: where V is expressed in millivolts.
where , AMPA ae ai JJ are the strengths of the fast synaptic couplings (the first index defines the input population, and the second indexthe output population), ma g is the membrane conductance of neurons from the population a . The synaptic coupling strengths are expressed as follows: where ai j is the amplitude of the instantaneous inhibitory postsynaptic potential on the neurons from the population a ; AMPA ae j is the amplitude of the instantaneous excitatory AMPA-postsynaptic potential onto the neurons in the population a ; NMDA ae j -the amplitude of the NMDA-postsynaptic current step onto the neurons from the population a ; ab K -number of inputs that a neuron from population a receives from neurons that belong to population b (a,b = e,i), ma  -membrane time constant the neurons from the population a .
The time constant of NMDA-synaptic input is much larger than the time constants that govern the dynamics of the population firing rates and the membrane voltages, i.e. , , , . As seen below, we use the time scale separation to analyze the stability of the fast and the slow subsystems independently, as well as to analyze the effect of the forced oscillations of the fast subsystem on the dynamics of the slow subsystem.

B. Modelling assumptions
In the model described above, we made the following assumptions: (1) fast AMPA and GABAA synapses are instantaneous; (2) NMDA currents are used as dynamic variables, instead of NMDAconductances; (3) NMDA-currents are non-saturating; (4) voltage and firing rate dynamics are described by linear differential equations, without taking the effect of spike-to-spike synchronization into account; (5) The assumption (4) is valid when the system operates in a subthreshold regime, in which the activity of the corresponding spiking network is noise-driven and irregular (in this work, we check that the parameters we use satisfy these conditions). Note that this regime is different from a supra-threshold one, where the corresponding spiking network would show regular oscillatory spiking. We an argue that the assumption (6) is valid if the fast synaptic weights are sufficiently small, which is also true in our case. While assumptions (1) and (5) may be problematic in general, they only affect the relationship between the external forcing and the forced oscillations of the fast system. In our analysis, we start from selecting the parameters of the forced oscillations (as opposed to deriving them from the parameters of the external forcing), so the aforementioned relationships are not important in our case. The assumptions (2) and (3) are crucial, as they affect the non-linear dependence of the NMDA-currents on the firing rates and the membrane voltages. We use these assumption to simplify our analysis.
We should also note that we the expression (3)   .Yet this would make our analysis unwieldy while having presumably small effect on our main conclusions.

C. The linearized system
The goal of our study is to investigate the impact of zero-mean periodic forcing on shift of the equilibrium. In the system (1), this shift can be accounted for by non-linearity of the functions ,0 We now describe the linearization procedure. Let us assume that the full system has a stable fixed point denoted as Below, we will refer to (6) as the "unforced system".
where the appropriate derivatives are:

D. Fixed points of the unforced system
Let us now reduce the system (6) to two algebraic equations with the excitatory rate i.e., if the input a u moves the firing rate towards a r , then the same input moves the membrane voltage towards () aa Vr . From (7), it follows that: Next, let us find the state of the fast subsystem ( , , , ) where 0 ,0 (1 ) ( Now we can conclude that the steady-state firing rates and NMDA-currents should satisfy: Using (11) As we mentioned previously, the first and the second equations of (15) define two curves on the ( , ) e NMDAe rI -plane, intersections of which correspond to the fixed points of the system (6). Further in this text, we will refer to these curves as the e r -curve and the NMDAe I -curve, respectively. We should note that these curves are not nullclines, although they intersect at the fixed points of the system.

E. Analysis of the linearized system with external periodic input: forced oscillations
Our goal is to understand under what condition external periodic forcing may change the mean activity of the neuronal populations in our model . We start by analyzing the linearized system (6) with periodic external forcing. Here we derive the amplitude and phase relations between the external periodic signal and the forced oscillations of the excitatory and inhibitory populations.
Let us apply an external periodic forcing to the system (6), with the circular frequency  and the

complex-valued amplitudes of the oscillatory inputs to the E-and I-populations equal to
Using the time scale separation, we can assume that: (1) the slow variables The dynamics of ( , ) ei rr, defined by the system (16), could be expressed in the matrix form: NMDAe me e e e e re ee r Using (20), we can find the external amplitudes As the next step, we want to express the oscillatory part of the total input osc a u as a function of ( , ) Solving (23) with respect to A a V yields: In our analysis, we found it convenient to parametrize our system by The latter pair of parameters could be derived using (21). The relation (24) will be required in the next section, where we express the forcing-induced shift of the time-averaged equilibrium.

F. System with external periodic input: forced shift of the time-averaged equilibrium
In NMDAe NMDAi II-subsystem due to its non-linearity.
First, let us introduce short notations for the functions that govern the dynamics of the NMDAcurrents (see (6) and (2) where ** , aa rV are given by (11).
To account for the effect of the external forcing, we write the equations for the slow dynamics with , aa rV not equal to ** , aa rV (as in (26)) , but oscillating about ** , aa rV (see (17): In order to analytically estimate the integral in (28), let us expand From the expansion (30), we get: where we assume that , , , Now we put (31) into (29), omitting the higher-order terms, and calculate the integral in (28). As we integrate over the period of oscillations, all time-dependent cosine terms cancel out, and we get an autonomous equation: We can consider as a slowly moving "center", around which the fast subsystem oscillates under the external forcing. From now on, we will refer to the -system described by (32), together with the variables that functionally depend on , as the timeaveraged forced system.
The first term in the square brackets in (32) (containing AA ea rV) reflects the fact that NMDAcurrent depends both on the presynaptic firing rate and the postsynaptic voltage, while the second term (containing 2 A a V ) reflects the fact that NMDA-current depends non-linearly on the postsynaptic voltage.
Our numerical results (see the following sections) show that, for the selected parameters, the first term is much larger than the second one. Thus, the forcing-induced shift of the time-averaged equilibrium is mainly related to the joint effect of the presynaptic firing rate oscillations and the postsynaptic voltage oscillations, occurring with a small phase lag between them. We where r  is the angle between A e r and A i r .
Taking into account (25) and (26), the steady-states of (32) should satisfy: which is the analog of (13) in the case of the time-averaged forced system. Now let us exclude where i r and NMDAi I are the same function as it was introduced in (14). Similar to the unforced case, the steady-states of the time-averaged forced system could be found as intersections of the -curve and the -curve on the -plane, which are defined by the first and the second equations of (37), respectively. Both

G. Stability analysis of the unforced and the periodically forced systems
Let us now investigate the stability conditions for the systems with and without the periodic forcing. We want to identify conditions under which the periodic forcing changes the time-averaged activity without destabilizing it. For the unforced system, we can apply time scale separation and analyze stability of the fast and the slow subsystems independently.

A. Parameter selection
In this section, we briefly describe the selection of the primary parameters of our model, as well as calculation of the parameters derived from the primary ones.
In this work, we are mainly interested in the effect of the recurrent synaptic weights on the system's ability to change its time-averaged steady state under external forcing. To separate this effect from possible influence of the unforced steady state change caused by changes in the synaptic weights, we a priori select the steady state of the unforced system and then vary the synaptic weights, automatically re-tuning the external input strengths for each combination of the weights in such way that the pre-selected steady-state is kept constant. More precisely, we fixed (i.e. set a priori) the total mean inputs  The aforementioned procedures of parameter derivation are described in more details in the Appendix.
To pursue our analysis and make clear the role of the non-linear NMDA synaptic currents, we considered two models. The first model (Model 1) contained NMDA currents on the excitatory neurons only ( The primary parameters are listed in the Table 1; the derived parameters, as well as characteristics of equilibria are listed in the Table 2.

B. Parameter validation
In this section, we demonstrate that the parameters we selected are physiologically plausible and lead to realistic behavior of the system.
Post-synaptic potentials produced by activation of the fast receptors are given by the synaptic weights . We suggest that these parameters fall into physiologically reasonable range of values.
As we use linearized versions of the gain functions in our model, the firing rates could, in principle, become negative, which is non-sense. However, we fix the steady-state and the amplitude of the induced oscillations in such way that it does not happen. Specifically, the minimal firing rate for the excitatory / inhibitory population during the oscillations ( are larger than usually observed experimentally, however, they still suggest that the system is far from unrealistic oversynchronization and all-or-none spiking pattern driven by oscillations. In summary, the system demonstrates hallmarks of the sparsely synchronous regime characterized by moderate periodic firing rate modulation and irregular fluctuation-driven spiking of individual neurons which is believed to be typical for oscillatory cortical networks.

IV. ANALYSIS OF SYSTEM WITH NMDA SYNAPSES ON THE EXCITATORY POPULATION ONLY (MODEL 1)
In this chapter, we analyze the system with NMDA-receptors located on the excitatory neurons only ( 0 NNMDAi J  ). First, we plot the e r -and NMDAe I -curves and discuss their relation to existence and stability of the equilibrium in the time-averaged forced system. Next, we confirm the predictions of the phase-plane analysis by a direct simulation of the system. Finally, we discuss relations between the forcinginduced shift of the equilibrium and the parameters of the model.

A. Phase-plane analysis
From (44) and (45), it is seen that the e r -and e r -curves are identical and both represented by a straight line. In the Fig. 1(a), the phase plane is shown: the solid black line represents the S disappear, and the dynamics of the forced system diverges.

B. Numerical simulations of the low-dimensional system
To confirm the predictions of the geometrical analysis, we performed numerical simulation of the system described by (16). Simulation was performed during

C. Dependence of the steady-state shift on parameters
In this section, we analytically describe how certain parameters of the system affect the forcedinduced shift of the time-averaged equilibrium.
Let us rewrite the system (45) in the following form:   As we have seen, the requirement of stability of the time-averaged forced system strongly limits the possible shift of the excitatory firing rate produced by the external forcing. However, the shift of the inhibitory firing rate could be made sufficiently large. Using (11) Consequently, the inhibitory firing rate shift could be derived from the excitatory firing rate shift:

V. ANALYSIS OF SYSTEM WITH NMDA SYNAPSES ON THE EXCITATORY AND INHIBITORY POPULATIONS (MODEL 2)
In this section, we analyze the Model 2 that contains NMDA receptors both on the excitatory and inhibitory neurons. First, we note that the excitatory-to-inhibitory (E-I) NMDA coupling is much weaker in the model than the excitatory-to-excitatory (E-E) NMDA coupling , so the main effect of the external forcing on the system is still excitatory. The role of the E-I NMDA coupling here is to make the system more robust, which allows the fast subsystem to be more excitable without having divergent dynamics in the presence of the external forcing. In the Model 2, the fast subsystem was made more excitable compared to the Model 1 by decreasing the E-I AMPA coupling .

A. Phase-plane analysis
The phase plane for the Model 2 is presented in the Figure 3. The legend is the same as in the Figure  1A (which represents the same picture for the Model 1), with the difference that the S . Thus, we expect that no saddle-fold bifurcation (and, consequently, no divergent dynamics) will occur if we make the fast subsystem more excitable (as we observed in the Model 1). In the following section, we demonstrate that it is indeed the case, and that the system loses its stability via Hopf bifurcation, and starts to generate very slow oscillations.
We should note that stability of 0

B. Numerical simulations of the low-dimensional system
We confirmed the predictions of the geometrical analysis for the Model 2 by numerical simulation of the system (16), using the same simulation parameters as we did for the simulation of the Model 1. The results of simulation are represented in the Figure 4; the legend is the same as in the Figure 2.
One can see that the simulation results are in the good agreement with the predictions; however the prediction error is a bit larger than for the Model 1. This error is caused by the fact that the prediction is obtained using time scale separation which assumes the infinite ratio between the slow and the fast time scales, while in our system this ratio is, obviously, finite. To prove that the error is caused by the finite time scale ratio, we increased  Fig. 6(a), on the fast inhibitory synaptic weightsin the Fig. 6(b). The borders of instability of the unforced slow subsystem are denoted by the black lines; of the unforced fast subsystemby the green lines; of the time-averaged forced systemby the red lines. periodic forcing without destabilizing it. The main purpose of this analysis is to define how strongly the equilibrium could be shifted by the external forcing without loss of stability under various parameter combinations. In the Fig. 7(a), the two-dimensional (2-D) bifurcation diagram for the time-averaged forced system is presented, with on the horizontal axis and on the vertical axis. In the Fig. 7(b), a region of the diagram is zoomed in. The region to the left of the blue line in the Fig. 7(a,b) corresponds to instability of the unforced equilibrium 0 e r , and it is irrelevant for our analysis. The black line corresponds to saddlenode bifurcations: two fixed points merge and disappear as the line AB or BC is crossed from the right to the left, and as the line BC is crossed from the left to the right. In the region above the green line in the Fig.  7(a,b), one of the fixed points is a focus, so the system demonstrates either damped or sustained slow oscillations. This focus is stable to the right of the red line CG and unstable to the left of it; the red line CG itself corresponds to the supercritical Hopf bifurcations. The Intersection of the AB and CG lines is denoted by E.
To provide deeper understanding of the system's behavior, we fixed at various levels and analyzed how the fixed points of the system depend on . The values of we used are shown in the Fig. 7(a,b) by horizontal lines and marked by numbers. The corresponding one-dimensional (1-D) bifurcation diagrams are shown in the Fig. 7(c) by the blue / grey lines [and marked by the same numbers as in the Fig. 7(a,b)], with on the horizontal axis, and the steady-state values of on the vertical axis. A region of interest of the Figure 5(c) is shown in more detail in the Fig. 7(d). Stable parts of the diagrams in the Fig. 7(c,d) are represented by solid lines, and the unstable partsby dashed lines. Superimposed transparent thick lines denote parts of the diagrams that correspond either to stable or unstable foci. A diagram parts have grey color if they are located in the range of values for which the unforced system is unstable. All points of saddle-node bifurcations are connected by the black line, and all points of Hopf bifurcationsby the red line; the points A,B,C,D,G have the same meaning as in the Fig.  7(a,b). Please note that point E in the Fig. 7(a,b) is not a codimension 2 bifurcation, because the lines AB and CG correspond to bifurcations that occur on the different branches of 1-D diagrams. In the Fig. 7(d), the bifurcation points on the two stable branches for the AMPA ie J value that corresponds to the point E, are denoted by E1 and E2.
For all values below the cusp point B [see lines 1 -6 in the Fig. 7(a)], the corresponding 1-D diagrams contain a stable branch that terminates by a saddle-node bifurcation [see the lowest parts of the diagrams 1 -6 in the Fig. 7 In summary, we can distinguish two regions of values, based on the location of the maximal that could be achieved by varying . First, for a value below the point E in the Fig. 7(a,b), the maximal is achieved at the value of that corresponds to the saddle-node bifurcation that terminates the lower stable branch of the 1-D diagram at its left side [see diagrams 1 -5 in the Fig. 7(c,d)].
All such values are located on the curve AE1 in the Fig. 7(d). Second, for a value above the point E in the Fig. 7(a,b), the maximal is achieved at the value of that corresponds to the Hopf bifurcation [see diagrams 6,7 in the Fig. 7(c,d)]. All such values are located on the red line E2G in the Fig. 7(d). Thus, the dashed line AE1E2G in the Fig. 7(d) defines maximal values that could be, in principle, achieved by the external forcing without destabilizing the system. As one can see, the largest of these values is reached at the point E2; this value corresponds to the firing rate shift (relatively to the unforced equilibrium) that is equal to approximately 20 Hz.

VI. CONCLUSION
Numerous theoretical concepts of neural processing (e.g., selectivity to certain information) are formulated in terms of average firing rates. At the same time, brain activity demonstrates collective oscillatory patterns that correlate with functional states, task requirements and behavioral features. In order to build a theory that reconciles rate-based neural coding with functional role of oscillations in computations and information routing, one should consider a non-linear mechanism that converts oscillatory power into tonic firing rate shifts. In this paper, we explored a potential role that NMDA-synapses with non-linear behavior could play in shifting of the average level of population activity in the presence of external oscillatory input. We considered an excitatory-inhibitory population model with the population firing rates, mean membrane voltages, and NMDA-currents as the dynamical variables. NMDA-current depended both on presynaptic firing rate and postsynaptic voltage. In order to delineate the NMDA-related effects from the effects of non-linearity of spike generation mechanism, we linearized the dependences of the firing rate and the mean voltage on the synaptic input about the unperturbed equilibrium of the system. We applied the time scale separation and the time-averaging method, which allowed us to analytically derive expressions for the mean firing rate shift, induced by the oscillatory input to the system. Under realistic model parameters, we found that such shift is mostly produced by joint effect of the presynaptic firing rate and the postsynaptic voltage oscillations occurring almost in phase, and not by the non-linear dependence of NMDA-current on the postsynaptic voltage. Our analytical predictions were confirmed by direct numerical simulations.
We considered two models. In the Model 1, NMDA receptors were located exclusively on the excitatory neurons, while in the Model 2, they were located both on the excitatory and the inhibitory neurons. Using phase plane analysis based on our analytical results, we geometrically proved that the oscillation-induced firing rate shift is strongly limited by stability requirements in the Model 1. Adding even very weak NMDA input to the inhibitory population allows to overcome this limitation and achieve pronounced firing rate shift (up to 20 Hz) without destabilizing the system. Finally, we explored the parameter space and found optimal regions of parameters (NMDA input weights to both populations, AMPA input weight to the inhibitory population), in which the strongest firing rate could be achieved under the same amplitude of the entrained oscillations.
It is clear that many neural processes could potentially link oscillatory activity to mean firing rate modulations. As the examples, one could propose non-linearity of spike generation, non-linear behavior of various slow voltage-dependent ion channels, transition between spiking and bursting, or synaptic plasticity. Furthermore, the character of such link should critically depend on microconnectivity patterns. Nevertheless, the NMDA-based mechanism that we proposed here is rather general and should play certain role in oscillatory-induced firing rate shifts alongside with other potential mechanisms in most neural configurations. However the extent to which the NMDA-related non-linearity is involved in these shifts is a subject of future experimental research.

Linearization coefficients and the CV
In our numerical computations, we assumed that the source of the input variance is mainly external, so this variance does not depend on the state of the system. Following this assumption, we simulated a single excitatory (and a single inhibitory) leaky integrate-and-fire neuron that received Gaussian white noise with the amplitude 2 0 a  and the tonic input a u . We probed several values of a u around

Population time constants
In order to find population time constants ra  and Va  , we considered 5000 N  uncoupled neurons whose dynamics are governed by (A1) with