A mean-field approach to the dynamics of networks of complex neurons, from nonlinear Integrate-and-Fire to Hodgkin-Huxley models

We present a mean-field formalism able to predict the collective dynamics of large networks of conductance-based interacting spiking neurons. We apply this formalism to several neuronal models, from the simplest Adaptive Exponential Integrate-and-Fire model to the more complex Hodgkin-Huxley and Morris-Lecar models. We show that the resulting mean-field models are capable of predicting the correct spontaneous activity of both excitatory and inhibitory neurons in asynchronous irregular regimes, typical of cortical dynamics. Moreover, it is possible to quantitatively predict the populations response to external stimuli in the form of external spike trains. This mean-field formalism therefore provides a paradigm to bridge the scale between population dynamics and the microscopic complexity of the individual cells physiology. NEW & NOTEWORTHY Population models are a powerful mathematical tool to study the dynamics of neuronal networks and to simulate the brain at macroscopic scales. We present a mean-field model capable of quantitatively predicting the temporal dynamics of a network of complex spiking neuronal models, from Integrate-and-Fire to Hodgkin-Huxley, thus linking population models to neurons electrophysiology. This opens a perspective on generating biologically realistic mean-field models from electrophysiological recordings.

In this article we present a general approach to determine the transfer function for complex models, from the Adaptive Exponential Integrate-and-Fire (AdEx) to the Hodgkin-Huxley (HH) 42 and the Morris-Lecar (ML) models. As a result, we obtain mean-field equations for the population 43 dynamics in AI regimes as observed in cortical regions for highly detailed models, creating a bridge 44 between electrophysiology at the microscopic scale and the details of the famous transfer function 45 first used by Wilson and Cowan as a sigmoid. 46 Finally, we test not only the ability of our mean-field models to describe spontaneous activity of the considered neuronal populations, but also their predictive power for network response to 48 external stimuli. We show that, provided the stimuli are fairly slow, the mean field model gives 49 good quantitative predictions. 50 2 Materials and Methods 51 We describe here the neuronal and network models used in this study. We also introduce mean-field 52 equations describing population dynamics and the template to estimate the transfer function that 53 we apply to all the neuronal models under consideration. Adaptive Exponential Integrate-and-Fire and Morris-Lecar). The dynamics of each node k follows: wherex and F (x) represent the neuronal state and dynamics, the latter depending on the 61 specific model (see Section 2.2). Note the notationx k which indicates that, in general, each neuron 62 is characterized by a vector of variables. The synaptic current impinging on the postsynaptic 63 neurons k, I syn , is modeled as: where E e = 0 mV (E i = −70 mV) is the excitatory (inhibitory) reversal potential and G (e,i) syn is 2.2 Single neuron models 74 We describe here the neuronal models used in the rest of the paper, starting from the Integrate- The dynamics of each of the AdEx neurons i is described by the following 2D (herex i = (v i , w i )) 78 differential equations (Brette and Gerstner, 2005) : where c m = 150pF is the membrane capacity, v i is the voltage of neuron i and, whenever v i > The dynamics of the Hodgkin-Huxley model (Hodgkin and Huxley, 1952) is given by the following 5-dimensional system of differential equations (Pospischil et al., 2008): with the gating functions, where v i is the voltage and (n i , m i , h i , p i ) are the corresponding gating variables of the i-th neuron.
We set the spike emission times t sp (k) for this model to time steps in which the membrane potential 93 v exceeded a voltage threshold of 10 mV. Unless stated otherwise, the membrane capacitance 94 c m = 200 pF/cm 2 , the maximal conductance of the leak current g L = 10 mS/cm 2 , the sodium 95 current g N a = 20 mS/cm 2 , the delayed-rectifier potassium current g K = 6 mS/cm 2 , the slow  The dynamics of the Morris-Lecar model (Morris and Lecar, 1981) is described by the system of 102 differential equations: where c m = 2µ F/cm 2 is the membrane capacitance, v i is the membrane potential in mV, N i and M ss are the fraction of open potassium and calcium channels, respectively. The current I 0 = 0.2 nA/cm 2 is a reference DC external current. Spike emission times are established in the same way as for the HH model. The maximal conductances for the leakage current (L), calcium (Ca) and potassium (K) were fixed to g L = 20 mS/cm 2 , g Ca = 80 mS/cm 2 and g K = 160 mS/cm 2 , respectively. The reversal potentials are E L = −50 mV for excitatory RS neurons and E L = −70 mV for inhibitory FS neurons, E Ca = 120 mV and E K = −84 mV. The quantities M ss and N ss are modeled as: with where V 1 = −1.2 mV, V 2 = 18 mV, V 3 = 2 mV, V 4 = 30 mV are tuning parameters that 104 determine the half activating voltage and slope of the activation curves for calcium and potassium 105 conductances. This choice of parameters is such that the ML neuron is set in a type II excitability 106 class, i.e. its response to a DC current is discontinuous and the neuron firing rate increases very 107 slowly with the injected current (data not shown). the asynchronous irregular regimes here investigated. The differential equations read: where µ = {e, i} is the population index (excitatory or inhibitory), ν µ the population firing rate  these assumptions the neuronal output firing rate F ν is given by the following formula: used α = 2. In the following section we introduce how the quantities µ V , σ V , τ V can be expressed 152 as functions of the presynaptic excitatory and inhibitory firing rates ν E and ν I .
153 Table 1: Fit Parameters AdEx Neurons (in mV ) where K i,e is the average input connectivity received from the excitatory or inhibitory population 159 (in our cases typically K e = 400 and K i = 100) and in our model τ e = τ i = τ (see Eq. 3).

160
The mean conductances will control the total input of the neuron µ G and therefore its effective dynamics of the currents coming into play at the spiking time (e.g. sodium channels dynamics or 164 the exponential term of the AdEx model). We thus consider, for all neurons, only the leakage term 165 and the synaptic input in order to estimate subthreshold moments. Accordingly, we can write the 166 equation for the mean subthreshold voltage: The final formulas for σ V and τ V follow from calculations introduced in Zerlaut et al. (2018) and they read: show in the following sections, the fitting procedure will account for discrepancies in the actual 171 evaluation of voltage moments by permitting an accurate prediction of neuron output firing rate. threshold was taken as a second order polynomial in the following form: where we introduced the quantity τ N V = τ V G l /c m . We evaluated {P } through a fit according to 179 simulations on single neurons activity setting first µ 0 We present here the results of a comparison between mean-field predictions and direct simulations.

228
Notice that the ML model shows a decrease of firing rate at frequencies higher than 8 Hz (i.e. no  In order to complete the comparison between the mean-field model and the network dynamics, 275 we study the response of the system to external stimuli. In particular, we consider an incoming where Θ is the Heaviside function, and T 1 and T 2 are the rise and decay time constants, 279 respectively. In Fig. 4 we report the comparison between the mean-field prediction and the network 280 dynamics. 281 Figure 4: Population response to external stimuli: AdEx, HH and ML Top panels show the raster plot for excitatory (green dots) and inhibitory (red dots) neurons in response to an external excitatory stimulus (black dashed line in lower panels). Lower panels show the corresponding population rate (noisy line) together with the mean and standard deviation over time predicted by the the second order mean field model (red for inhibition and green for excitation). Superimposed the result obtained for the mean-field at the first order (black dots), which are almost coincident with results at the second order. Left column is obtained for the AdEx model, middle column for the HH model and right column for the ML model. Parameters are the same as in Fig. 3 and the external input (see Eq. (26)) has parameters A = 2 Hz, T 1 = 100 ms, T 2 = 150 ms for AdEx and HH and A = 2 Hz, T 1 = 100 ms, T 2 = 150 ms for ML, with t 0 = 2 s.
By looking first at the AdEx and HH models, we observe that both mean-field models under 282 investigation compare favorably with their corresponding network dynamics. We also verified, as 283 it has been shown in (di , that the faster the input dynamics is, the worse the 284 agreement becomes. Indeed, for the Markovian hypothesis to hold, we need the time scale T to be 285 much larger than the autocorrelation time in the spontaneous activity T ∼ τ m ∼ 10 ms.

286
Considering now the case of the ML model, we observe that by looking at Fig. 2 a relatively   287 strong input would bring single neurons to a depolarization block, which appears at relatively low 288 activity levels. According to this difference with respect to AdEx and HH models, we would expect 289 the population dynamics to show different properties in response to external perturbations. Indeed, 290 as reported in Fig. 4c, the response to an external stimulus is very different from the one observed 291 in the HH and AdEx models. In fact, in this case the excitatory stimuli turns out to inhibit both 292 population activities. This anti-correlation between population input and output is well captured 293 in its time course also by the mean-field model. This result shows that also for a more complex Finally, we compare the results of the first and second order mean field on average population 298 rates. In Fig. 4 we superimpose the continuous green (red) line for excitatory (inhibitory) rate 299 obtained with the second order mean file with the results obtained with the first order (black dots). 300 We observe that the two quantities almost overlap (the difference is too small to be appreciated 301 at this scale). Nevertheless, it is worth noticing that the second order mean field permits to 302 obtain non-trivial information on the population dynamics and its fluctuations in time, with good 303 quantitative predictions of the covariance of population rates (see the histograms in Fig. 3 and 304 shadows in Fig. 4)).

306
In this manuscript, we have reviewed a formalism to derive mean-field models from networks of 307 spiking neurons and we have applied it to different complex neuronal models. The key to derive 308 such "biologically realistic" mean-field models is to be able to obtain the transfer function of 309 complex neuronal models. The approach we followed used a mean-field formalism based on a It is important to note that we limited here to "simple" firing patterns, i.e. neurons fire tonically 319 in response to an external stimulus. In this setup the transfer function is well defined as the neuron's 320 firing rate defines completely the spiking pattern. In cases where neurons exhibit different kind 321 of activity, e.g. bursting, a different approach needs to be employed (see (Ostojic and Brunel,322 2011)). Nevertheless, in the context of tonic neuronal activity the method is shown to be able to 323 capture the response function of highly realistic models. We have studied here the predictions of level, as the case of focal seizures. 345 We also reported that, in the framework of the Asynchronous regimes here considered, correc- current, is a stimulating perspective for future works and will permit to obtain mean-field models 361 for realistic neuronal models beyond the asynchronous irregular regime.

362
Moreover, beyond the input-output transfer function used here, a more complex transfer func-363 tion has been used in order to take into account other features of neuron response, e.g. response 364 in frequency (Ostojic and Brunel, 2011). The addition of variables to account for a richer spiking 365 pattern is an interesting direction, in case one is interested in modeling brain regions characterized 366 by non-tonic firing of neurons (e.g. bursting cells in the thalamus). The general framework pre-367 sented here could be extended in this direction, as it has been done to account for spike frequency 368 adaptation yielding slow oscillations at the population scale.

369
Another possible extension is to apply the same formalism to complex neuronal models that 370 include dendrites. A first attempt has been made in this direction (Zerlaut and Destexhe, 2017b) by 371 considering simple "ball and stick" neuron models, where some analytic approximation is possible.

372
In principle, it should be possible to apply this approach to models based on morphologically-373 reconstructed neurons, and calculate the transfer function of such models (work in progress).

374
This will lead to mean-field models based on morphologically realistic neuronal models. However,