Exact reduction methods for networks of neurons with complex dynamic phenotypes

Collective dynamics of spiking networks of neurons has been of central interest to both computation neuroscience and network science. Over the past years a new generation of neural population models based on exact reductions (ER) of spiking networks have been developed. However, most of these efforts have been limited to networks of neurons with simple dynamics (e.g. the quadratic integrate and fire models). Here, we present an extension of ER to conductance-based networks of two-dimensional Izhikevich neuron models. We employ an adiabatic approximation, which allows us to analytically solve the continuity equation describing the evolution of the state of the neural population and thus to reduce model dimensionality. We validate our results by showing that the reduced mean-field description we derived can qualitatively and quantitatively describe the macroscopic behaviour of populations of two-dimensional QIF neurons with different electrophysiological profiles (regular firing, adapting, resonator and type III excitable). Most notably, we apply this technique to develop an ER for networks of neurons with bursting dynamics.


I. INTRODUCTION
For decades, theoretical neuroscientists have been using mean-field theory to reduce the description of neural circuits composed of many interacting neurons to a lowdimensional system that describes the macroscopic dynamical states of the network. Different techniques have been employed in recent years based on approximations of the underlying network dynamics [1,4,5,15,20]. This approach is quite powerful as it generates a reduced picture of the neural population that can be used to study how brain functions arise from the collective behaviour of spiking neurons. Recently, an exact reduction method based on the Ott Antonssen ansatsz [17] allows to derive exact mean-field equations for heterogenaous and globally coupled networks of quadratic integrate-and-fire (QIF) neurons [14] . This methodology has been employed to study networks with electrical synapses [13], delays [19], working memory [21], to analitically estimate cross-frequency-couplings [6] and recently to study brain activity at the whole brain level [9]. Moreover, a recent work shows that it is possible to derive mean field equations also for sparse networks in presence of noise [10]. However, this approach is limited to networks of onedimensional quadratic integrate-and-fire (QIF) neurons that cannot account for complex spiking and bursting dynamics. Two-dimensional quadratic integrate-and-fire models (e.g., the Izhikevich neuron model [11]) reproduce a wide variety of spiking and bursting behaviours. Com-ing up with reduced descriptions of networks that cover the various Izhikevich neuron dynamical profiles has been a largely open question. An exact mean-field reduction of such neuron models will enable us to derive the macroscopic dynamics of populations of neurons with any kind of spiking properties. In this work, we propose a new reduction method for large networks of conductance-based interacting Izhikevich neurons. We start by presenting the two-dimensional neuron model and then show how a separation of time scales of the variables describing the state of the neurons allows us to explicitly solve the continuity equation of the system, which represents the evolution of the state of the neural population. By doing so, we obtain a system of two coupled variables, the firing rate and the mean voltage, which together describe the evolution of the macroscopic system. We study the accuracy of the neural mass by comparing its predictions with numerical simulations of finite-size networks of neurons with different spiking properties. This approach opens the possibility of generating realistic mean-field models from electrophysiological recordings of individual neurons and can be used to relate the biophysical properties of neurons with emerging behaviour at the network scale.

II. POPULATION MODEL OF IZHIKEVICH NEURONS
We derive a mean-field model for populations of coupled Izhikevich neurons. Each neuron i from a population W is described by a fast variable representing the membrane potential, V (mV ), and a slow variable representing the recovery current, u (pA): where the onset of an action potential is taken into account by a discontinuous reset mechanism: The parameters are as follows: C m stands for the membrane capacitance, V r is the resting potential, V t the threshold potential, a is a scaling factor, α the time constant of the recovery variable u, β modulates the sensitivity of the recovery current to subthreshold oscillations, and I i is the total current acting on neuron i. We consider I i = η i + I ext + I syn,i . The parameter η i represents a background current. To account for the network heterogeneity, the parameter η i is randomly drawn from a Lorentzian distribution with half-width ∆ centered at η, g(η) = 1 π ∆ (η−η) 2 +∆ 2 . I ext is an external current acting on population W (identical to all neurons). I syn,i is the total synaptic current acting on neuron i given by: where E Z r is the reversal potential of the synapse, and s W Z the synaptic conductance. If we assume that all neurons of population W are connected to all neurons of population Z, the synaptic conductance s W Z can be described according to the following equation: where δ is the Dirac mass measure and t k f is the firing time of neuron k. The parameter τ s represents the synaptic time constant, N Z is the number of neurons of population Z, and p W Z is the coupling strength of population Z on population W .
A. Adiabatic approximation of the two-dimensional QIF neuron model We exploit the time scales difference between the dynamics of the membrane potential V and the recovery variable u to reduce the dimensionality of the neural network. If the time scale of the recovery variable is much slower than the other variables, we can invoke an adiabatic approximation by considering that all neurons of population W receive a common recovery variable u W . This results in the modified Izhikevich QIF model: where V W is the mean membrane potential of the population W , described as follows: Note that we have incorporated the resetting mechanism of the variable u W in the last term of equation (5). The onset of an action potential is now described by: From now on we will consider equation (4) written in terms of the parameters b = a(−V r − V t ) and c = aV r V t : The main consequence of the adiabatic approximation is the reduction in the number of state variables describing a neuron in the population from (V W i , u W i ) to (V W i ). This is a crucial step in our method since it enables us to solve the continuity equation of the system analytically, as we demonstrate in the next section.

III. MEAN-FIELD REDUCTION
Having reduced the dimensionality of the Izhikevich neuron model, in the adiabatic approximation we have that the effective current acting on neuron i is I i → I i + u W . As a result, we can employ the reduction method developed in [14] for a one-dimensional QIF neuron model, thus obtaining a macroscopic description of a population of Izhikevich neurons. The reduced model obeys the following differential equations (see appendix A for more details): Here r W is the population firing rate, v W the mean membrane potential across neurons and u W the mean slow recovery variable (e.g. adaptation)

A. Numerical simulations
The Izhikevich two-variable QIF model can, with the adequate choice of parameters, reproduce many of the key features of firing patterns observed in neurons, such as tonic spiking, subthreshold oscillations, and bursting [11]. We here compare the neural mass model of Eq.s 8 with direct simulations of QIF neurons in different parameter regimes characterized by different intrinsic firing dynamics of neurons. Figure 1 illustrates a comparison of the dynamics of the full network of Izhikevich QIF neurons with its corresponding reduced system. Regarding the full system, each population is made up of N = 3000 neurons. The neurons are described by the two-dimensional QIF model 1, with the respective parameters specified in Table I. The firing rate is calculated according to: . For the reduced system, the firing rate is calculated according to equation A16a. The reduced description closely follows the firing activity of the full network for all populations.

B. Limitations of reduction formalism
One crucial assumption of the mean-field reduction formalism here presented is the slow dynamics of the recovery variable u. However, every time a neuron reaches V peak , the membrane potential V is reset to V reset and the recovery variable u is instantaneously increased by u jump , adding a discontinuity to the system. Given that the adiabatic approximation made in section II A relies on the assumption that the variable u is the slowest variable in the system and that therefore we can consider that all the neurons in the population receive a variable u with approximately the same value, adding an instantaneous jump invalidates the approximation made. The higher the jump, the more evident this is. In the examples considered in Figure 1, u jump is sufficiently small for the reduction to work with minimal errors. This means the reduced system here derived can be used to study the activity of populations with any of the spiking dynamics portrayed in Figure 1. Still, it may not be adequate to describe the activity of certain neural populations, such as rat spiny projection neurons of the neostriatum and basal ganglia, for which the value of u jump used to describe their dynamics induces imprecisions [12].
In Figure 2 we compare the full and reduced system for a population of uncoupled spiny projection neurons (u jump = 150 pA) [12]. We then systematically decrease the value of u jump we see how that changes the accuracy between dynamics of the full and reduced system. All the neurons receive an external current as described in Figure 1 (D). We see that there is not a perfect agreement between the full and reduced system for a population of spiny projection neurons (left panel). Decreasing the value of u jump notably improves the agreement between the full and reduced system significantly, confirming that the high u jump is at the origin of the mismatch observed. For u jump = 150, one way to improve the representation of the population activity would be to decrease ∆. By decreasing the variance ∆ of the intrinsic variable η, we decrease the heterogeneity of the network. As a result, we can consider again that all the neurons in a population W are receiving the same variable u at any given time (simulations not shown).

The particular case of bursting neurons
A critical point of the derivation of our reduced meanfield model is the assumption that the firing rate of a population is defined as the flux at infinity. In other Results were obtained using the Izhikevich two-dimensional QIF neuron model [11] with the adequate choice of parameters (see Table I  words, we consider V peak → ∞. Similarly, we assume that V reset → −∞. While moving V peak towards infinity does not change the intrinsic spiking properties of the neurons that constitute the population, moving V reset in the direction of −∞ changes the microscopic behavior of bursting neurons. Figure 3 (B) depicts the phase portrait of an intrinsically bursting neuron. Starting at point A, we are on the V-nullcline, where by definition dV dt = 0, and the dynamics is going to be governed by the u-component. Since we are on the left of the u-nullcline, the trajectory follows a downward flux. As u slowly decreases, we reach point B below the V-nullcline, and the fast dynamics in the V direction pushes the system towards V peak , at which point the system is reset to V reset . This last process repeats while u slowly increases until it reaches point C, where a voltage reset takes the system to a point above the V-nullcline. In this region, the flux is directed towards the left, which brings the system back to point A. By decreasing the value of V reset , we lose the bursting dynamics, and the neuron model now shows tonic spiking instead (see Figure 3 (C)). One way to preserve the bursting dynamics of the microscopic system would be to move the V and u-nullclines by the same amount as the V reset (Figure 3 (D)). We do so by decreasing the values of V r and V t (remember that b = −a(V r + V t ) and c = aV r V t ). From Figure 3 (E), we see that by adopting this change the full and reduced system activity have approximately the same shape, but that they do not perfectly agree. It is important to note, however, that this method presents important faults: it implies that at V reset → −∞, the resting and threshold potential, V r and V t , should also move to −∞. This is not only a problem from the biological point of view, but it can also invalidate some mathematical results adopted during the derivation of the mean-field reduction; namely, when solving explicitly the integrals that define the firing rate and mean voltage of the population (equations (A8) and (A9)).

IV. DISCUSSION
In this work, we presented a reduction formalism that allows us to predict the collective dynamics of large networks of conductance-based interacting spiking neurons with different spiking properties. Starting with a population of Izhikevich neurons, we considered an adiabatic approximation of a slow recovery variable, which enabled us to reduce the dimension of variables that describes the state of the network. We then applied the Lorentzian ansatz to solve the continuity equation and reduce our full network to a low-dimensional macroscopic system. This mean-field formalism provides a paradigm to bridge the scale between population dynamics and the microscopic complexity of the physiology of the individual cells, opening the perspective of generating biologically realistic mean-field models from electrophysiological recordings for a variety of neural populations.
A similar adiabatic approach appears in [5] and [15]. Di Volo and colleagues propose a mean-field model of spiking neurons with recovery variable by calculating the transfer function (i.e. neurons' stationary firing rate in response to external spike trains) in a semi-analytical way. This approach, however, assumes that neuron dynamics has a stationary firing rate in response to external spike trains, and it does not allow to study populations of neurons whose transfer function is not completely defined with only one variable (i.e. neurons' stationary firing rate) [5]. Similar to our approach, Nicola and Campbell use moment closure and a steady-state approximation of the recovery variable to write an expression for the population firing rate, defined as the integral of the population density function. However, they cannot apply the Lorentzian ansatz to solve the integral because they don not consider the heterogeneity of the population. Therefore, for some types of networks it won't be possible to be evaluated explicitly the firing rate integral [15]. In the context of QIF neurons with a spike frequency adaptation, an adiabatic approach has been employed in [8], while in [3] moment closure is used to get a neural mass model. Nevertheless, in these reduced models the neurons are not intrinsically bursting as in the case of the Izhikevic model here considered. Other adiabatic approaches in this field are appearing, e.g. to consider slow variables modeling conductance-based ion exchange [2].
Sufficient requirements for our approach to be valid are that the recovery variable u is the slowest in the system and that u jump is relatively small. This means that even though it is possible to describe any class of spiking dynamics, the reduced model might be unable to describe the activity of specific neural populations, such as spiny projection neurons of the neostriatum and basal ganglia whose models require a rather high value of u jump [12].
It is important to note that even though the original Izhikevich neuron model for tonic bursting neurons fills all the requirements, the mean-field system seems to be inadequate to describe the population's behavior. Since, in the Izhikevich two-dimensional QIF model, the bursting mechanism depends on the position the system acquires in the phase-plane (V ,u) upon reset, when moving V reset to inf ty will alter the behavior of the microscopic system. Therefore, despite having a good agreement between the full and reduced system, the population at study is no longer a population of tonic bursting neurons but of tonic spiking neurons. A solution found was to move the u and V-nullclines with V reset , so that an action potential will reset the system to the same position in the phase-plane relative to the nullclines and preserve the bursting mechanisms of the original model. We do so by decrease the resting and threshold potential, V r and V t by the same amount as V reset . The full and reduced system of the resultant tonic bursting neurons does not perfectly agree, but it accurately reproduces the oscillatory behavior of the population. In other words, when the system receives a big enough external input I ext , both the full and reduced system show damped oscillations but with different frequencies. Therefore, the mean-field description may still be useful to study certain features of a population of bursting neurons and qualitative behavior. However, it is important to note that the approach taken for the case of the bursting neurons presents some fundamentals problems. Namely, it implies that both the reset, resting and threshold potential are set to −∞. An alternative solution is to consider the two-dimensional theta neuron model with a slow recovery variable, which with the appropriate choice of parameters can produce bursting [7], and apply the steps as for the derivation of a two-dimensional QIF. In the theta neuron model, the system evolves along a circle and V ∈ [−∞, +∞] maps to θ ∈ [0, 2π]. A deeper analysis of this approach is necessary to prove its validity. the net fraction of trajectories per time unit that crosses the value V , we can write the continuity equation which expresses the conservation of the number of neurons. Note that in integrate-and-fire models, the number of trajectories is not conserved at V = V reset and V = V peak . By taking V reset and V peak to infinity, we ensure that the boundary conditions are the same and that the number of trajectories is conserved [16]. According to the Lorentzian ansatz (LA) [14], solutions of the continuity equation (A1) for a population of QIF neurons converge to a Lorentzian-shaped function with half-width x(η, t) and center at y(η, t) of the form: We discuss the validity of the LA here applied in appendix B. Here, x(η, t) and y(η, t) are statistical variables that represent the low dimensional behavior of the probability density function ρ. Adopting the LA, we obtain the low dimensional system: that can be written in the complex form as: with w(η, t) = ix(η, t) + y(η, t)

The macroscopic variables: firing rate and mean voltage
The firing rate is obtained by summing the flux for all η at V = V peak . Taking V peak → ∞ the firing rate of a population W is defined as follows The mean voltage of the population is obtained by integrating the probability density function ρ for all V and η values: Adopting the solution for the continuity equation (A2) and inserting it into equations A6 and A7 we have that the phenomenological variables x and y relate with the firing rate, r, and mean voltage, v, as follows: To avoid indeterminancy of the improper integral, we resort to the Cauchy principal value to evaluate the integral A9 (p.v.
In the case of a Lorentzian distribution, the principal value is given by x (x−x0) 2 +σ 2 dx = x 0 . We then have that the mean voltage is given by: As previously mentioned, in the mean-field limit, the probability distribution function g(η) is given by The distribution g(η) has poles at η − i∆ and η + i∆, and can be written as The integrals in equations A8 and A11 are evaluated by closing the integral contour in the complex η plane and using the residue theorem. We then have that the firing rate and mean potential relate to the Lorentzian coefficients x and y according to the following expression: r W (t) = a C m π x(η ± i∆, t) (A12) v W (t) = y(η ± i∆, t) (A13) Given equations A12 and A13 and noting that we have that the continuity equation reduces to the low-dimensional macroscopic dynamical system: Since the firing rate always has to be non-negative, we needed to evaluate the closed integral contour containing the pole of g(η) in the lower half of η plane, i.e., η − i∆. Until now, we considered the integral contour in both the upper and lower half of the η. This is because the Lorentzian variables x and y have no physical meaning. Therefore we could not make any conclusions regarding which contour to consider when using the residue theorem to solve (3.15) and (3.18) until now. We have that an exact mean-field reduction of a population of interacting conductance-based Izhikevich two-dimensional QIF neurons is given by : Therefore, in the original variable V X , our system converges to: After some algebraic manipulations, we recover the LA (A2) ρ(V W |η, t) = 1 π x(η, t) (V W − y(η, t)) 2 + x(η, t) 2 (B8) The LA ansatz solves the continuity equation exactly, making the system amenable to theoretical analysis. In section III A, we show that these solutions agree with the numerical simulations of the original QIF neurons, further validating the application of the LA.