# A Mean Field to Capture Asynchronous Irregular Dynamics of Conductance-Based Networks of Adaptive Quadratic Integrate-and-Fire Neuron Models

Christoffer G. Alexandersen, Chloé Duprat, Aitakin Ezzati, Pierre

Houzelstein, Ambre Ledoux, Yuhong Liu, Sandra Saghir, Alain Destexhe, Federico Tesler, Damien Depannemaecker

► To cite this version:

Christoffer G. Alexandersen, Chloé Duprat, Aitakin Ezzati, Pierre Houzelstein, Ambre Ledoux, et al.. A Mean Field to Capture Asynchronous Irregular Dynamics of Conductance-Based Networks of Adaptive Quadratic Integrate-and-Fire Neuron Models. Neural Computation, 2024, 36 (7), pp.1433-1448. 10.1162/neco_a_01670. hal-04649199

HAL Id: hal-04649199

https://hal.science/hal-04649199v1

Submitted on 16 Jul 2024

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Mean Field to Capture Asynchronous Irregular Dynamics of Conductance-Based Networks of Adaptive Quadratic Integrate-and-Fire Neuron Models

## Christoffer G. Alexandersen

Mathematical Institute, University of Oxford, OX2 6GG, Oxford, U.K.

## Chloé Duprat

chloe.duprat@epfedu.fr

Paris-Saclay University, Institute of Neuroscience, CNRS, 91400 Saclay, France, and Institut de Neurosciences des Systèmes, Aix-Marseille University, INSERM, 13005 Marseille, France

## Aitakin Ezzati

aitakin.EZZATI@univ-amu.fr

Institut de Neurosciences des Systèmes, Aix-Marseille University, INSERM, 13005 Marseille, France

## Pierre Houzelstein

pierre.houzelstein@gmail.com

Group for Neural Theory, LNC2, INSERM U960, DEC, École Normale

Supérieure—PSL University, 75005 Paris, France

## Ambre Ledoux

ambre.ledoux35@gmail.com

Paris-Saclay University, Institute of Neuroscience, CNRS, 91400 Saclay, France

## Yuhong Liu

yuhong.liu@uni-mainz.de

Institute of Physiological Chemistry, Johannes Gutenberg University of Mainz,

55128 Mainz, Germany, and Institute of Experimental Epileptology and Cognition

Research, University of Bonn Medical Center, 53127 Bonn, Germany

## Sandra Saghir

sandra.saghir@campus.tu-berlin.de

Department of Software Engineering and Theoretical Computer Science, Technische

Universität Berlin, 10623 Berlin, Germany

<footnote>Christoffer G. Alexandersen, Chloé Duprat, Aitakin Ezzati, Pierre Houzelstein, Ambre Ledoux, Yuhong Liu, and Sandra Saghir are equally contributing first authors. Federico Tesler and Damien Depannemaeker are equally contributing last authors.</footnote>
Alain Destexhe
Federico Tesler

ftesler@gmail.com

Paris-Saclay University, Institute of Neuroscience, CNRS, 91400 Saclay, France

Damien Depannemaecker

damien.DEPANNEMAECKER@univ-amu.fr

Institut de Neurosciences des Systèmes, Aix-Marseille University, INSERM,
13005 Marseille, France

Mean-field models are a class of models used in computational neuro-
science to study the behavior of large populations of neurons. These mod-
els are based on the idea of representing the activity of a large number of
neurons as the average behavior of mean-field variables. This abstraction
allows the study of large-scale neural dynamics in a computationally effi-
cient and mathematically tractable manner. One of these methods, based
on a semianalytical approach, has previously been applied to different
types of single-neuron models, but never to models based on a quadratic
form. In this work, we adapted this method to quadratic integrate-and-
fire neuron models with adaptation and conductance-based synaptic in-
teractions. We validated the mean-field model by comparing it to the
spiking network model. This mean-field model should be useful to
model large-scale activity based on quadratic neurons interacting with
conductance-based synapses.

# 1 Introduction

Modeling brain activity over different scales is a relevant challenge. Mul-
tiple models and approaches have been proposed over the years, going
from subcellular to whole-brain scales, to serve various purposes and ap-
plications. In order to model the mesoscopic scale in particular, one option
is to build phenomenological neural-mass models describing observations
made at this scale. Another alternative is a bottom-up approach where the
dynamics of the mesoscopic scale are derived by developing a mean-field
model of the microscopic scale (i.e., of spiking neural network models).
The mean-field approximation is a powerful tool for modeling the behav-
ior of large populations of neurons, enabling multiple applications (Depan-
nemaecker et al., 2023). Over the past decade, many mean-field approaches
based on different spiking models have been proposed (El Boustani & Des-
texhe, 2009; di Volo et al., 2019; Zerlaut et al., 2018; Carlu et al., 2020;
Chen & Campbell, 2022; Montbrió et al., 2015; Cakan & Obermayer, 2020;
Bandyopadhyay et al., 2021). One of these methods, based on a semiana-
lytical approach (El Boustani & Destexhe, 2009; Zerlaut et al., 2018), was
successfully applied to many different single-neuron models (di Volo et al.,
2019; Carlu et al., 2020), but never to a quadratic neuron model. In this work, we aimed to adapt and apply this method to a quadratic neuron model proposed by Izhikevich (2003).

Mean-field or neural-mass models are appropriate to model clinically recorded signals such as fMRI, EEG, or MEG (Tesler et al., 2022; Tesler, Linne et al., 2023) because they provide a simplified representation of the complex electrical and synaptic activity of large populations of neurons, helping to bridge scales (Depannemaecker et al., 2021). These models are based on the idea that the activity of a large group of neurons can be described by the average electrical activity of the group without having to consider the individual activity of each neuron.

One of the most widely used models in computational neuroscience is the quadratic integrate-and-fire neuron model (QIF). In this model, the membrane potential of a single neuron is described by a quadratic differential equation, and spikes are generated when the membrane potential reaches a given threshold. In an extension of this model, the adaptive-quadratic-integrate-and-fire neuron model (aQIF), a second slower variable describes the adaptive behavior and enables the system to capture a large repertoire of electrophysiological patterns (Izhikevich, 2003).

In this work, we build a model of cortical columns based on a balanced spiking network of aQIF neurons, composed of a population of excitatory regular-spiking (RS) neurons and a population of inhibitory fast-spiking (FS) neurons interacting through conductance-based synapses. This sparse network exhibits asynchronous irregular dynamics (Brunel, 2000) as observed in awake brain states. We build the corresponding mean field based on a previously developed master equation formalism (El Boustani & Destexhe, 2009; di Volo et al., 2019; Zerlaut et al., 2018), and we compare its dynamics to the spiking neural network. We show that the mean-field can correctly capture the dynamics of the network for both constant and time-varying inputs. In addition, we show that our semianalytical model remains valid for a wide range of cell parameters, which guarantees the robustness and generality of this formalism.

## 2 Methods

### 2.1 Spiking Network Model.
To build the spiking network model, we consider a system made by two neuronal populations: regular spiking excitatory cells (RS) and fast spiking inhibitory cells (FS). Each cell in the network is described by the quadratic point neuronal model proposed by Izhikevich (2003; equations 2.1 and 2.2):

$$ \frac{dv}{dt} = 0.04v^2 + 5v + 140 - u + I_{\text{syn}}, \qquad (2.1) $$

$$ \frac{du}{dt} = a(bv - u), \qquad (2.2) $$
where the first equation represents the membrane potential (being $I_{syn}$ the incoming synaptic currents) and the second equation represents a slow adaptation variable. These equations can be rewritten in a more general form:

$$ \tau_{Iz} \frac{dv}{dt} = g_{Iz}(v - E_{Iz})^2 - u + I_{syn}, \qquad (2.3) $$

$$ \tau_u \frac{du}{dt} = bv - u, \qquad (2.4) $$

When an action potential is emitted (i.e., the membrane potential crosses a threshold), the system is reset as in the equation 2.5 (and maintained at this state during a refractory period $t_{ref}$):

$$ \text{if } v \ge v_D \text{ then } \begin{cases} v \leftarrow c \\ u \leftarrow u + d \end{cases} \qquad (2.5) $$

With this formulation, $g_{Iz}$ could be thought of as a conductance, $E_{Iz}$ as a leak reversal potential, and $V_D$ and $c$ as the value of the membrane potential. Following the original model (Izhikevich, 2003), we keep these parameters dimensionless.

Neuronal interactions are mediated through synaptic inputs $I_{syn}$,

$$ I_{syn} = g_E(E_E - V) + g_I(E_I - V), \qquad (2.6) $$

where $g_{E,I}$ can be understood as the conductance of the excitatory and inhibitory synapses, respectively, and $E_{E,I}$ as the corresponding reversal potential. We model the conductances $g_{E,I}$ as a decaying exponential function that takes kicks of amount $Q_{E,I}$ at each presynaptic spike,

$$ g_{E,I} = Q_{E,I} \sum_{N_{pre}} \Theta(t - t_{sp}) e^{-\frac{t - t_{sp}}{\tau_e}}, \qquad (2.7) $$

where the sum goes through all presynaptic spikes, $\Theta$ is the Heaviside function, $\tau_e = \tau_i = 5$ is the decay timescale of excitatory and inhibitory synapses, and $Q_E = 1.5$ ($Q_I = 5$) is the excitatory (inhibitory) quantal increment (the change in conductance generated by a single spike). All the parameters used for simulations are indicated in Table 1. For the refractory period we used $t_{ref} = 5 \times 10^{-3}$. We assume that $\tau_u \gg \tau_{Iz}$, making the adaptation a slow variable with respect to the membrane potential.

**2.2 Mean-Field Model.** To build the mean-field model of our system, we follow the formalism proposed by El Boustani and Destexhe (2009; di Volo et al., 2019). This formalism provides a second-order mean field
<table><caption>Table 1: Parameters of the Mean-Field Model and the Network Model of a QIF Neurons.</caption><thead><tr><th colspan="3">Excitatory Neurons/Regular-Spiking Neurons</th></tr><tr><th>Parameter</th><th>Description</th><th>Value</th></tr></thead><tbody><tr><td>N<sub>e</sub></td><td>Number of excitatory neurons</td><td>8000</td></tr><tr><td>p</td><td>Probability of connection</td><td>0.05</td></tr><tr><td>K<sub>e</sub></td><td>Average number of excitatory synapses per neuron</td><td>N<sub>e</sub> × p</td></tr><tr><td>g<sub>Iz,e</sub></td><td>Leakage conductance of excitatory neurons</td><td>0.01</td></tr><tr><td>E<sub>Iz,e</sub></td><td>Reversal leakage potential of RS neurons</td><td>-65</td></tr><tr><td>τ<sub>Iz,e</sub></td><td>Time constant of membrane dynamics for RS neurons</td><td>1</td></tr><tr><td>T</td><td>Time constant of mean-field dynamics</td><td>5.10<sup>-3</sup></td></tr><tr><td>b<sub>RS</sub></td><td>Constant for voltage-related adaptation current for RS neurons</td><td>0</td></tr><tr><td>c<sub>RS</sub></td><td>After-spike reset voltage parameter for RS neurons</td><td>-65</td></tr><tr><td>d<sub>RS</sub></td><td>Spike-triggered adaptation parameter for RS neurons</td><td>1.5</td></tr><tr><td>τ<sub>u,e</sub></td><td>Time constant of adaptation variable for RS neurons</td><td>10</td></tr><tr><td>E<sub>e</sub></td><td>Reversal potential for excitatory synapses</td><td>0</td></tr><tr><td>Q<sub>e</sub></td><td>Quantal increment in excitatory conductance</td><td>1.5</td></tr><tr><td>τ<sub>e</sub></td><td>Decay time constant of the excitatory synaptic conductance</td><td>5.10<sup>-3</sup></td></tr></tbody><tbody><tr><th colspan="3">Inhibitory Neurons/Fast-Spiking Neurons</th></tr><tr><th>Parameter</th><th>Description</th><th>Value</th></tr><tr><td>N<sub>i</sub></td><td>Number of inhibitory neurons</td><td>2000</td></tr><tr><td>p</td><td>Probability of connection</td><td>0.05</td></tr><tr><td>K<sub>i</sub></td><td>Average number of inhibitory synapses per neuron</td><td>N<sub>i</sub> × p</td></tr><tr><td>g<sub>Iz,i</sub></td><td>Leakage conductance of FS neurons</td><td>0.04</td></tr><tr><td>E<sub>Iz,i</sub></td><td>Reversal leakage potential of the inhibitory neurons</td><td>-60</td></tr><tr><td>τ<sub>Iz,i</sub></td><td>Time constant of membrane dynamics for inhibitory neurons</td><td>1</td></tr><tr><td>T</td><td>Time constant of mean-field dynamics</td><td>5.10<sup>-3</sup></td></tr><tr><td>b<sub>FS</sub></td><td>Constant for voltage related adaptation current for inhibitory neurons</td><td>0</td></tr><tr><td>c<sub>FS</sub></td><td>After-spike reset voltage parameter for inhibitory neurons</td><td>-55</td></tr><tr><td>d<sub>FS</sub></td><td>Spike-triggered adaptation parameter for inhibitory neurons</td><td>0</td></tr><tr><td>E<sub>i</sub></td><td>Reversal potential for inhibitory synapses</td><td>-80</td></tr><tr><td>Q<sub>i</sub></td><td>Quantal increment in conductance</td><td>5.0</td></tr><tr><td>τ<sub>syn,i</sub></td><td>Decay time constant of the inhibitory synaptic conductance</td><td>5.10<sup>-3</sup></td></tr></tbody></table>

allowing us to derive a set of differential equations that describe the evolution of the mean firing rate $v_μ$ of each population, the covariance $c_{λη}$ between populations $λ/η$, and the average adaptation of the excitatory population $U$ (we assume that only excitatory RS population is affected by adaptation; see the parameters in Table 1), with $μ, λ, η = e, i$ (for excitatory and inhibitory, respectively):

$$T \frac{dv_μ}{dt} = (F_μ - v_μ) + \frac{1}{2} c_{λη} \frac{\partial^2 F_μ}{\partial v_λ \partial v_η}, \qquad (2.8)$$
$$T \frac{dc_{\lambda\eta}}{dt} = \delta_{\lambda\eta} \frac{F_{\lambda}(1/T - F_{\eta})}{N_{\lambda}} + (F_{\lambda} - v_{\lambda})(F_{\eta} - v_{\eta}) + \frac{\partial F_{\lambda}}{\partial v_{\mu}} c_{\eta\mu} + \frac{\partial F_{\eta}}{\partial v_{\mu}} c_{\lambda\mu} - 2c_{\lambda\eta}, \quad (2.9)$$

$$\frac{\partial U}{\partial t} = dv_e + (b\mu_V(v_e, v_i, U) - U)/\tau_u, \quad (2.10)$$

where the Einstein index notation (summation over repeated indices) is used and $\mu_V$ is the mean membrane potential, which we describe in further detail below. The parameter $T$, the time constant of the firing rate equations and covariance equations, is related to a main assumption used in this mean-field derivation: the network dynamics is considered to be Markovian within a time resolution $T$ (El Boustani & Destexhe, 2009). In this work, we set $T$ equal to the maximum neuronal firing rate $T = v_{\max}^{-1} = t_{ref}$, a common criterion used in previous work (El Boustani & Destexhe, 2009; Zerlaut et al., 2018).

The functions $F_e$ and $F_i$ correspond, respectively, to the transfer function of the excitatory and inhibitory neurons (i.e., each neural subtype's output firing rate when receiving excitatory and inhibitory inputs with rates $v_e$ and $v_i$). They are a function of the firing rates and of the adaptation: $F_{e,i}(v_e + v_{ext}, v_i, U)$, where $v_{ext}$ is the firing rate of an external drive, corresponding to the Poissonian external input in the spiking network model. Following Zerlaut et al. (2018), we can use the following general template function to write the transfer function for each neuronal type:

$$F_v = \frac{1}{2\tau_v} erfc \left( \frac{V_{thre}^{eff} - \mu_v}{\sqrt{2}\sigma_v} \right), \quad (2.11)$$

where erfc is the Gauss error function; $V_{thre}^{eff}$ is an effective neuronal threshold; and $\mu_V$, $\sigma_V$, and $\tau_V$ are the mean, standard deviation, and correlation decay time of the neuronal membrane potential.

To estimate the effective threshold, we write $V_{thre}^{eff}$ using a second-order polynomial expansion:

$$V_{thre}^{eff}(\mu_V, \sigma_V, \tau_V) = P_0 + \sum_{x \in \{\mu_V, \sigma_V, \tau_V\}^N} P_x \times \left( \frac{x - x^0}{\delta x^0} \right) + \sum_{x, y \in \{\mu_V, \sigma_V, \tau_V\}^2} P_{xy} \times \left( \frac{x - x^0}{\delta x^0} \right) \left( \frac{y - y^0}{\delta y^0} \right), \quad (2.12)$$

The coefficient of the expansion ($P_{xy}$) will be determined via a fit of equation 2.11 over the numerically obtained transfer function from simulation
of spiking cells for each neuronal type. The parameters $x^0$ and $\delta x^0$ (normalization constants) are defined according to the characteristic values of the associated variables for the range of inputs and parameters used. For these last constants, we use $\mu_V^0 = -45$, $\delta\mu_V^0 = 2$, $\sigma_V^0 = 4$, $\delta\sigma_V^0 = 5$, $\tau_V^0 = 0.005$, $\delta\tau_V^0 = 0.005$.

Considering asynchronous irregular regimes (Brunel, 2000), we make the assumption that the input spike trains follow Poissonian statistics (di Volo et al., 2019; El Boustani & Destexhe, 2009; Zerlaut et al., 2018) to calculate the averages ($\mu_{Ge,Gi}$) and standard deviations ($\sigma_{Ge,Gi}$) of the conductances, described in equation 2.13:

$$ 
\begin{aligned}
\mu_{Ge}(v_e, v_i) &= v_e K_e \tau_e Q_e, \\
\sigma_{Ge}(v_e, v_i) &= \sqrt{\frac{v_e K_e \tau_e}{2}} Q_e, \\
\mu_{Gi}(v_e, v_i) &= v_i K_i \tau_i Q_i, \\
\sigma_{Gi}(v_e, v_i) &= \sqrt{\frac{v_i K_i \tau_i}{2}} Q_i.
\end{aligned} 
\qquad (2.13) 
$$

In these equations, $K_e$ and $K_i$ are the average input connectivity received from the excitatory or inhibitory population respectively. As in the spiking network, $\tau_e = \tau_i = \tau_{syn}$ are synaptic time constants and $Q_e$ and $Q_i$ are the quantal increments of the conductances, respectively, for the excitatory or inhibitory population.

We can calculate the mean subthreshold membrane potential value by taking the stationary solution of equation 2.3. The quadratic form of equation 2.3 in $v$ gives rise to two solutions, of which only one is stable and used for the analysis, and the other is discarded. We thus obtain equation 2.14, which differs from the form obtained for other point neuron models (Zerlaut et al., 2018; di Volo et al., 2019; Carlu et al., 2020). Then, applying the approach described in previous work (Zerlaut et al., 2018), we determine $\sigma_v$ and $\tau_v$ as follows:

$$ \mu_V(v_e, v_i) = \frac{\left( (2g_{iz}E_{iz} + \mu_{Ge} + \mu_{Gi}) - \sqrt{(2g_{iz}E_{iz} + \mu_{Ge} + \mu_{Gi})^2 - 4g_{iz}(g_{iz}E_{iz}^2 + \mu_{Ge}E_e + \mu_{Gi}E_i - U)} \right)}{2g_{iz}}. \qquad (2.14) $$

Following Zerlaut et al. (2018), we can obtain $\sigma_v$ and $\tau_v$ by computing the power spectrum density of the membrane fluctuations. Within the assumption that in the asynchronous irregular state the spiking activity follows a Poissonian process, the power spectrum density of the fluctuations resulting from the sum of events $PSPs(t)$ (single postsynaptic potential) at
frequency $K_s v_s$ can be obtained from shot noise theory (Daley & Vere-Jones, 2008),

$$P_V(f) = \sum_{s \in \{e, i\}} K_s v_s \| \hat{P} \hat{S}_p(f) \|^2, \qquad (2.15)$$

where $PSP_s$ is the solution of equation 2.1 around the voltage $\mu_V$ for a single input spike of type s and $\hat{P}SP_s(f)$ is its Fourier transform. Then we can write (Zerlaut et al., 2018):

$$\sigma_v^2 = \int_R P_V(f), \qquad (2.16)$$

$$\tau_v = \frac{1}{2} \frac{\int_R P_V(f)}{P_V(0)}. \qquad (2.17)$$

Calculating $P_V$ from equation 2.15 and replacing in equations 2.16 and 2.17, we get

$$\sigma_V(v_e, v_i) = \sqrt{K_e v_e \left( \frac{2A_e B_e \tau_{s_e}^3}{\tau_{Iz_e}^2} + \frac{\tau_{s_e}^3 B_e^2}{8\tau_{Iz_e}} + \frac{B_e^2 \tau_{s_e}^3}{8\tau_{Iz_e}^2} \right) + K_i v_i \left( \frac{2A_i B_i \tau_{s_i}^3}{\tau_{Iz_i}^2} + \frac{\tau_{s_i}^3 B_i^2}{8\tau_{Iz_i}} + \frac{B_i^2 \tau_{s_i}^3}{8\tau_{Iz_i}^2} \right)} \qquad (2.18)$$

$$\tau_V = \frac{0.5 \left( K_e v_e \left( \frac{B_e^2 \tau_{s_e}^4}{2\pi \tau_{Iz_e}^2} \right) + K_i v_i \left( \frac{B_i^2 \tau_{s_i}^4}{2\pi \tau_{Iz_i}^2} \right) \right)}{K_e v_e \left( \frac{2A_e B_e \tau_{s_e}^3}{\tau_{Iz_e}^2} + \frac{\tau_{s_e}^3 B_e^2}{8\tau_{Iz_e}} + \frac{B_e^2 \tau_{s_e}^3}{8\tau_{Iz_e}^2} \right) + K_i v_i \left( \frac{2A_i B_i \tau_{s_i}^3}{\tau_{Iz_i}^2} + \frac{\tau_{s_i}^3 B_i^2}{8\tau_{Iz_i}} + \frac{B_i^2 \tau_{s_i}^3}{8\tau_{Iz_i}^2} \right)} \qquad (2.19)$$

where

$$\begin{align*}
A_e &= g_{iz}(\mu_{ve} - E_{iz})^2, \\
A_i &= g_{iz}(\mu_{vi} - E_{iz})^2, \\
B_e &= Q_e(E_e - \mu_{ve}), \\
B_i &= Q_i(E_i - \mu_{vi}). \tag{2.20}
\end{align*}$$

We point out that even if the model is dimensionless, we used units in the section 3 for the quantities that can be directly compared to experimental measurements, as Izhikevich (2003) did.

### 3 Results

In this section, we implement and validate the mean-field model for the adaptive quadratic integrate-and-fire model (Izhikevich, 2003) described in
<figure><img src="image_2.png" /><figcaption>Figure 1: Validation of the fitted semianalytical transfer function of the mean-field model. (a) Firing rates of a single aQIF cell with varying excitatory and inhibitory input. (b) Comparison of the firing rates of the fitted transfer function (dashed line) to the corresponding rates obtained from a single excitatory aQIF cell (dots) for varying external input. The transfer function is calculated using equation 2.11, where the corresponding stationary value of the adaptation *U* is used for each pair of excitatory and inhibitory inputs. (c) Comparison of the predicted mean-membrane voltage obtained from equation 2.14 (dashed lines) with the corresponding mean potential of a single excitatory aQIF neuron (dots) for varying external input. (d) Mean excitatory (black) and inhibitory (gray) conductances for the excitatory population in aQIF network simulations (dots) compared to the corresponding prediction from the mean field according to equation 2.13 (dashed lines) for varying external input.</figcaption></figure>

the previous sections. The first step to implement the mean-field formalism is estimating the semianalytical transfer function (TF) described in equations 2.12 and 2.11. The parameters of the TF are estimated by fitting the template transfer function with the output firing rate obtained numerically from a single Izhikevich neuron for varying inputs rates $v_e$ and $v_i$. The output of the single neuron is shown in Figure 1a. In Figure 1b, we show the comparison of the semianalytical TF obtained from the fit with the corresponding output rate of the single neuron. We see that the TF accurately captures the output rate obtained from the numerical simulations.
Two other relevant quantities predicted from the mean-field formalism are the membrane potential and the synaptic conductances. We display in Figures 1c and 1d the values obtained numerically for these quantities, together with the mean-field prediction as a function of the excitatory input rate, showing a good match between the numerical value and the one expected from the mean-field.

**3.1 Spontaneous Activity and Second-Order MF Evaluation.** Once the transfer function has been obtained and validated, we continue with the analysis of the mean field response and its comparison with the corresponding network simulations. We begin by analyzing the response of the mean field to constant external excitatory drives. We show in Figure 2a the results of the firing rates of the mean field together with the results obtained from the network as a function of the external excitatory drive. In addition, we show in Figures 2b and 2c the distribution of firing rates obtained from the network, together with the distribution predicted by the second-order mean field (equations 2.8 and 2.9).

As described in section 2, one of the key features of the Izhikevich model is the inclusion of the slow adaptation variable. Thus, capturing the effect of the adaptation on the system is of great relevance for the validation of the mean field. To this purpose, we show in Figure 2d the response of the mean field to a constant input as a function of the adaptation parameter *d*, and in Figures 2e and 2f the firing rate distribution together with mean-field prediction for different values of *d*. As we can see, the mean-field correctly estimates the impact of the adaptation on the firing rate of the network.

**3.2 Time-Varying Inputs.** In the previous section, we tested our formalism for a constant external input. We now turn to study the response of the mean field for a time-varying input. In particular, we test the mean-field response to stimuli of different amplitudes and speeds. We show in Figure 3 the results of our simulations together with the response of the network for a gaussian-shaped stimulus of various widths and amplitude. As we see in the figure, the mean field model can correctly capture the response of the network to the different inputs.

**3.3 Model Robustness and Parameter Exploration.** In principle, the estimation of the semianalytical transfer function (see equations 2.11 and 2.12) should remain valid under variation of the neuronal parameters, as these are explicitly taken into consideration within the formulation. This provides a large flexibility to our formalism and makes it analyzing different regimes emerging from the network suitable. To test the validity of the mean-field model under different parameters, we performed a parameter exploration and compared the results of the network with the predictions of the mean field. The results of this analysis are shown in Figure 4. We show the firing
<figure><img src="image_3.png" /><figcaption>Figure 2: Validation of the second-order mean-field model through comparison of the mean-field prediction and the corresponding result from the spiking-network. (a) Mean firing rates of the excitatory (green) and inhibitory (red) populations, with fixed adaptation *d* = 15 for different values of the constant external input *v*<sub>ext</sub>. Dots show the results obtained from simulation of the full network (receiving an external Poissonian input with the frequency *v*<sub>ext</sub>), with error bars of 1 standard deviation. The dashed lines show the average activity predicted from the mean field reduction, the shaded area spanning over ±1 predicted standard deviation. (b), (c) Probability distributions of the average excitatory (green) and inhibitory (red) population firing rates for two different constant inputs (indicated in panel a). The histograms show the results from full network simulation, obtained by binning a large time realization of the activity. The shaded line is the gaussian distribution with mean and standard deviation predicted by the corresponding mean-field reduction. (d) Mean firing rates of the excitatory and inhibitory populations, with fixed external input *v*<sub>ext</sub> = 10 Hz and for different values of the adaptation parameter *d*. Colors and style are the same as in panel a. (e)–(f) Probability distributions of the average excitatory and inhibitory population firing rates for two different adaptation values (indicated in panel d). Colors and style are the same as in panels b and c.</figcaption></figure>

rates obtained from the network and the mean field predictions for both
excitatory and inhibitory neurons, together with the adaptation variable *u*<sub>e</sub>.
In particular, we explored the parameters *E*<sub>Iz</sub> and *g*<sub>Iz</sub> that regulate the ex-
citability of the system. As shown in the figure, the mean field can capture
the behavior of the network for a wide range of parameters (one order of
magnitude in *g*<sub>Iz</sub>), although a discrepancy appears for low values of *E*<sub>Iz</sub> and
high *g*<sub>Iz</sub>.
<figure><img src="image_4.png" /><figcaption>Figure 3: Mean-field response to time-varying inputs. We show the response of the mean field (darker colors) for different time-varying inputs and the corresponding network simulation (lighter colors). The applied input is shown in dashed blue. We see that the mean-field can correctly capture the variation in the mean-firing rates driven by fast and slow inputs of different amplitude. (a) The input shown represents a fast and strong signal. (b) A fast but weaker signal. (c) A slow and strong signal. (d) A slow but weak signal.</figcaption></figure>

# 4 Discussion

In this letter, we have derived a mean-field model of populations of neurons described by the quadratic integrate-and-fire model, interacting with conductance-based synapses. We discuss the limitations and usefulness of this model relate it to previous work.

A major advantage of the mean-field approximation is that it allows large populations of neurons to be studied in a computationally efficient manner, while still capturing the important features of the dynamics. By capturing the average behavior of the population, it provides a computationally efficient and mathematically tractable way to study the dynamics of large-scale neural systems. It is worth noting that our analysis is limited to simple firing patterns, where neurons fire consistently in an asynchronous-irregular pattern in response to an external stimulus, in which case the validity of our formalism is guaranteed. A more in-depth study of the capacity of this formalism to capture other dynamical behavior of the spiking network, such as bursting, would be interesting for the future, and some initial analysis in this direction has been recently presented (Overwiening et al., 2023). As the model of a single neuron presents properties similar to those previously used (di Volo et al., 2019; Goldman et al., 2020, 2023; Stenroos et al.,
<figure><img src="image_5.png" /><figcaption>Figure 4: Robustness of the mean-field model. We analyzed the validity of the mean field for an extensive parameter exploration. In the figure, we show the stationary values of the firing rates (panels a–c for excitatory and d–f for inhibitory). The adaptation variable (panels g–i) resulted from the computation of the mean-field model (panels a, d, and g) and the spiking network (b, e, h) as a function of the parameters *E*<sub>1z</sub> and *g*<sub>1z</sub>, which regulate the excitability of the system. Panels (c, f, and i) show the difference between the mean-field calculations and the spiking network results. The original values used to build the transfer function are indicated with a pink square. We see that the mean field remains valid for a wide range of parameters, showing the robustness of the formalism.</figcaption></figure>

2023), we can expect that oscillatory dynamics such as slow waves must
exist. In this work, we focused on the basic dynamics for which this mean-
field approach was designed (El Boustani & Destexhe, 2009; asynchronous-
irregular regimes), and we have shown that this mean field can capture
this type of dynamics for different input types. Other approaches using the
same type of neuron model have been developed (Chen & Campbell, 2022).
These have the advantage of not being limited in temporal resolution by
the timescale of *T*. However, they do not account for network size effects or probability of connection (an all-to-all connection is assumed, which is far from biologically plausible).

It must be noted that a previous mean-field approach was proposed for the quadratic integrate-and-fire model (Montbrió et al., 2015), which was more recently extended to adaptive quadratic integrate-and-fire neurons (Chen & Campbell, 2022). To consider adaptation, the mean-field had to use a particular mathematical technique (called the Lorentzian ansatz) to enable the inclusion of adaptation in the mean field. In our approach, we directly integrate adaptation in the formalism and in the transfer function (di Volo et al., 2019), which is simpler. Another key difference is that these previous models work at the thermodynamic limit (when the number of neurons tends to infinity), while our approach is finite size and includes the network size in its parameters (El Boustani & Destexhe, 2009). Hence, our approach can also be seen as complementary to these previous works.

Finally, it is clear that the mean-field approximation is an important tool for studying large populations of neurons. By capturing the average behavior of the population, this approach provides a computationally efficient and mathematically tractable way to study the dynamics of large-scale neural systems up to entire brain regions (Tesler, Kozlov et al., 2023; Lorenzi et al., 2023) and the whole brain scale (Goldman et al., 2020, 2023). Future studies could characterize, as done for previous derivation (Kusch et al., 2023) and extend this approach to explore other possible dynamics, such as bursting, while also considering the limitations and advantages of other modeling approaches that use the same type of neuron model.

## 5 Conclusion

Modeling brain activity across different scales is a complex task, and there are various methods available to achieve this. Neural mass models and mean-field approaches have emerged as powerful tools to describe the dynamics of large populations of neurons and have been successfully applied to model clinically recorded signals such as fMRI, EEG, and MEG. In this work, we have adapted a mean-field method to the adaptive quadratic integrate-and-fire neuron model, a widely used model in computational neuroscience. We have built a model of cortical columns based on a balanced spiking network and compared its dynamics to the mean-field approximation. Our results demonstrate that the mean-field approach captures the asynchronous-irregular dynamics of the spiking neural network for different input types, making it a computationally efficient and mathematically tractable way to study the behavior of large-scale neural systems.

## Acknowledgments

This research resulted from a collaboration with students during the Fall School 2022 of the European Institute of Theoretical Neuroscience
(www.eitn.org). This work was funded by the European Union's Horizon 2020 Framework Program for Research and Innovation under Specific Grant Agreement 945539 (Human Brain Project SGA3) and the Centre National de la Recherche Scientifique (CNRS, France).

**Data Availability**

The code with the implementation of the spiking network and mean-field simulations is available in the following repository: https://github.com/sandrasaghir/Izhikevich-MF/tree/master.

**References**

Bandyopadhyay, A., Rabuffo, G., Calabrese, C., Gudibanda, K., Depannemaecker, D., Ivanov, A., ... Petkoski, S. (2021). Mean-field approximation of network of biophysical neurons driven by conductance-based ion exchange. 10.1101/2021.10.29.466427

Brunel, N. (2000). Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons. *Journal of Computational Neuroscience*, 8(3), 183–208.
10.1023/A:1008925309027

Cakan, C., & Obermayer, K. (2020). Biophysically grounded mean-field models of neural populations under electrical stimulation. *PLOS Computational Biology*, 16(4), e1007822. 10.1371/journal.pcbi.1007822

Carlu, M., Chehab, O., Porta, L. D., Depannemaecker, D., Héricé, C., Jedynak, ... di Volo, M., (2020). A mean-field approach to the dynamics of networks of complex neurons, from nonlinear integrate-and-fire to Hodgkin–Huxley models. *Journal of Neurophysiology*, 123(3), 1042–1051. 10.1152/jn.00399.2019

Chen, L., & Campbell, S. A. (2022). Exact mean-field models for spiking neural networks with adaptation. *Journal of Computational Neuroscience*, 50, 445–469.
10.1007/s10827-022-00825-9

Daley, D. J., & Vere-Jones, D. (2008). *An introduction to the theory of point processes vol. 2: General theory and structure*. Springer.

Depannemaecker, D., Destexhe, A., Jirsa, V., & Bernard, C. (2021). Modeling seizures: From single neurons to networks. *Seizure*, 90, 4–8. 10.1016/j.seizure.2021.06.015

Depannemaecker, D., Ezzati, A., Wang, H. E., Jirsa, V., & Bernard, C. (2023). From phenomenological to biophysical models of seizures. *Neurobiology of Disease*, 182, 106131. 10.1016/j.nbd.2023.106131

di Volo, M., Romagnoni, A., Capone, C., & Destexhe, A. (2019). Biologically realistic mean-field models of conductance-based networks of spiking neurons with adaptation. *Neural Computation*, 31(4), 653–680. 10.1162/neco_a_01173

El Boustani, S., & Destexhe, A. (2009). A master equation formalism for macroscopic modeling of asynchronous irregular activity states. *Neural Computation*, 21(1), 46–100. 10.1162/neco.2009.02-08-710

Goldman, J. S., Kusch, L., Aquilue, D., Yalçınkaya, B. H., Depannemaecker, D., Ancourt, K., ... Destexhe, A. (2023). A comprehensive neural simulation of
slow-wave sleep and highly responsive wakefulness dynamics. *Frontiers in Computational Neuroscience*, **16**, 2022. 10.3389/fncom.2022.1058957

Goldman, J. S., Kusch, L., Yalcinkaya, B. H., Depannemaecker, D., Nghiem, T.-A. E., Jirsa, V., & Destexhe, A. (2020). *Brain-scale emergence of slow-wave synchrony and highly responsive asynchronous states based on biologically realistic population models simulated in the virtual brain*. bioRxiv. 10.1101/2020.12.28.424574

Izhikevich, E. (2003). Simple model of spiking neurons. *IEEE Transactions on Neural Networks*, **14**(6), 1569–1572. 10.1109/TNN.2003.820440

Kusch, L., Depannemaecker, D., Destexhe, A., & Jirsa, V. (2023). *Dynamics and bifurcation structure of a mean-field model of adaptive exponential integrate-and-fire networks*. bioRxiv. 10.1101/2023.12.09.570909

Lorenzi, R. M., Geminiani, A., Zerlaut, Y., De Grazia, M., Destexhe, A., Gandini Wheeler-Kingshott, C. A., . . . D'Angelo, E. (2023). A multi-layer mean-field model of the cerebellum embedding microstructure and population-specific dynamics. *PLOS Computational Biology*, **19**(9), e1011434. 10.1371/journal.pcbi.1011434

Montbrió, E., Pazó, D., & Roxin, A. (2015). Macroscopic description for networks of spiking neurons. *Physical Review X*, **5**(2).

Overwiening, J., Tesler, F., Guarino, D., & Destexhe, A. (2023). *A multi-scale study of thalamic state-dependent responsiveness*. bioRxiv.

Stenroos, P., Guillemain, I., Tesler, F., Montigon, O., Collomb, N., Stupar, V., . . . Barbier, E. L. (2023). *How absence seizures impair sensory perception: Insights from awake FMRI and simulation studies in rats*. bioRxiv. 10.1101/2023.12.02.567941

Tesler, F., Kozlov, A., Grillner, S., & Destexhe, A. (2023). *A multiscale model of striatum microcircuit dynamics*. bioRxiv. 10.1101/2023.12.28.573546

Tesler, F., Linne, M.-L., & Destexhe, A. (2023). Modeling the relationship between neuronal activity and the BOLD signal: Contributions from astrocyte calcium dynamics. *Scientific Reports*, **13**(1). 10.1038/s41598-023-32618-0

Tesler, F., Tort-Colet, N., Depannemaecker, D., Carlu, M., & Destexhe, A. (2022). Mean-field based framework for forward modeling of LFP and MEG signals. *Frontiers in Computational Neuroscience*, **16**, 10.3389/fncom.2022.968278

Zerlaut, Y., Chemla, S., Chavane, F., & Destexhe, A. (2018). Modeling mesoscopic cortical dynamics using a mean-field model of conductance-based networks of adaptive exponential integrate-and-fire neurons. *Journal of Computational Neuroscience*, **44**(1), 45–61. 10.1007/s10827-017-0668-2