Incorporating slow NMDA-type receptors with nonlinear voltage-dependent magnesium block in a next generation neural mass model: derivation and dynamics

We derive a next generation neural mass model of a population of quadratic-integrate-and-fire neurons, with slow adaptation, and conductance-based AMPAR, GABAR and nonlinear NMDAR synapses. We show that the Lorentzian ansatz assumption can be satisfied by introducing a piece-wise polynomial approximation of the nonlinear voltage-dependent magnesium block of NMDAR current. We study the dynamics of the resulting system for two example cases of excitatory cortical neurons and inhibitory striatal neurons. Bifurcation diagrams are presented comparing the different dynamical regimes as compared to the case of linear NMDAR currents, along with sample comparison simulation time series demonstrating different possible oscillatory solutions. The omission of the nonlinearity of NMDAR currents results in a shift in the range (and possible disappearance) of the constant high firing rate regime, along with a modulation in the amplitude and frequency power spectrum of oscillations. Moreover, nonlinear NMDAR action is seen to be state-dependent and can have opposite effects depending on the type of neurons involved and the level of input firing rate received. The presented model can serve as a computationally efficient building block in whole brain network models for investigating the differential modulation of different types of synapses under neuromodulatory influence or receptor specific malfunction.


Introduction
"N-methyl-D-aspartate (NMDA) and dopamine (DA) receptors and their interactions control an incredible variety of functions in the intact brain and, when abnormal, these interactions underlie and contribute to numerous disease states.These receptor interactions are relevant in such diverse functions as motor control, cognition and memory, neurodegenerative disorders, schizophrenia, and addiction."(VanDongen, 2008) Our work has branched out of an effort to develop a biophysical mechanistic framework for incorporating dopaminergic neuromodulation in whole brain network models (Sanz Leon et al., 2013), particularly with the aim to create personalized models for Parkinson's disease (PD) patients, in the spirit of what has been done for the case of epilepsy (Jirsa et al., 2017).PD is associated with loss of dopamine in the dorsal Striatum due to degeneration of dopaminergic neurons of the Substantia Nigra Pars Compacta.The resulting dopamine depletion in the Basal ganglia (BG) is accompanied by changes in firing rates, as well as patterns of firing and levels of synchronization of neurons in the BG subcortical circuit (Bergman, 2021).However, the dynamic alterations are not merely a local phenomenon but are instead distributed across the macroscopic BG-thalamocortical networks (Tinkhauser et al., 2018); the increase in beta burst prevalence is accompanied by greater concurrence in time of beta burst activity across the network, during which highly stable phase locking is observed between the activity of cortical regions and BG nuclei (Yu et al., 2021).This motivates us to develop needed mathematical formulation for incorporating dopaminergic Action Editor: Jonathan D. Victor neuromodulatory processes in the Virtual Brain Modeling platform (Sanz Leon et al., 2013), in which parcellated brain areas are represented as a set of nodes on a network, each endowed by a dynamical system describing the aggregate neuronal activity in a given brain region.A first step towards this aim is to devise a way to capture the effects of dopamine on a mesoscopic (neuronal population) level.
We take as a starting point the seminal work in (Humphries et al., 2009) which argued that dopamine action on medium spiny neurons of the striatum can be captured in spiking neuron models as a scaling in the maximal conductances of NMDA, AMPA and GABA receptors.A similar framework was also suggested in (Durstewitz et al., 2000) for modeling dopamine action on pyramidal cells of the prefrontal cortex.More specifically, in the striatum, it is reported that dopamine enhances NMDA receptor (NMDAR) currents in neurons with D1-type dopamine receptors and attenuates AMPA receptor (AMPAR) currents in neurons with D2-type dopamine receptors, and that the balance between the two is at the core of the proper functioning of the striatum (Humphries et al., 2009;Lindahl & Kotaleski, 2016).Dopamine is also reported to enhance NMDAR currents via D1-receptors in the prefrontal cortex, while attenuating AMPAR currents and simultaneously enhancing GABA receptor (GABAR) currents (Durstewitz et al., 2000).As such, the effect of dopamine is, in fact, statedependent and can switch from a predominantly net inhibitory effect (in low-activity states) to a net excitatory effect (in high-activity states) (Durstewitz et al., 2000).This special feature of dopamine action is facilitated by a unique property of NMDAR, which is its voltage-dependent nonlinear magnesium (Mg 2+ ) block, a different mechanism than that governing the voltage-gated channels that generate the action potential.At the resting membrane potential, extracellular Mg 2+ binds tightly to a site in the pore of the channel, preventing any ionic current flow, even in the presence of glutamate.It is only when the membrane is already depolarized enough (such as due to the opening of AMPAR channels) that the Mg 2+ block is removed and ionic currents can flow (Jahr & Stevens, 1990;Nowak et al., 1984).As such, NMDARs are often referred to as "coincidence detectors", allowing maximal current flow when two conditions are met: glutamate is present and the cell is already depolarized.
This Mg 2+ block-induced nonlinear voltage-dependence of NMDAR conductivity, which is at the core of the role that NMDAR plays in neuronal dynamics, is often either omitted, for mathematical convenience, from traditional neural mass models (Wilson & Cowan, 1972) or reduced to a linearization around a local working point of mean membrane voltage (Brunel et al., 2001).More recently, a novel approach was put forward to derive exact mean field models of populations of quadratic-integrate-and-fire neurons by invoking the Lorentzian ansatz for describing the distribution of membrane voltage values (Montbrió et al., 2015).The resulting models, coined "next generation neural mass models", offer the advantage of retaining the mean membrane voltage of the population as an explicit dependent variable, along side the mean firing rate, in the final mean field equations (Coombes, 2023).In the work we present here, we show that this allows us to incorporate in the derivation the full nonlinearity of the NMDAR conductance.This is done by rewriting the expression for the nonlinear NMDAR conductance as a piece-wise quadratic polynomial to maintain the applicability of the Lorentzian ansatz.This consequently permits the analysis of the effect of NMDAR across the full range of dynamical regimes of the system.
We consider a population of what is often referred to as an Izhikevich spiking neuron, which was derived as a canonical model for spiking neurons and, for different values of its four parameters, can reproduce different spiking and bursting behaviour of known types of neurons (Izhikevich, 2018).The Izhikevich neuron is basically a quadratic-integrate-and-fire (QIF) neuron with an added membrane recovery (adaptation) variable that provides negative feedback (Izhikevich, 2003).While the next generation neural mass models that exploited the Lorentzian ansatz where first derived for populations of simple QIF neurons (Byrne et al., 2017;Coombes & Byrne, 2019;Montbrió et al., 2015), the same approach was soon thereafter shown to be applicable for the case of Izhikevich neurons (Chen & Campbell, 2022) and the closely related QIF neurons with slow adaptation (Ferrara et al., 2023).We will here build on these efforts and extend this approach to include NMDAR type synaptic currents with nonlinear voltage-dependent conductance in the derivation of a mean-field model for a population of Izhikevich type neurons.
In the next section, we present the equations governing the population of neurons and the main steps involved in the derivation of the corresponding mean-field model.We then present the results of the numerical analysis of the dynamics for two example cases, that is, with parameter sets corresponding to two different types of neurons: 1) a population of excitatory regular spiking (cortical) neurons (Izhikevich, 2003) and 2) a population of inhibitory striatal neurons (Humphries et al., 2009).We generate bifurcation diagrams that show transitions in stability of constant firing rate states, as parameters vary, and emergence of oscillatory limit cycle solutions.We also present sample time series of simulations displaying the different types of possible oscillatory solutions.All results are contrasted to the case of omitting the nonlinear Mg 2+ block from the NMDAR term, that is, the case of assuming linear NMDAR conductance-type current with the same slow synaptic decay time constant.The analysis shows that the NMDAR nonlinearity can lead to a voltage-dependent shift in the bifurcation diagram and a change in the relative size of the different operating regimes, with possible disappearance of the high constant firing rate states.In conjunction, in the oscillatory regime, the NMDAR nonlinearity can lead to alterations in oscillation amplitudes, frequency power profile, as well as the level of neuronal synchrony within the population.Ultimately, the presented mean-field model exhibits a rich range of dynamical behaviors and can serve as a building block for computationally efficient mechanistic whole brain models the explicitly incorporate the differential effects of neuromodulatory forces as a heterogeneous modulation of AMPAR, NMDAR and GABAR synaptic currents.

Derivation of mean-field equations
We consider a population of all-to-all coupled Izhikevich neurons with conductance type synapses.The single neuron equations are of the form (Izhikevich, 2018): where V is the membrane potential, U is the recovery cur- rent, C m is the membrane capacitance, V r is the resting membrane potential, V t is a threshold potential, and k are scaling factors, u is a time constant for the recovery variable, and I is an external applied current.The equations are complemented with the following spike reset condition: if whenever the membrane potential grows beyond the set peak value, V peak , it is reset back to the value, V reset , and the recovery current is augmented by a constant value U jump .
The synaptic currents can include excitatory AMPAR and NMDAR currents as well as inhibitory GABAR currents, such that: (Jahr & Stevens, 1990).
) correspond to the reversal potential and ampli- tude of conductance for GABAR, AMPAR and NMDAR, respectively.Ignoring the rise time for synaptic channel opening, the conductance amplitudes are governed by the following equations: with subscript s ∈ {A, N, G} , referring to AMPAR, NMDAR and GABAR, respectively; syn,s is the synaptic decay time constant, G jump,s is the increment in conductance per spike received by a neuron, and N is the total number of neurons in the population.After non-dimensionalization, the equations take the following simpler form: We have here included the i term, as a background current that introduces a heterogeneity among the neurons of the population.The scaled variables and parameters are related to the dimensional ones as follows: The expression for the Mg 2+ voltage dependent nonlinear NMDAR current was, historically, obtained empirically by numerical fitting of experimental data (Jahr & Stevens, 1990;Nowak et al., 1984), so we can take a step back and instead rewrite it as a piece-wise polynomial.It can be seen in Fig. 1 that f N (v) can be well approximated as: with v cut = −1.2; the range of values for v was chosen to be wide enough to cover the range of V presented in the origi- nal empirical fit (Nowak et al., 1984).Given the reversal potential for NMDA is E N = 0, e N = 1 , then the coefficients of the polynomial fit will only vary with the choice of V r .For the fit in Fig. 1, we used V r = −82.66mV (as in the first example case presented in Section 3).
The polynomial fit, as well as the discrete interval boundaries, were obtained by minimizing square errors while ensuring smoothness at the boundaries (that is, ensuring continuity in f N (v) and f N v ); the fitting was done using the "scipy.optimize.leastsq"python function and the goodness of fit was not found to be sensitive to the choice of V r value.
Given that the time scale of the u dynamics is much slower than that of v , we can invoke the adiabatic approximation and consider that all neurons experience a common mean adaptation current, u = ⟨u i ⟩ .We follow the same rigorous derivation presented in (Chen & Campbell, 2022), so for brevity, we only provide the main steps here and refer the reader to the latter reference for more elaborate details.
In the limit of N → ∞ , the distribution of the membrane potential is governed by the following continuity equation: We assume a Lorentzian distribution for each of the membrane voltage and the heterogeneous background current: x( , t) and Δ correspond to the widths of the respective dis- tribution, while y( , t) and correspond to the respective centers.
The different terms of the continuity equation can be expanded as follows: Plugging the above expressions in the continuity equation and balancing coefficients of v 2 & v , respectively, we get: We then define w = x + iy , such that −iw 2 = i y 2 − x 2 + 2xy and ẇ = ̇x + i ̇y , then the above ̇x and ̇y equations can be com- bined into: The firing rate can be obtained as the flux probability at v peak , as v peak → ∞ ; for a given : Here, we have used the fact that as Then the mean firing rate is: For a given value of , the variable y can be computed as the Cauchy principal value of ∫ (v| , tc)vdv .Then, the mean membrane potential can be expressed as: The residue theorem can be applied in a closed contour integral in the lower half complex plane containing the pole of g( ), * = − iΔ , then: Fig. 1 Left: ( ⋅⋅ red) piece-wise approximation of (-blue) f(v).Right: ( ⋅⋅ red) smooth approximation of (-blue) p 2 (v); for parameter values for the example case of excitatory regular spiking neurons reported in the next section Here, we have neglected the ġ N contribution since the NMDAR dynamics is an order of magnitude slower than the v dynamics.
Then, evaluating the equation for ẇ at = − iΔ , we arrive at the equations governing the mean firing rate and mean membrane potential: the final mean field equations can more compactly be written as: In addition, the equations governing the mean adaptation current and mean synaptic conductances take the following form: Here, the r in the g s equation is the mean firing rate of the population in the case of recurrent synapses, or the mean firing rate received by the population in the case of an external input signal.
The coefficient p 2 appearing on the right-hand side of the ̇v equation is defined piecewise as in Eq. ( 2), but to avoid dis- continuity in the mean-field equations, we replace it with the following smooth approximate function: with = 0.15.
The right panel of Fig. 1 shows the small error that we incur by making this step.We also replace the piece-wise v with its smooth equivalent by taking the derivative of the continuous f N (v) , such that: Finally, we note that there is a well-known link between QIF neurons and what is referred to as a -neuron which is described by a phase variable and can be used to describe excitable and spiking behavior.For a population of -neurons, the Kuramoto order parameter is defined as: where R is a measure of the degree of synchrony within the network and Ψ is the average phase of the population, R = 0 & R = 1 correspond to a fully asynchronous and fully synchronized states, respectively (Byrne et al., 2017).A conformal transformation has been reported that allows one to move from the plane of the mean firing rate and mean membrane potential of a QIF population to the complex plane of the Kuramoto order parameter.Defining W = r + iv , then the Kuramoto order parameter can be computed as: This allows us to have insight on the internal state of synchrony within the population and probe how it is affected with variation in the dynamical regime or parameter values, all the while staying on the level of a mesoscopic mathematical description of the population dynamics.

Numerical analysis
We numerically analyze the dynamics of the presented neural mass model for two example cases with typical sets of parameters.The bifurcation analysis as well as the numerical integration of the system of equations was performed using PyDSTool (Clewley et al., 2007), a time step of 0.1 was used for the timeseries simulation, and the power spectral density was obtained using the scipy.signal.periodogrampython function.

The case of excitatory regular spiking neurons
We analyze the resulting mean-field equations for the case of a population of excitatory regular spiking neurons with recurrent AMPAR synapses.
The equations for the Izhikevich regular spiking neuron are often written in the following form in the literature: The part 0.04V 2 + 5V + 140 came from fitting the spike initiation dynamics of a cortical neuron so that the membrane potential has mV scale, and the time has ms scale (Izhikevich, 2003).After factorization, this will correspond to the following parameter values for the original form of Eq. ( 1a): For the corresponding non- dimensional system, this corresponds to = 0.488.
Note that here we shift the slow current variables to add a V r term in the u equation for non-dimensionalization, such that we have a Accord- ingly, we add the non-zero external current term in the v equation to correct for the shift in u.
We consider the case in which the population receives external excitatory input that drives both AMPAR and NMDAR synapses.The system of equations becomes: (from now on, we drop the bar from ) J rec,A is the strength of the recurrent AMPAR synapses, whereas For r input = 0 , we fix Δ & J rec,A and analyze the stabil- ity of the equilibrium points as is varied and obtain the bifurcation diagrams shown in Fig. 2. The plots show the value of the mean firing rate and mean membrane potential corresponding to the fixed-point solution of Eqs.(4a)-(4d) (Fig. 3).For small enough , the system has a stable equilibrium point as a low constant firing rate that loses stability through a Hopf bifurcation giving rise to a limit cycle solution of multiscale oscillations (of similar nature as that shown in Fig. 4).For intermediate values of , the constant firing rate fixed point regains stability at high firing rate values, before losing stability again at a higher critical value of and leading way to a stable limit cycle solution of the fast spiking-like oscillation type (as that shown in Fig. 5).We then fix Δ, , J rec,A , J A & J N and vary r input .Figure 3 shows the bifurcation diagrams for three configurations: J N = 0 (AMPAR only), J N = 3 with linear NMDAR (the nonlinear Mg 2+ block is omitted), and J N = 3 with nonlin- ear NMDAR.In all cases, we have the following different dynamical regimes: stable fixed point with a low constant  4a)-( 4d) with = 0.01, Δ = 0.002, J rec,A = 6, J A = 6, J N = 3, r input = 0.06 ; the mean firing rate vs. time (r, top left), population synchrony measure vs. time ( R sync , top right), total synaptic current vs. time ( I synaptic , bottom left); power spectral density of the total synaptic current (bottom right); in all the plots, (-blue) is the case of nonlinear NMDA and (--green) is the case of linear NMDA firing rate, multiscale oscillation limit cycle, stable fixed point with high constant firing rate, and fast spiking-like oscillation limit cycle.With AMPAR only synapses, the latter regime persists for larger values of r input , that is, for the two cases of linear and nonlinear NMDAR, the curve is shifted to the left such that the fast spiking-like limit cycle loses stability at a significantly smaller value of r input .It can also be seen that in the presence of nonlinear NMDAR synapses, the regime of stable high firing rate fixed point is significantly reduced and shifted towards higher levels of r input when compared to the case of linear NMDAR.
In Fig. 4, we show a sample simulation time series for the mean firing rate, the synchrony measure R sync , and the total synaptic current I synaptic for sample parameter values in the regime of multiscale oscillations.
It can be seen that in the presence of nonlinear NMDAR, the oscillation is larger in amplitude with augmented slow and fast components.Figure 4 also shows the power spectrum density of the total synaptic currents indicating a downward shift in the main slow frequency and an upward shift in the dominant high frequencies.
Figure 5 shows another sample simulation data for parameter values in the regime of fast spiking-like oscillation.Here, the presence of nonlinear NMDA, instead, causes a smaller oscillation amplitude as compared to the case of linear NMDA, while still shifting frequency components of the total synaptic current towards lower values.

The case of inhibitory striatal neurons
The second example case is that of a population of inhibitory striatal medium spiny neurons (MSNs), the specific parameters are taken from (Humphries et al., 2009) in which an Izhikevich neuron model was fit to capture neurocomputational properties of MSNs as well as the effect of dopamine on MSNs as a linear scaling in the AMPAR and NMDAR maximal conductances.The allto-all coupled population of MSNs has recurrent inhibitory GABAR synapses and receives external excitatory AMPAR and NMDAR synapses.The resulting mean-field equations are: For the corresponding non-dimensional system, we have: = 0.629 and the coefficients for the f N (v) piece- wise fit are as follows:a 0 = 0.0306, Figure 6 shows the bifurcation diagram for the system of Eqs.(5a)-(5d), with r input = 0, as is varied.The system has a stable constant low firing rate stable equilibrium that loses stability through a Hopf bifurcation for large enough when fast spiking-like oscillations emerge.In Fig. 7, we fix Δ, , J rec,G , J A & J N and vary r input .It can be seen that the presence of nonlinear NMDAR allows for a regime of stable constant high firing rate solutions that do not exist for the case of linear NMDAR.
Moreover, the sample simulation time series in Fig. 8 illustrates how nonlinear NMDAR leads to oscillations with a larger amplitude for r input values that are below those of the high firing rate regime, along with a downward shift in the main frequencies of oscillation, as compared to the case of linear NMDAR.Whereas for r input values that are above those of the high firing rate regime, nonlinear NMDAR has the opposite effect, leading to much smaller amplitude oscillations than those of the case of linear NMDAR with the same r input value, along with a smaller downward shift in the frequencies, as shown in Fig. 9.
While an exhaustive parameter space exploration is beyond the scope of this work, it is worth noting that the effect of Nonlinear NMDA, as compared to the Linear case, was found to be qualitatively robust for a wide range of Δ & values (results not shown here).

Discussion
We presented a next generation neural mass model of Izhikevich-type neurons that includes NMDAR synaptic currents with nonlinear Mg 2+ block, along with APMAR and GABAR conductance-based synapses.We have shown how rewriting the nonlinear Mg 2+ block expression as a piece-wise polynomial maintains applicability of the Lorentzian ansatz in the mean-field model derivation.We have then performed numerical analysis of the resulting equations for two example cases.The first example involved a population of excitatory regular spiking neurons with recurrent AMPAR synapses and input AMPAR & NMDAR synapses.It was shown that the nonlinear NMDAR causes power spectral density of the total synaptic current (bottom right); in all the plots, (-blue) is the case of nonlinear NMDA and (-green) is the case of linear NMDA a shift in the constant high firing rate regime towards higher values of input and causes an amplification in amplitude of the multiscale oscillation solutions, while increasing the timescale separation by pushing the slow dominant frequencies downwards and the fast ones upwards.This is concurrent with larger intermittent dips in the population synchrony measure.Whereas, in the regime of fast spikinglike oscillations, nonlinear NMDAR exerts the opposite influence of decreasing the amplitude of oscillations in mean firing rate and those of the dips in population synchrony.The second example involved a population of inhibitory striatal medium spiny neurons with recurrent GABAR synapses and input AMPAR & NMDAR synapses.In this case, the nonlinearity in the NMDAR current gives rise to a constant high firing rate regime that is not present for the linear NMDAR case.This high firing rate regime separates two regimes of fast spiking-like oscillations, in which the effect of nonlinear NMDAR acts in opposite directions, for the lower input values, the amplitude of the fast oscillations is amplified due to the nonlinearity, while for the higher input values, the nonlinearity causes a significant decrease in the size of the oscillations and an associated curbing of the intermittent dip in population synchrony.
In summary, the nonlinearity of NMDAR synapses exerts a complex state-dependent effect on the population dynamics, that can manifest in opposite directions depending on neuronal type and level of input signal.This is in line with what is known on the complex state-dependent effect of dopamine action that utilizes the modulation of NMDAR maximal conductance as one of its routes of action.The presented model offers a computationally efficient building block for exploring the repercussions of complex NMDAR action, and associated neuromodulatory effects, on a mesoscopic and macroscopic brain level.By explicitly representing the different types of synaptic currents, not only with different timescales of action but also with the nonlinear state-dependent NMDAR feature, the model allows for the interrogation of mechanistic hypotheses on the role of imbalance in receptor or neuromodulator activity underlying mechanisms of brain disorders.That is, the formulation power spectral density of the total synaptic current (bottom right); in all the plots, (-blue) is the case of nonlinear NMDA and (-green) is the case of linear NMDA suggested here will permit going beyond the common simplifying assumptions of phenomenological sigmoidal input-output response functions, to a more biophysically grounded approach that expresses neuromodulatory action as a modulation of maximal conductance of the different synaptic currents involved.More specifically, the sigmoidal "firing rate vs. input current" response function, commonly used in classical mean-field models, is best suited to describe constant population firing rates in response to an input current that drives a near-constant deviation of membrane potential from its resting value (Abbott & Chance, 2005), a condition which is violated in the dynamical regime of oscillatory firing rate solutions representing complex burstlike activity as observed in the BG nuclei (Bergman, 2021).Moreover, the pharmacological sigmoidal function, recently proposed in the literature to account for neuromodulatory currents in whole brain network models (Joshi et al., 2017;Kringelbach et al., 2020), can only capture effects that enter as an additive current in target regions and falls short of capturing state-dependent parametric (i.e.multiplicative) effects as those reported in the case of dopaminergic action.
Our formulation allows for coupling multiple populations of different types of neurons, such as those of the BG nuclei, and opens the path to realistically investigate the effect of different dopamine levels as manifesting in a scaling of different maximal conductances to match known distribution of different receptor types (Lindahl & Kotaleski, 2016).Building the BG circuit using the proposed mean-field model will permit its embedding in the larger full brain network through connectome based whole brain models (Breakspear, 2017;Sanz Leon et al., 2013), and thus facilitate the study of emergent effects on a whole brain level, while preserving computational efficiency by avoiding spiking neurons simulation that would be required for detailed BG modeling (Humphries et al., 2006(Humphries et al., , 2018;;Meier et al., 2022).
The presence of a conformal mapping to go from mean firing rate and mean membrane potential to the Kuramoto order parameter variables offers the added advantage of probing alterations in average internal synchrony within regions while remaining on the level of mesoscopic description, this is particularly of relevance to disorders like PD in which the pathology manifests not only in changes in mean firing rate but also in that of coherence within the target neuronal populations (Bergman, 2021).
More broadly, NMDARs activity plays a central role in brain plasticity, learning and memory, the disruption of which is associated with impairment seen in a wide range of pathologies, such as Alzheimer's disease, amyotrophic lateral sclerosis, Huntington's disease, Parkinson's disease, schizophrenia and major depressive disorder (Adell, 2020;Anticevic et al., 2012).We hope that the model proposed here will contribute to better representation and inclusion of nonlinear NMDAR action in the growing efforts on multiscale brain modeling, aiming to disentangle mechanisms of brain disorders, by bridging the gap between microscopic and macroscopic neuronal dynamics (D'Angelo & Jirsa, 2022;Deco & Kringelbach, 2014;Jancke et al., 2022;Joshi et al., 2017;Rolls et al., 2008;Shine et al., 2021).

Fig. 10
Time series for the mean firing rate of the spiking neuron network; blue is for the nonlinear NMDA case; green line is for the linear NMDA case.Neuron parameters are the same as in the case of excitatory regular spiking neurons presented in the Section 3.1, with = 0.01, Δ = 0.002, J rec = 6, J A = 6, J N = 3, r input = 0.06 Fig. 11 Time series for the mean firing rate of the spiking neuron network (-blue) compared to the mean firing rate timeseries from the integration of the mean-field model (--red).Neuron parameters are the same as in the case of excitatory regular spiking neurons presented in the Section 3.1, with = 0.1, Δ = 0.002, J GABA = 1, J AMPA = 1.4,J NMDA = 0.7, r input = 0.06 J A & J N are the strengths of the input AMPAR & NMDAR synapses and the corresponding reversal potentials are E A = E N = 0 , such that e A = e N = 1 .The synap- tic time constants for AMPAR & NMDAR are taken to be 6ms & 160ms, respectively, corresponding to A = 19.83& N = 529.For V r = −82.66, the coefficients for the f N (v) piecewise fit are as follows: (we set v cut = −1.0). a 0 = 0.027, a 1 = 0.106, a 2 = 0.089 , b 0 = −0.396,b 1 = 1.559, b 2 = 1.158, c 0 = 1.038, c 1 = −1.018 ,along with v 0 = 0.582, v 1 = 1.112.

Fig
Fig. 3 Bifurcation diagram for Eqs.(4a)-(4d), with Δ = 0.002, = 0.01, J rec,A = 6 .The plots show the value of mean firing rate (r, left) and mean membrane potential (v, right) of the fixed points of the system as a function of the external input firing rate

Fig. 12
Fig. 12 Time series for the mean firing rate of the spiking neuron network; blue is for the nonlinear NMDA case, green line is for the linear NMDA case.Parameter values as in The case of inhibitory striatal neurons presented in the Section 3.2, with = 0.1, Δ = 0.002, J GABA = 1, J AMPA = 1.4,J NMDA = 0.7, r input = 0.08

Fig. 13
Fig. 13 Time series for the mean firing rate of the spiking neuron network (-blue) compared to the mean firing rate timeseries from the integration of the mean-field model (--red).Parameter values as in The case of inhibitory striatal neurons presented in the Section 3.2, with = 0.1, Δ = 0.002, J GABA = 1, J AMPA = 1.4,J NMDA = 0.7, r input = 0.08

Fig. 14
Fig. 14 Time series for the mean firing rate of the spiking neuron network (-blue) compared to the mean firing rate timeseries from the integration of the mean-field model (--red).Parameter values as in the case of inhibitory striatal neurons with