1. Neuroscience

Biophysically inspired mean-field model of neuronal populations driven by ion-exchange mechanisms

  1. Giovanni Rabuffo  Is a corresponding author
  2. Abhirup Bandyopadhyay
  3. Carmela Calabrese
  4. Kashyap Gudibanda
  5. Damien Depannemaecker
  6. Lavinia Mitiko Takarabe
  7. Sourin Chatterjee
  8. Maria Luisa Saggio
  9. Mathieu Desroches
  10. Anton Ivanov
  11. Marja-Leena Linne
  12. Christophe Bernard
  13. Spase Petkoski  Is a corresponding author
  14. Viktor Jirsa  Is a corresponding author
  1. Aix Marseille University, Institut de Neurosciences des Systèmes, UMR 1106, INSERM, France
  2. MathNeuro, INRIA, University of Montpellier, France
  3. MCENteam, Basque Center for Applied Mathematics (BCAM), Spain
  4. Faculty of Medicine and Health Technology, Tampere University, Finland
https://doi.org/10.7554/eLife.104249.3
Share this article
Cite this article
  1. Giovanni Rabuffo
  2. Abhirup Bandyopadhyay
  3. Carmela Calabrese
  4. Kashyap Gudibanda
  5. Damien Depannemaecker
  6. Lavinia Mitiko Takarabe
  7. Sourin Chatterjee
  8. Maria Luisa Saggio
  9. Mathieu Desroches
  10. Anton Ivanov
  11. Marja-Leena Linne
  12. Christophe Bernard
  13. Spase Petkoski
  14. Viktor Jirsa
(2026)
Biophysically inspired mean-field model of neuronal populations driven by ion-exchange mechanisms
eLife 14:RP104249.
https://doi.org/10.7554/eLife.104249.3
  2. Download BibTeX
  3. Download .RIS

eLife Assessment

This manuscript presents a useful mean-field model for a network of Hodgkin-Huxley neurons retaining the equations for ion exchange between the intracellular and extracellular space. The mean-field model derived in this work relies on approximations and heuristic arguments that, on the one hand, allow a closed-form derivation of the mean-field equations, but also raise questions about their justifications and the degree to which the results agree with experiments as well as direct numerical simulations. While the revised manuscript is much improved, reviewers continue to question the methodology for reducing model dimensionality and therefore the evidence for the utility of this approach remains incomplete at present.

https://doi.org/10.7554/eLife.104249.3.sa0

Abstract

Whole-brain simulations are a valuable tool for gaining insight into the multiscale processes that regulate brain activity. Due to the complexity of the brain, it is impractical to include all microscopic details in a simulation. Hence, researchers often simulate the brain as a network of coupled neural masses, each described by a mean-field model. These models capture the essential features of neuronal populations while approximating most biophysical details. However, it may be important to include certain parameters that significantly impact brain function. The concentration of ions in the extracellular space is one key factor to consider, as its fluctuations can be associated with healthy and pathological brain states. In this paper, we develop a new mean-field model of a population of Hodgkin–Huxley-type neurons, retaining a microscopic perspective on the ion-exchange mechanisms driving neuronal activity. This allows us to maintain biophysical interpretability while bridging the gap between micro- and macro-scale mechanisms. Our model is able to reproduce a wide range of activity patterns, also observed in large neural network simulations. Specifically, slow-changing ion concentrations modulate the fast neuroelectric activity, a feature of our model that we validated through in vitro experiments. By studying how changes in extracellular ionic conditions can affect whole-brain dynamics, this model serves as a foundation to measure biomarkers of pathological activity and provide potential therapeutic targets in cases of brain dysfunctions like epilepsy.

Introduction

Large-scale brain dynamics can be studied in silico with network models (Deco et al., 2011). Local activity can be represented by neural mass models (Deco et al., 2008), which coupled together through synapses, time delays, and noise (Deco and Jirsa, 2012; Petkoski and Jirsa, 2019; Perl et al., 2023; Forrester et al., 2024) allow the emergence of whole-brain activity that can be linked to empirical neuroimaging data (Sanz-Leon et al., 2015). In the context of large-scale simulations, these models have been employed for the study of resting-state brain activity (i.e., in the absence of any stimulus or task) in several mammalian species (Melozzi et al., 2019; Courtiol et al., 2020; Rabuffo et al., 2021), for the analysis and identification of chaos in brain signals (Baladron et al., 2012; Bandyopadhyay and Kar, 2018), the study of epileptic seizure genesis and propagation (Spiegler et al., 2011; Touboul et al., 2011; Nevado-Holgado et al., 2012; Hocepied et al., 2013; El Houssaini et al., 2015; Proix et al., 2017), among other applications.

At the mesoscopic level, the observable properties of a neuronal ensemble are generally explained by statistical physics formalism of mean-field theory (Wilson and Cowan, 1972; David and Friston, 2003; Moran et al., 2007; Wong and Wang, 2006). Mean-field models demonstrated a predictive value for studying the mesoscopic dynamics of neuronal populations (Izhikevich and Edelman, 2008), providing statistical descriptions of neuronal networks (Wilson and Cowan, 1972; Deco et al., 2008; Petkoski and Stefanovska, 2012; Montbrió et al., 2015; Montbrió and Pazó, 2020; Laing, 2017; Byrne et al., 2017; Devalle et al., 2018), which can be used to address questions related to network-level mechanisms (Bandyopadhyay and Kar, 2018; Petkoski and Stefanovska, 2012; Benitez Stulz et al., 2024). In general, neural mass models have a low enough number of parameters to be tractable and provide general intuitions regarding mechanisms underlying complex neuronal activity (Jirsa et al., 2014; Jirsa et al., 2017; Deco et al., 2021; Amunts et al., 2020; Depannemaecker et al., 2021; Carlu et al., 2020). For example, statistical population measures, such as the firing rate, can be used to assess mesoscopic dynamics (Shriki et al., 2003; Deco et al., 2011; Roxin et al., 2005; Luke et al., 2013; Jirsa et al., 2014; Sanz-Leon et al., 2015; Ashwin et al., 2016; Carlu et al., 2020; Hashemi et al., 2025).

Recently, a class of these models, called next-generation neural mass models (Byrne et al., 2020), has been developed based on an analytical approach introduced by Montbrió et al., 2015 that allowed for the exact derivation of mean-field parameters for a population of quadratic integrate-and-fire neurons. These can be linked to EEG/MEG oscillations (Byrne et al., 2022), including epipeltic seizures (Gerster et al., 2021), and have been used to study various aspects of the whole-brain dynamics such as the low-dimensional manifold of the resting state (Fousek et al., 2024; Gudibanda et al., 2025), aging (Lavanga et al., 2023) and neural signatures of consciousness (Breyton et al., 2024). Number of works dealt with the introduction of biologically realistic aspects in the mostly phenomenological neural mass model derived in Montbrió et al., 2015. These included short-term synaptic plasticity (Taher et al., 2020; Gast et al., 2021; Taher et al., 2022), spike frequency adaptation (Gast et al., 2020; Ferrara et al., 2023), spike timing-dependent plasticity (Duchet et al., 2023), synaptic delay (Devalle et al., 2018), random connectivity and noise (Goldobin, 2021; Klinshov and Kirillov, 2022; Goldobin et al., 2021), as well as an extension of the conductance-based neurons with a recovery variable (Guerreiro et al., 2022; Chen and Campbell, 2022; Gast et al., 2023).

Although it is practically inconvenient to reproduce the entire complexity of a neural mass, including all known biophysical parameters, it may be important to include parameters that can have widespread and general effects on neuronal activity (Amunts et al., 2020). The concentrations of Na+, K+, Ca2+, and Cl− ions in the extracellular space are key parameters to consider. Extracellular ion concentrations change dynamically in vivo as a function of the brain state, for example, between arousal and sleep (Ding et al., 2016; Rasmussen et al., 2017). Changes in extracellular ion concentrations can induce fluctuations in the spontaneous activity (Marchetti et al., 2005), in particular in the awake state (Rasmussen et al., 2019; Krishnan et al., 2018) as well as to switch one brain state to another (Ding et al., 2016). In particular, the extracellular potassium concentration [K+]ext plays a central role. Transient changes in [K+]ext can have large effects on cell excitability and spontaneous neuronal activity (Amzica et al., 2002; Ransom et al., 2000; Cressman et al., 2009), a result consistently reported in modeling studies (Bazhenov et al., 2004; Park and Durand, 2006; Fröhlich et al., 2008). Increases in [K+]ext are tightly controlled by astrocytes (Orkand et al., 1966), which can efficiently pump out [K+]ext and distribute it via their syncytium to prevent hyperexcitability (Breslin et al., 2018; Nwaobi et al., 2016; Crunelli et al., 2015; Hansson et al., 2000; Bedner and Steinhäuser, 2013). The saturation or lack of efficiency of these buffering mechanisms is often linked to a pathological state. Detailed single-neuron models demonstrate that continuous increases in [K+]ext can lead to different firing states, from tonic to bursting, to seizure-like events and depolarization block (Cressman et al., 2009; Depannemaecker et al., 2022b). In such detailed models, the buffering action of astrocytes is represented by a parameter named [K+]bath (Cressman et al., 2009; Breslin et al., 2018; Ullah et al., 2009; Nwaobi et al., 2016; Depannemaecker et al., 2022b). This parameter may also account for the potassium concentration in the perfusion bath registered in vitro. Our goal is to extend the use of ion concentration variables at the neural mass level, to be incorporated into whole-brain modeling.

In this study, we applied a mathematical formalism to estimate the mean-field behavior of a large neuronal ensemble, taking into account the ion exchange between the intracellular and extracellular space. Large networks of biophysically realistic neurons display complex behavior and – likely – do not satisfy the typical conditions required for deriving the mean-field dynamics (e.g., the Lorentzian Ansatz; Ott and Antonsen, 2008). In this work, relying on approximations and heuristic arguments, we derive a new neural mass model of a large population of an all-to-all coupled network of heterogeneous Hodgkin–Huxley (HH)-type neurons operating near synchrony (Figure 1). While the derivation is not exact, our model captures the mean-field behavior of connected neuronal populations operating in a synchronous regime, also displaying emergent dynamics not present in the single-neuron mode, as we show comparing our neural mass model outcome with the simulated activity from a large number of connected neurons. Our model faithfully characterizes the slow modulation of local field potential (LFP) fluctuations depending on ion concentrations, which we confirm through in vitro experiments. Considering different parameter configurations, we identify the mesoscopic states of the connected neuronal population, in various dynamical regimes that can be linked to different healthy and pathological states, and be of use in large-scale brain simulations. This approach establishes a link between the biophysical description at the cellular scale and the dynamics observable at the mesoscopic scale, enabling the study of the influence of changes in extracellular ionic conditions on whole-brain dynamics in health and disease.

Figure 1
Download asset Open asset
Biophysically inspired neural mass model.

Schematic diagram of the ion channel mechanism in extracellular and intracellular space in the brain. A biophysical model of a single neuron consists of three compartments (left panel): the intracellular space (ICS; in red), the extracellular space (ECS; in dark gray), and the external bath (EB; in light gray). The ion exchange across the cellular space occurs through the ion channels: Na+ gets inside the ICS (yellow channel), K+ gets out (green channel), the flow of Cl− can be bidirectional (purple channel); for the pump (blue), Na+ gets out and K+ gets into the ICS. A population of interacting neurons sharing the same [K+]bath concentration forms a local neural mass (middle panel), for which we model the mean-field equations in this work. Brain network model (right panel) with the activity of each brain region represented by neural masses.

Results

Biophysically inspired mean-field model

In this work, we derived a mean-field model describing the activity of a neural mass regulated by ion-exchange mechanisms at the cellular level (Figure 1). This model establishes a computationally accessible baseline that allows large-scale brain activity to be understood in terms of a few key biophysical details regulating the micro-scale mechanisms. The mean-field model equations (Equation 38) were derived by approximating a locally homogeneous network of HH-type neurons operating near synchrony, in the thermodynamic limit of an infinitely large population (see Methods section). By locally homogeneous, we mean that all neurons in the population are assumed to share the same extracellular and intracellular ionic environment and are connected with identical coupling rules, allowing us to treat the population as uniform with respect to ion dynamics and connectivity.

The single-neuron equations (Equation 1) were previously derived in Depannemaecker et al., 2022b as a simplified version of HH neurons which includes three compartments: an intracellular space, extracellular space, and an external bath in communication through ion fluxes (Figure 1, left). The decoupled single-neuron equations exhibit a range of activity patterns including bursting behavior, tonic spiking, seizure-like events, sustained ictal activity, and depolarization block (Figure 2a). As we will describe in more detail in the next sections, also the mean-field model exhibits these activity patterns, together with new complex behaviors emergent from the network interactions at the neural mass level. The variables describing a single neuron are characterized by a fast compartment, including membrane potential and gating variable fluctuations, and a slow compartment, which includes the slow-fluctuating ion concentrations (Figure 2b). The mean-field model inherits the single-neuron biophysical parameters (Table 1), which can be tuned to match different neuron types. For example, excitatory or inhibitory neurons can be characterized by tuning the reversal potential to high (e.g., Esyn=0 mV) and low (e.g., Esyn=−80 mV) values, respectively. It was previously established that a system of all-to-all coupled neuronal equations can be solved exactly in the thermodynamic limit (i.e., infinite neurons limit) if the single-neuron membrane potential equation is a quadratic function and if the instantaneous distribution of membrane potentials of neurons in a population is described by a Lorentzian (Montbrió et al., 2015). In our case, for any fixed value of the gating and potassium variables, the membrane potential equation resembles a cubic function (Figure 2c).

Figure 2 with 1 supplement see all
Download asset Open asset
Single Hodgkin–Huxley-type neuron model.

(a) Different patterns of electrophysiological activities previously identified in Depannemaecker et al., 2022b are also reproduced in our parameter setting by varying the potassium concentration in the external bath [K+]bath. The membrane potential V is measured in mV and the ion concentrations in mmol/m3. (b) Phase space trajectory of the seizure-like event simulation ([K+]bath=15.5). Fast oscillations occur in a fast subsystem identified by the membrane potential V and the gating variable n. The oscillations of the slow subsystem, here captured by the potassium concentration in the extracellular space [K+]ext, enable the transition to bursting. (c) Fixing the value of the state variables n,Δ[K+]int, and [K+]g as constants, the membrane potential equation resembles a cubic function for different values of [K+]bath. We can model this function as a step-wise quadratic approximation, corresponding to two parabolas with vertices at coordinates (c−,I−) and (c+,I+) and curvature R− and R+, respectively (d). The two parabolas meet at an intersection point V⋆ where the membrane potential equation changes curvature. (e) At each time, we assume that the membrane potential of a neuronal population is distributed according to a Lorentzian centered at y=y(η,t) and with width x=x(η,t), for each value of the excitability η (Lorentzian Ansatz). In the case depicted, the cubic function meets the zero for V<V⋆, the neuronal population is described by the Lorentzian distribution in blue in the steady-state solution, and the neuronal dynamics is governed by the positive parabola according to the continuity equation. In the case where the derivative of the membrane potential crosses zero for V>V⋆ (e.g., if the cubic function is shifted up by adding a constant current to the membrane potential derivative), the population is described by the red distribution in the steady state, and the continuity equation is governed by the negative parabola equation. Cases where the cubic function meets the zero in more than one point are not well described by this approximation, see Steady-state solution and Lorentzian Ansatz section.

Table 1
List of parameters and their values used for the single-neuron simulation.
ParameterSymbolValue
Membrane capacitanceCm1nF
Gating time constantτn4ms
Chloride conductancegCl7.5 nS
Maximal potassium conductancegK22 nS
Maximal sodium conductancegNa40 nS
Potassium leak conductancegK,l0.12 nS
Sodium leak conductancegNa,l0.02 nS
Intracellular volumeωi2160 μm3
Extracellular volumeωext720 μm3
Intra-/extracellular volume ratioβ=ωi/ωext3
Conversion factorγ0.04 mol/C
Diffusion rateε0.001 mHz
Maximal Na/K pump currentρ250 pA
Initial concentration of extracellular K[K+]0,ext4.8 mmol/m3
Initial concentration of intracellular K[K+]0,int130 mmol/m3
Initial concentration of extracellular Na[Na+]0,ext138 mmol/m3
Initial concentration of intracellular Na[Na+]0,int16 mmol/m3
Initial concentration of extracellular Cl[Cl−]0,ext112 mmol/m3
Initial concentration of intracellular Cl[Cl−]0,int5 mmol/m3

Therefore, to proceed with the mean-field reduction, we make several key assumptions: (1) The cubic-like profile can be described by a step-wise quadratic function corresponding to two parabolas with opposite curvature (Figure 2d). (2) The membrane potential distribution of HH-type neurons is described at each time by a Lorentzian (Figure 2e), which is a good approximation in situations where neurons are operating close to synchrony (e.g., Figure 3a), but is not a suitable description for other dynamical states better described by a bimodal distribution (e.g., Figure 1a, orange distribution). (3) We assume that the potassium concentrations, both intracellular (Δ[K+]int) and extracellular (through the buffering variable [K+]g), are homogeneous across the neuronal population. This is justified physiologically by the rapid redistribution of ions through diffusion and electrochemical gradients, which enforce near-instantaneous equilibration at the mesoscopic scale. As such, assigning separate compartments to each neuron is neither practical nor biologically meaningful in this context. (4) We assume that the gating variable n, which governs potassium conductance, can be treated as a population-averaged variable. This allows us to describe the neuronal ensemble using a reduced set of collective (mean-field) variables.

Figure 3 with 1 supplement see all
Download asset Open asset
Mean-field model versus neuronal network across dynamical regimes.

(a) Example raster plot of a population of N=3000 all-to-all coupled HH-type neurons displaying sub- and supra-threshold dynamics that can be well described by a Lorentzian distribution. (b) For several [K+]bath values we simulated the activity of a population of N=3000 all-to-all coupled HH-type neurons across dynamical regimes (parameters in Supplementary file 1). We compare the mean membrane potential Vpop and external potassium [K+]ext of such population (in blue) with the results obtained using the mean-field model equations (in green). To properly match the slow timescale of the population, we defined an effective value of the potassium concentration in the bath [K+]batheff, represented in panel (c) (the dashed line represents the identity).

Using these key approximations, we derived the primary result of this work: a closed form for the neural mass equation (Equation 38) expressed in terms of mean-field variables.

Comparison with neural network simulations

The accuracy of the mean-field approximation is first validated by comparing it to the simulation of a large network of coupled HH-type neurons. The simulation of HH-type single-neuron dynamics (Equation 1) driven by an ion-exchange mechanism, has revealed that the parameter [K+]bath has a significant impact on the dynamics (Figure 2a). Thus, for different [K+]bath we simulated a network of N=3000 coupled neurons described by Equation 1, coupled as in Equation 9 with synaptic strength J=1 and heterogeneous excitability distributed with mean η¯=0 and width Δ=1.

Using the mean-field equations (Equation 38), we show that the mean-field model membrane potential VMF matches qualitatively the average membrane potential Vpop of the neuronal population (Figure 3b). We stress that VMF was obtained by solving 5 coupled equations Equation 38, while Vpop was obtained by averaging the solution of 4 × 3000 equations (Equation 1) coupled as in Equation 9. This result shows that the mean-field model can be used to simulate several regimes of activity and emulate the average dynamics of a large neuronal population in regimes where the membrane potentials’ distribution is unimodal and can be reasonably approximated by a Lorentzian. These regimes include activities such as spike trains, tonic spiking, bursting, seizure-like events, status epilepticus-like events, and depolarization block (such as in Figure 3b). Notice that in such regimes the population average displays dynamics qualitatively similar to the single-neuron simulation (Figure 2b). However, due to the network effects in the population, there is a frequency shift compared to the single-neuron dynamics (for example, the bursts in the [K+]bath=7.5 regime are ∼3 times faster in the population (Figure 3b) than in the single neuron (Figure 2a)).

A frequency shift is also present in the mean-field model compared to the population dynamics. However, the definition of an effective [K+]batheff can adjust this frequency shift (Figure 3c). In Figure 3c, we report the [K+]batheff values obtained via visual inspection of the match with the population dynamics for two values of synaptic strength J=0.1 and J=1, and used in (Figure 3b). Also note that the gating variable n is treated as microscopic in the neural network, while in the derivations for the mean field it is considered as a mesoscopic and identical for the whole population. This is likely responsible for some of the discrepancies between the two modalities (Taher et al., 2020; Taher et al., 2022; Gast et al., 2021).

Comparison with in vitro experiments

We rearranged Figure 4 to better highlight the range of dynamics involving extracellular potassium concentration and bursting activity. The neural mass model (Figure 4a) reproduces a periodic behavior of the extracellular potassium concentration, with fast voltage bursts riding on top (we used the parameters Cm=16,τn=8,ϵ=0.0001,γ=0.00025). A similar pattern is observed in the in vitro recordings (Figure 4c), although we emphasize that these are AC-coupled LFP traces. As such, slower components of the membrane potential – including jumps during bursting – are filtered out and not visible in the LFP recordings. We do not attempt to classify the bursting pattern in either data or simulations. To compare with experimental timescales, we also simulate a network of N=3000 HH-type neurons (Figure 4b). The parameters were adjusted for computational efficiency, yielding a shorter bursting period, but preserving the key qualitative feature: modulation of fast activity by slow potassium dynamics. More complex patterns can emerge in vitro as well. For example, in Figure 4d, we observe a progressive slowing of the burst frequency, which our deterministic model does not capture. We hypothesize that such behavior may arise from slow parameter drifts or noise-driven transitions between metastable regimes – effects that are not included in the present model, but could be explored by introducing noise or varying parameters such as [K+]bath or the nullcline geometry (see e.g., Saggio et al., 2020). Finally, Figure 4e shows an emergent regime from the mean-field model with isolated bursts in the up state, qualitatively resembling some features of the in vitro activity.

Figure 4
Download asset Open asset
Comparison of numerical results and in vitro experiments.

(a) Neural mass model showing slow periodic fluctuations in extracellular potassium concentration (green) with voltage bursts (blue) riding on top. (b) Network simulation of N=3000 coupled HH-type neurons exhibiting a similar bursting pattern at a shorter timescale (due to rescaled parameters). (c) In vitro recording showing LFP; blue, AC-coupled and extracellular potassium concentration (green), exhibiting slow periodic fluctuations and bursting activity. (d) In vitro trace showing more complex dynamics – after a shift to a high-potassium state (marked by arrows), the burst frequency slows down progressively. (e) Simulation from the mean-field model showing an emergent bursting regime with isolated events in the up state, not seen at the single-neuron level.

Bifurcation analysis: emergent network states and multistability

A previous study has established that varying the external potassium concentration ([K+]bath) induces significant changes in the dynamical behavior of the single-neuron model (Depannemaecker et al., 2022b). Within a specific range of [K+]bath, the system exhibits multistability, a phenomenon crucially analyzed using bifurcation analysis. This approach allows us to assess system stability and qualitative behavior independently of initial conditions.

We here applied multiple-timescale bifurcation analysis to the neural mass model (Equation 38). Indeed, we treated the variables Δ[K+]int and [K+]g as slow compared to the other variables (due to the relative magnitude of parameters ε, γ, and ωi) and we considered the so-called fast subsystem. This corresponds to the limiting system in which slow variables’ dynamics are frozen and slow variables become parameters that are varied to discover the corresponding bifurcation structure (see Figure 5a). This analysis unveiled complex regimes within the neural mass model, including distinct regions of multistability. The oscillatory behavior observed in the fast subsystem occurs within the range bounded by the Saddle-Homoclinic (SH), the Hopf2, and the Fold of Limit Cycles FLC1 curves (gray area).

Figure 5
Download asset Open asset
Fast subsystem bifurcation diagram in two parameters.

(a) The bifurcation diagram shows the behavior of the fast subsystem of system (Equation 38) to the slow variables Δ[K+]int and [K+]g, which are parameters in the fast-subsystem limit. Each curve represents parameter values for which specific bifurcations occur, obtained through a numerical continuation algorithm. Saddle Node (SN) bifurcations are shown in blue, Saddle Homoclinic (SH) bifurcations in black, Fold Limit Cycle (FLC) bifurcations in green, and Hopf bifurcations in red. Colored dots indicate codimension-2 points. In particular, a Bogdanov–Takens (BT) point and a Saddle-Node-Loop (SNL) point also identified in the single-neuron model, and epileptor model. Note that this diagram is quite involved and we do not claim that the present version is complete, however it is representative of the complexity of the bifurcation structure of the mean-field model. (b1, b2) Comparison of the fast subsystem bifurcation diagrams for the single-neuron model and the neural mass model. The left side shows a zoomed-in version of the diagram from panel (a) to facilitate comparison with the single-neuron model diagram on the right. Both diagrams use the same color coding as described above. This comparison highlights the similarities and differences in the bifurcation structures of these two models and indicates emergent structures due to interactions within the network.

The qualitative dynamics, illustrated in Figure 3b, are influenced by temporal variations in the slow variables and the specific bifurcations encountered in the fast subsystem diagram. They occur in the surroundings of the Saddle-Node-Loop (SNL) point (marked with a yellow star in Figure 5b), which is the region of local topological equivalence between the single neuron and the neural mass model. For instance, in our previous study on single-neuron dynamics (Depannemaecker et al., 2022b), we have shown that seizure-like events arise from transitions between fixed point and limit cycle dynamics and vice versa through SN and SH bifurcations (El Houssaini et al., 2015). Similarly, in the neural mass model, seizure-like events (Figure 3b at [K+]bath=15.5) also result from transitions involving these bifurcations. Throughout the burst phase, the slow variables remain constrained between the SN2 and SH curves.

However, further away from the SNL point, the topological equivalence between single neuron and neural mass models is broken. In summary, our neural mass model displays a richer bifurcation structure than the single neuron one (Figure 5b, right versus left panels), with new complex activity regimes, which is also supported by the simulation of new behaviors as in Figure 4e.

While these results confirm that the neural mass model is not simply mirroring the single neuron one, a detailed numerical exploration of the activity regimes for a population of all-to-all coupled HH-like neurons is required to confirm the accuracy of the neural mass model in representing population behavior emerging from network interactions.

Coupled neural masses

Finally, we show that the presented mean-field model can be used to perform large-scale network simulations of brain activity, for example, to be used in the context of brain stimulation or epilepsy. To do so, we coupled six neural masses via long-range structural connections with random weights (Figure 6a). Each population (network node) is described by the mean-field model derived in this article (Equation 38). Nodes A, B, C, D, E, and F are tuned to a healthy regime with [K+]bath=5.5 which do not show pathological bursts when decoupled from other external inputs. Node D is tuned to a pathological regime, with [K+]bath=15.5. We run a simulation of the system by setting the global coupling parameter to G=0 (Figure 6b) and G=100 (Figure 6c). In the first case, the system is effectively decoupled and, as expected, the only node displaying pathological activity is node D. However, in the second case, when G is increased, the system is coupled and the pathological fluctuations of node D differentially affect the activity of all other nodes at the whole-network level. The evoked activity in downstream regions is induced by a network re-organization phenomenon, as the local parameters were not modified in these regions.

Figure 6
Download asset Open asset
Network simulation of structurally connected neural mass models and propagation of pathological bursts.

(a) Structural connectivity for six all-to-all connected nodes A, B, C, D, E, and F with random weight allocation. Each node is described by a neural mass model derived as the mean-field approximation of a large population of Hodgkin–Huxley (HH)-type neurons. (b) When decoupled (global coupling G=0), all the nodes operate in a ‘healthy’ regime with the potassium concentration in their bath set to low values [K+]bath=5.5, except for node D which is tuned into a pathological regime [K+]bath=15.5 characterized by the spontaneous presence of bursts. (c) When the global coupling is sufficiently increased, the pathological value of [K+]bath in node D generates bursts that diffuse through the connectome mimicking the spreading of a seizure.

This simulation serves as a proof of concept to illustrate how local pathological activity can propagate through a network depending on the strength of coupling. We used a single representative realization of randomly weighted structural connectivity. While we did not perform a systematic exploration of different realizations or coupling strengths, we observed that the qualitative behavior – namely, the emergence of network-wide bursting beyond a critical coupling threshold – remains robust across similar setups. The model is compatible with empirical connectome data and can be readily extended to simulations using realistic brain network architectures.

Limitations of the model

The mean-field model derived in this work relies on approximations and heuristic arguments that, on the one hand, allow a closed-form derivation of the mean-field equations (Equation 38), and on the other hand restrict its validity to a limited regime of activity corresponding to quasi-synchronous neuronal populations. Therefore, rather than an exact mean-field representation, the model provides a qualitatively realistic description of a mesoscopic population of connected neurons driven by ion-exchange dynamics. Moreover, the discrepancy between the two modalities would have likely been smaller if for the neural network we also adopted a gating variable that is mesoscopic and identical across the spiking neurons, as in similar works (Taher et al., 2020; Taher et al., 2022; Gast et al., 2021). However, here we demonstrate the validity of the mean-field approximation even for the more natural, microscopic representation of the gating variable in the neural network.

The approximation of the membrane potential dynamics as a step-wise quadratic function and the assumption of Lorentzian distributed membrane potential in the population allowed us to apply the Ott-Antonsen Ansatz and to develop the mean-field approximation of a heterogeneous network of biophysical neurons driven by ion-exchange dynamics coupled all-to-all via conductance-based coupling. With these approximations, the mean-field equations do not show a one-to-one correspondence with the neural network simulations, except in regimes where the distribution of the membrane potentials of the neurons is well described by a Lorentzian (Figure 3a). In fact, it is not guaranteed that such a population of HH-type neurons allows for a closed formulation of the thermodynamic limit. For example, the distribution of the membrane potentials resembles a Lorentzian only in some cases (e.g., when all the neurons are unimodally distributed below or above the threshold; Figure 4b, bottom left), while it can deviate from it in other moments (e.g., when a subpopulation of neurons transition suprathreshold while another population remains subthreshold; Figure 4b, bottom right).

The mean-field model could be improved, for example by allowing for a bimodal distribution of the membrane potential (e.g., a double-Lorentzian approximation), although this condition is not guaranteed to allow a closed-form derivation of the mean-field equations. Furthermore, the parabola coefficients c−,c+,R−,R+ were fixed as constants, however, these coefficients could be made functions of the slow variables and the gating variable, which might unveil new dynamical regimes and extend the validity of the thermodynamic limit beyond the regimes described in this work. Also, in the case of constant values, an in-depth exploration of the parameter space is required to fully characterize the model and its bifurcation structure. Other limiting assumptions are the moment closure condition (Equation 19) and the assumptions that the functions (Equation 3) averaged across the neuronal population can be expressed as functions of the average membrane potential V and gating variable n (which is only true in the cases where the functions (Equation 3) can be reasonably approximated as linear functions in a range of V and n).

The definition of the firing rate in the thermodynamic limit was approximated assuming a firing threshold for membrane potential going to infinity as in Equation 25, however, more refined definitions (e.g., accounting for a heterogeneous firing threshold) can be considered (Gast et al., 2022). In this work, we showed that the model can retrieve emergent network behaviors similar to those observed in neuronal network simulations and in vitro. This correspondence is not quantitative because the network size is finite, and, in experimental observations (in vivo/in vitro networks), the connectivity profile is not all-to-all connected. Despite its simplicity, the model displays desired properties such as bistability in healthy and pathological regimes (Figure 5) and sensitivity to external coupling (Figure 6). This approach, taking into account key biophysical details, offers a first step in considering the role of the glia in neural tissue excitability. Following this direction, other ions, such as calcium should be taken into consideration, as well as other effects such as plastic synapses, random connectivity, noise, adaptation, spike-timing-dependent plasticity, as already discussed in the Introduction.

Discussion

A quest of modern neuroscience is understanding and explaining how the brain operates for healthy individuals, and how it deviates from its healthy state in case of different diseases and disorders. Most brain imaging techniques register the collective electrical and metabolic fluctuations of large neuronal populations. Because these fluctuations emerge from the electrochemical interactions of many neural cells, they cannot be fully understood at the microscopic level (e.g., as the mechanism of a single neuron). Instead, they must be studied at a mesoscopic level, where the emergent behavior of the entire population can be observed. Also, besides the recent progress (Markram et al., 2015), it is still not possible to compute the behavior of a few billion such neurons (which comprises a mammalian brain), and even if it was, it remains debatable what knowledge we would gain from such a complex system (Frégnac, 2017). Hence, an appropriate way to model recorded brain activities is to write mathematical descriptions for large neuronal populations (Sanz-Leon et al., 2015), using formalism from statistical physics and mathematics. In this case, one can measure population properties such as the mean membrane potential, rather than the activity of an individual neuron. Some of these observables acquire meaning only in the context of mean-field analyses, for example, temperature or pressure in the case of molecular ensembles, order parameter (Kuramoto, 1984) and mean ensemble frequency (Petkoski et al., 2013) for coupled oscillators, or more specifically the firing rate in the case of neuronal population (Gerstner and Kistler, 2002). However, to the best of our knowledge, no theory to date can explain the behavior of a large population of neurons (brain areas) from the perspective of the driving mechanism of the neuronal activities, that is the ion-exchange and transportation dynamics at the cellular level. Different pathological trials and experiments have already revealed that changes in ion concentration in brain regions could result in different brain dysfunctions but the pathway of these phenomena is still unclear.

In this study, we developed a biophysically inspired mean-field model for a network of all-to-all connected, heterogeneous HH-type neurons, which are characterized by ion-exchange mechanism across the cellular space. Unlike phenomenological or reduced models, the HH framework allows us to retain explicit ion-exchange dynamics, which are essential for linking membrane behavior to extracellular potassium fluctuations. This level of biophysical detail is crucial for modeling pathological regimes such as seizure onset and propagation.

The intermediate approximation into the step-wise quadratic function allowed us to apply the analytical formalism to the complex HH-type equations. The further assumption that the membrane potentials of neurons in a large population are distributed according to a Lorentzian (Lorentzian ansatz) makes the mean-field approximation analytically tractable. From the simulation of the network behavior of such neurons, it is evident that the mean-field model captures the dynamic behavior of the network when the population is highly synchronized. This assumption might find applications in the study of seizure dynamics. Future work is needed to explore the non-synchronous cases.

The bifurcation analysis reveals that the mean-field model is capable of qualitatively capturing emergent neuronal network states observed in numerical simulations, as well as in vitro (e.g., Figure 4). The significance of such qualitative analysis lies in its relevance to a recently proposed classification of seizure dynamotypes, which categorizes seizures based on their observable characteristics and dynamic composition (Saggio et al., 2020). This classification was supported by a model-based bifurcation analysis using the ‘Epileptor’ model as a neural mass model (Jirsa et al., 2014). Our model complements the Epileptor by allowing the interpretation of parameters in terms of measurable microscopic quantities, specifically ion concentrations, rather than phenomenological parameters. Therefore, our model allows us to translate the classification of seizure dynamics from electrophysiology recording with mesoscopic scale resolution (S/M/EEG) to experiments with microscopic details, as we have shown in-vitro. While our study contributes to both the qualitative agreement between our model and in-vitro seizure patterns and the improved interpretability of phenomenological models, a full bifurcation analysis, beyond the one presented in this work is required to align our results with the dynamical patterns of the Epileptor model.

The derived mean-field model relates the slow-timescale biophysical mechanism of ion exchange and transportation in the brain to the fast-timescale electrical activities of large neuronal ensembles. In fact, as demonstrated via brain imaging of different modalities (Rasmussen et al., 2020), ionic regulation acts on a relatively slow time scale on the fast electrical activity of neurons.

Analyzing the model, we found the potassium concentration in the external bath (parameter [K+]bath) to be a significant determinant of the dynamics, monitoring the biophysical state of the neuronal population (i.e., a brain region). For low ion concentrations in the external bath, the mean-field model demonstrates the existence of healthy brain dynamics. Increasing the excitability in terms of the ion concentration leads to the appearance of spike trains, tonic spiking, bursting, seizure-like events, status epilepticus-like events, and depolarization block, similar to the Epileptor model (Jirsa et al., 2014; El Houssaini et al., 2015; El Houssaini et al., 2020).

For several values of the potassium concentration in the external bath, or given an external input current, we showed that our model can be tuned into a bistable regime characterized by the coexistence of high and low firing rate states, which is the hallmark of many mean-field representations of linear (Di Volo et al., 2019; Zerlaut et al., 2018) or quadratic integrate and fire neurons (Coombes and Byrne, 2019; Montbrió et al., 2015), and rate models (Wong and Wang, 2006). Bistability is a desired feature for mean-field models. For example, bistability might provide a mechanism for the dynamic occurrence of so-called up (high firing) and down (low firing) states, as observed both in vitro (Plenz and Kitai, 1998; Cossart et al., 2003) and in vivo under several conditions such as quiet waking, anesthesia, slow-wave sleep and during perceptual task across several species (Steriade et al., 1993; Luczak et al., 2007; Vyazovskiy et al., 2011; Engel et al., 2016; Jercog et al., 2017). In fact, when bistable neural masses are coupled through a connectome and driven by a fine-tuned stochastic input noise, the simulated activities in brain regions can spontaneously jump between high and low firing rate states, which provides a mechanism for dynamic Functional Connectivity (Rabuffo et al., 2021), as observed in large-scale brain recordings (Preti et al., 2017).

Thus, the derived mean-field model links the presence of high and low firing rate states as well as the spiking and bursting neuroelectric behaviors with the biophysical state of the neuronal ensemble (brain regions) in terms of the ion concentrations across the cellular space. The effect of constant stimulus current is also analyzed and it is observed that even within the healthy regime, several stimulations, which could either correspond to external stimuli or inputs from some other brain regions, could generate transient spiking and bursting activity (Figure 5d). This result is particularly interesting in the case of brain stimulation. For example, several kinds of brain stimulation protocols in epileptic patients have already been found to generate epileptic seizures in pathological practice (Fisher and Velasco, 2014; Kahane and Depaulis, 2010). Our results demonstrate that these phenomena could be reproduced and tracked in a biophysical-inspired modeling approach. Therefore, we assume that the derived model could potentially be applied to improve predictive capacities in several types of brain disorders, particularly epilepsy.

Using conductance-based coupling between six neuronal masses we also demonstrated that a hyper-excited population of neurons can propagate bursting and spiking behavior to otherwise healthy populations. This could lead the path to understanding how brain signals propagate as a coordinated phenomenon depending on the distribution of biophysical quantities and structural as well as architectural heterogeneity on the complex network structure of the connectome. Until now, mean-field models used in large-scale network simulations, like the Virtual Brain (Sanz-Leon et al., 2015), did not take into consideration the extracellular space. Our mean-field model is a first step toward the integration of biophysical processes that may play a key role in controlling network behavior, as shown at the spiking network level (Depannemaecker et al., 2022a).

In this work, we focused on [K+], given its known role in neuronal activity. Changes in K+ occur when the brain alternates between arousal and sleep (Ding et al., 2016). Although a causal relationship is not clearly established, seizures and spreading depression, which is assumed to underlie certain forms of migraine (Tottene et al., 2009; Vinogradova, 2018), are associated with large (>6 mM) and very large (>12 mM) concentrations of K+ (Tottene et al., 2009; Hertz and Chen, 2016). The invariant increase of K+ concentration during seizure (Heinemann et al., 1986), provokes a saturation of potassium buffering by the astrocytes. In fact, the extracellular concentration of K+ is tightly controlled by astrocytes (Breslin et al., 2018; Nwaobi et al., 2016; Hansson et al., 2000). Most large-scale simulations do not integrate glial cells, which make up half of the brain cells (Hansson et al., 2000; Breslin et al., 2018). Their functions are altered in most, if not all, brain disorders, in particular, epilepsy (Bedner and Steinhäuser, 2013; Crunelli et al., 2015). Our approach indirectly accounts for the astrocytic control of extracellular K+. Other ion species also vary during arousal/sleep and seizures, in particular, Ca2+ (Ding et al., 2016; Pocock and Kettenmann, 2007; Auld and Robitaille, 2003; Fernández-Chacón et al., 2001). A decrease in Ca2+ will decrease neurotransmission and thus change cell-to-cell communication (Pocock and Kettenmann, 2007; Auld and Robitaille, 2003; Fernández-Chacón et al., 2001). Future studies are needed to integrate Ca2+ in neural mass models. This paper serves as an introductory exploration of our novel approach, aiming to establish its feasibility and potential applications for modeling large ensembles of biophysically realistic neurons. A comprehensive and systematic investigation to address the limitations discussed in section ‘Limitations of the model’ and Figure 2—figure supplement 1 will be undertaken in future studies.

To date, the Epileptor model (Jirsa et al., 2014) is a gold standard for inferring the location of the epileptogenic zone, with direct applications in clinics (Wang et al., 2023; Jirsa et al., 2023). Nonetheless, the Epileptor model remains phenomenological, and parameters (e.g., epileptogenicity) are non-interpretable in terms of measurable quantities. By addressing the biophysical information at the neuronal level, our mean-field formalism allows keeping biophysical interpretability while bridging between micro- to macro-scale mechanisms. The mean-field model derived in this work aggregates a large class of brain activities and repertoire of patterns into a single neural mass model, with direct correspondence to biophysically relevant parameters. This approach could serve as a computational baseline to address core questions in epilepsy research. In particular, how to identify the multiscale mechanisms implicated in epileptogenicity and propagation of seizures. This would eventually lead to establishing pathologically measurable bio-markers for large-scale brain activities and consequently offering therapeutic targets for different brain dysfunctions.

Materials and methods

Immature hippocampus in vitro preparation

Request a detailed protocol

In accordance with French law, ex vivo experiments do not require approval from an ethics committee. All procedures performed on live animals (deep sedation and euthanasia by decapitation prior to brain extraction) were carried out at the Centre for Scientific Functional Exploration (CEFOS) according to protocols approved under agreement No. J1305505, issued by the French Ministry of Agriculture, and in compliance with the ARRIVE guidelines. Experiments were performed on intact hippocampi taken from FVB NG mice of both sexes between postnatal (P) days 5 and 7 (P0 was the day of birth); in total, seven recordings were obtained from four mice (distinct animals/hippocampi). The brain was rapidly extracted from the skull and transferred to oxygenated (95% O2/5% CO2) ice cold (4°C) artificial cerebrospinal fluid (aCSF) containing (in mM): NaCl 126; KCl 3.5; CaCl2 2.0; MgCl2 1.3; NaHCO3 25; NaHPO4 1.2; glucose 10 (pH = 7.3). After at least 2 hr of incubation at room temperature in aCSF, hippocampi were transferred to the recording chamber. To ensure a high perfusion rate (15 ml/min) we used a closed-circuit perfusion system with recycling driven by a peristaltic pump: 250 ml of solution was used per hippocampus and per condition. The pH (7.3) and temperature (33°C) were controlled during all experiments. After 30 min baseline recording in Mg2+ containing aCSF solution, the media was switched to one without added Mg2+ ion. In this condition, the extracellular concentration of Mg2+ may be influenced by other constituents of the aCSF, possibly near 0.08 mM (Mody et al., 1987). Therefore, we use the term low-Mg2+ aCSF rather than zero-Mg2+ aCSF.

Magnesium removal as a mean to influence the external potassium dynamics

Request a detailed protocol

The model of epileptic discharges presented in our study was first introduced over 20 years ago (Quilichini et al., 2002) and has since become a well-established paradigm for screening potential antiepileptic drugs and research on the mechanism of epileptic seizure (Quilichini et al., 2003).

The membrane of hippocampal neurons is equipped with N-methyl-D-aspartate type glutamate receptors (NMDARs). These receptors have a very high affinity for glutamate and can, in principle, be activated by ambient glutamate present at low concentrations in the brain extracellular fluid. Under normal physiological conditions, this activation does not occur because extracellular magnesium ions (Mg2+) block the NMDAR channel at membrane potentials more negative than about –50 mV; this voltage-dependent block prevents receptor activation at rest. When extracellular magnesium is removed, the block is relieved, allowing NMDARs to be activated, leading to neuronal depolarization toward the action potential threshold (Coan and Collingridge, 1985). In addition, as a divalent cation, Mg2+ interacts with the negatively charged neuronal membrane, contributing to the stabilization of the resting membrane potential. Lowering extracellular magnesium concentration disrupts this effect, resulting in membrane depolarization (Isaev et al., 2012).

Consequently, magnesium removal not only facilitates NMDAR-dependent depolarization, but also directly depolarizes neurons. This depolarization increases the driving force for outward potassium currents through K+ channels, meaning that variations in Mg2+ can indirectly influence external potassium dynamics during neuronal activity.

Electrical activity monitoring

Request a detailed protocol

Extracellular glass electrode filled with low Mg2+ aCSF was placed into the mid CA1 region of the ventral hippocampus. LFPs were amplified with a MultiClamp700B (Molecular Devices, San Jose, CA) amplifier for DC-coupled recording, then digitized with a Digidata 1440 (Molecular Devices, San Jose, CA), stored on the hard drive of the personal computer and displayed using PClamp 9 software (Molecular Devices, San Jose, CA).

External potassium concentration measurements

Request a detailed protocol

Potassium-selective microelectrodes were prepared using the method described by Heinemann and Arens, 1992. In brief, electrodes were pulled from double-barrel theta glass (TG150-4, Warner Instruments, Hamden, CT, USA). The reference barrel was filled with 154 mM NaCl solution. The silanized ion-sensitive barrel tip (5% trimethyl-1-chlorosilane in dichloromethane)was filled with a potassium ionophore I cocktail A (60031 Fluka distributed by Sigma-Aldrich, Lyon, France) and backfilled with 100 mM KCl. Measurements of [K+]ext-dependent potentials were performed using a high-impedance differential DC amplifier (kindly provided by Dr. U Heinemann; Heinemann and Arens, 1992) equipped with negative capacitance feedback control, which permitted recordings of relatively rapid changes in [K+]ext (time constants 50–200 ms). The electrodes were calibrated before each experiment and had a sensitivity of 42–63 mV/mM. The tip of the electrode was placed into the CA1 region of the ventral hippocampus at 200–300 µm distance from the LFP recording electrode at 150–200 µm depth that, at P5–P6 age, corresponds to the Stratum Radiatum.

Biophysical neuron model

Request a detailed protocol

The membrane potential of a single neuron in the brain is related to an ion-exchange mechanism in intracellular and extracellular space. The concentrations of potassium, sodium, and chloride in the intracellular and extracellular space along with the active transport pump (Na+/K+ pump) in the cell membrane of neurons generate input currents that drive the electrical behavior in terms of the evolution of its membrane potential. The ion-exchange mechanisms in the cellular microenvironment, including local diffusion, glial buffering, ion pumps, and ion channels, have been mathematically modeled based on conductance-based ion dynamics to reflect the ‘healthy’ and seizure behaviors in single neurons (Hodgkin and Huxley, 1952; Cressman et al., 2009; Ullah et al., 2009; Depannemaecker et al., 2022b). Moreover, conductance-based couplings between the spiking neurons have been already implemented in neural mass models (Capone et al., 2019; Guerreiro et al., 2022; Coombes and Byrne, 2019; Di Volo et al., 2019; Chen and Campbell, 2022), but without an extracellular exchange mechanism. The mechanisms of ion exchange in the intracellular and extracellular space of the neuronal membrane are represented schematically in Figure 1.

This biophysical interaction and ion-exchange mechanism across the membrane of a neuronal cell can be described as an HH-type dynamical process, represented by the following dynamical system as described in Depannemaecker et al., 2022b.

(1) dVdt=−1Cm(ICl+INa+IK+Ipump),dndt=n∞(V)−nτn,dΔ[K+]intdt=−γωi(IK−2Ipump),d[K+]gdt=ϵ([K+]bath−[K+]ext).

This model represents the ion-exchange mechanism of a single conductance-based neuron in terms of membrane potential V, the potassium conductance gating variable n, intracellular potassium concentration variation Δ[K+]int, and extracellular potassium buffering by the external bath [K+]g. This mechanism considers ion exchange through the chloride, sodium, potassium, voltage-gated channels, intracellular sodium, and extracellular potassium concentration gradients and leak currents. The intrinsic ionic currents, the sodium-potassium pump current, and potassium diffusion regulate the different ion concentrations. The Nernst equation was used to couple the neuron’s membrane potential with the concentrations of the ionic currents. This mechanism gives rise to a slow–fast dynamical system in which the membrane potential V and potassium conductance gating variable n constitute the fast subsystem and the slow subsystem is represented in terms of the variation in the intracellular potassium concentration Δ[K+]int and extracellular potassium buffering by the external bath [K+]g (Figure 2b). In Equation 1, the input currents due to different ionic substances and pump are represented as follows Cressman et al., 2009; Depannemaecker et al., 2022b:

(2) INa=(gNa,l+gNam∞(V)h(n))(V−26.64log⁡([Na+]ext[Na+]int)),IK=(gK,l+gKn)(V−26.64log⁡([K+]ext[K+]int)),ICl=gCl(V+26.64log⁡([Cl−]ext[Cl−]int)),Ipump=ρ11+exp⁡(10.5−0.5[Na+]int)11+exp⁡(5.5−[K+]ext).

The gating functions for the conductance are modeled as

(3) n∞(V)=11+exp⁡(−30−V5),m∞(V)=11+exp⁡(−22−V15),h(n)=1.1−11+exp⁡(−8(n−0.4)).

Notice that, compared to Depannemaecker et al., 2022b, these equations were reparametrized to account for changes in timescales (see also Table 1). In this model the concentration of chloride ions inside and outside of the membrane is invariant, that is, [Cl−]ext=[Cl−]0,ext and [Cl−]int=[Cl−]0,int. The extracellular and intracellular concentrations of sodium and potassium ions are represented in terms of the state variables (Δ[K+]int,[K+]g) as follows Depannemaecker et al., 2022b:

(4) Δ[Na+]int=−Δ[K+]int,Δ[Na+]ext=−βΔ[Na+]int,Δ[K+]ext=−βΔ[K+]int,[K+]int=[K+]0,int+Δ[K+]int,[Na]int=[Na+]0,int+Δ[Na+]int,[Na+]ext=[Na+]0,ext+Δ[Na+]ext,[K+]ext=[K+]0,ext+Δ[K+]ext+[K+]g.

The biophysically relevant values of the parameters could be obtained from several previous studies, both from in vivo and in vitro experiments (Cressman et al., 2009). The parameters used for the single-neuron simulations are shown in Table 1, and allow the reproduction of several bursting and oscillatory regimes by varying the parameter [K+]bath (Figure 2a).

Network of coupled biophysical neurons

Request a detailed protocol

We aim to develop a mean-field model for a heterogeneous population of the all-to-all coupled biophysical neurons described by Equation 1. Our strategy consists of using a series of approximations and heuristic arguments to obtain a putative functional form for the mean-field model, which is capable of simulating the mean activity of a large population of HH-type neurons operating near synchrony. The basic mechanism of synaptic transmission can be described by the arrival of a spike at the presynaptic terminal, resulting in an influx of calcium which depolarizes the presynaptic terminals and thus stimulates the release of neurotransmitters into the synaptic cleft. The neurotransmitters play a crucial role in opening the postsynaptic channels, which causes the flow of ionic currents across the membrane. This flow creates an (inhibitory or excitatory) post-synaptic potential in the efferent neuron. We take this into account by adding a synaptic input current in the membrane potential equation at neuron i using a conductance-based model

(5) Isyn,i=gsyn,i(t)(Esyn−Vi),

where gsyn,i is the synaptic conductance, Vi is the postsynaptic potential, and Esyn is the synaptic reversal potential at which the direction of the net current flow reverses. In the present work, we fixed Esyn=0 mV, corresponding to excitatory neurons. Lower reversal potential values (e.g., Esyn=−80 mV) model inhibitory neurons. Typically, at the arrival time tjk of the kth spike from neuron j, the synaptic conductance quickly rises and consequently undergoes an exponential decay with some rate constant. This mechanism is captured by defining the conductance dynamics gsyn,i∝Jsi(t), where J is a constant conductance value and the adimensional dynamical variable si(t) is the solution of the following equation (see, e.g., Byrne et al., 2017):

(6) Q^si(t)=ssynN∑j=1N∑kWijδ(t−tjk),

where ssyn=1 is the synaptic strength, Wij is the synaptic weight which here we assume to be equal to 1 for all neurons i and j, and the operator Q^ can be defined with various levels of accuracy as

(7) Q^=1,Q^=(1+τaddt), Q^=(1+τaddt)(1+τbddt), …

Here, for simplicity, we select the first model, which leads to

(8) si=1N∑j=1N∑kδ(t−tjk).

Under the assumption of instantaneous synaptic transmission and homogeneous all-to-all coupling, the synaptic activation variable si is the same for all neurons and corresponds to the population firing rate, which we denote by r.

In the coupled system, the membrane potential equation for neuron i is defined by:

(9) CmdVidt=−(ICl,i+INa,i+IK,i+Ipump,i)+Jsi(Esyn−Vi)+ηi

where we fix the capacitance to Cm=1nF, and the term ηi represents a heterogeneous noise current distributed according to a Lorentzian distribution with half-width Δ and location of the center at η¯

(10) g(η)=1πΔ(η−η¯)2+Δ2.

Continuity equation

Request a detailed protocol

In the continuous formulation, in the limit of large population N→∞, the density of neurons in a phase space point (V,n,Δ[K+]int ,[K+]g) at time t and excitability η is described by the population density function ρ(t,V,n,Δ[K+]int ,[K+]g,η,t), and the continuity equation holds

(11) ∂∂tρ(t,V,n,Δ[K+]int ,[K+]g,η,t)+∇⋅J(t,s,V,n,Δ[K+]int ,[K+]g,η,t)=0,

where J is the flux along the V, n, Δ[K+]int , and [K+]g directions. Since our system displays a fast and a slow subsystem (e.g., Figure 2b), we treat the variables Δ[K+]int and [K+]g as slowly varying parameters. In other words, we assume they are nearly constant over the timescale of the fast variables V and n. These slow variables are in addition considered to be mesoscopic, meaning they are identical for every neuron in the population. In this setup, the slow variables Δ[K+]int and [K+]g determine the long-term behavior of the system, while the equations for V, and n describe the rapid responses around this slow evolution. Furthermore, we consider that the potassium concentrations are homogeneous across the neuronal population, and we assume this holds also for the other ion concentrations.

Therefore, we write the population density function for the fast subsystem as

(12) ρ(t,V,n,η)=ρn(t,n∣V,η)ρV(t,V∣η)g(η)

which was here expressed in the conditional form independent of the potassium variables. The continuity equation reads

(13) ∂∂tρ(t,V,n,η)+∇⋅J(t,V,n,s,η)=0,

where the flux in the V and n directions is

(14) J(t,V,n,s,η)=(JV(t,V,n,s,η)Jn(t,V,n))=(GV(V,n,s,η)Gn(V,n))ρ(t,V,n,η),

with

(15) GV(V,n,s,η)=V˙0+Isyn+η,Gn(V,n)=n∞(V)−nτn,

and

(16) V˙0=−1Cm(ICl+INa+IK+Ipump),Isyn=JCms(Esyn−V).

We note that, under the assumption of globally shared gating and ion concentration variables across the neuronal population, the resulting mean-field equations can also be derived using simpler methods as proposed by Guerreiro et al., 2022. In this work, we follow the more general formalism of Chen and Campbell, 2022, which makes the role of key approximations (e.g., moment closure, vanishing flux at boundaries) explicit. This also facilitates potential generalizations to settings with partial heterogeneity or dynamic gating distributions.

Following (Chen and Campbell, 2022), we obtain a modified version of the continuity equation by integrating the continuity equation (Equation 13) with respect to n. Using the normalization condition on the marginal density of n, the first term gives

(17) ∫∂n∂∂tρ(t,V,n,η)dn=∂∂t(ρV(t,V∣η)∫∂nρn(t,n∣V,η)dn)g(η)=∂∂tρV(t,V∣η)g(η).

Similarly, integrating the second term in the continuity equation (Equation 13) with respect to n we obtain

(18) ∫∂n∇⋅J(t,V,n,s,η)dn=∂∂V[GV(V,⟨n∣η⟩,s,η)ρV(t,V∣η)]g(η)+∫∂n∂∂n[Gn(V,n)ρn(t,n∣V,η)]dnρV(t,V∣η)g(η)=∂∂V[GV(V,⟨n∣η⟩,s,η)ρV(t,V∣η)]g(η)+[Gn(V,n)ρn(t,n∣V,η)]n=+∞ρV(t,V∣η)g(η)−[Gn(V,n)ρn(t,n∣V,η)]n=−∞ρV(t,V∣η)g(η)=∂∂V[GV(V,⟨n∣η⟩,s,η)ρV(t,V∣η)]g(η),

where we assumed that the flux along n vanishes on the boundary ∂n. Next we assume a first-order moment closure condition for the variable n (Chen and Campbell, 2022), justified by the numerical simulations of the full network (see Figure 3—figure supplement 1) which show that for most of the neurons (close to 99% for the value of Δ same as in the other simulations) the mean of the population is well capturing the behavior of the single neurons (Nicola and Campbell, 2013). Finally, putting together these factors and assuming that n can be treated as a collective variable for each neuron (see Limitations of the model section) we arrive to

(19) ⟨n∣V,η⟩≃⟨n∣η⟩.

From here we obtain a modified version of the continuity equation for the fast subsystem

(20) ∂∂tρV(t,V∣η)=−∂∂V[GV(V,⟨n∣η⟩,s,η)ρV(t,V∣η)].

The validity of the first moment closure, Equation 19, as in Chen and Campbell, 2022, is supported by the numerical simulations, which show that, both, during the silent regime and when seizure-like events occur, ni for most neurons track the network averaged ⟨n∣V,η⟩ . In particular, it is less than 2 % of the neurons that fire while the mean is low, and vice versa, Figure 1. In less synchronized scenarios (larger Δ or smaller J), however, this value would increase, but the mean would always capture the qualitative behavior of the population.

Step-wise quadratic approximation

Request a detailed protocol

For a decoupled system described by Equation 1, considering the potassium variables as slowly varying parameters, and for any value of the gating variable n, the equation for the membrane potential V˙0 (Equation 16) has the profile of a cubic-like function (Figure 2c). Here, we approximate this function as

(21) V˙0≃{R−(V−c−)2+I− ;      ∀V0≤V⋆R+(V−c+)2+I+ ;      ∀V0>V⋆

corresponding to two parabolas with opposite curvature (R−>0,R+<0), centered at c− and c+ and shifted by I+ and I−, respectively (Figure 2d). The parameter V⋆ defines the intersection point of the two parabolas. In general, the coefficients of the parabolas (c±,I±,R±) would be functions of n(but also of Δ[K+]int and [K+]g). In the expression (Equation 20), this dependence is reduced to c±=c±(⟨n∣η⟩), I±=I±(⟨n∣η⟩), and R±=R±(⟨n∣η⟩). Here, to further simplify our problem, we consider these as free parameters, that will appear in the mean-field model, and that can be tuned to fit the average activity of a large neuronal network (see Limitations of the model section).

Steady-state solution and Lorentzian Ansatz

Request a detailed protocol

The steady-state solution for Equation 20 corresponds to

(22) ∂∂V[GV(V,⟨n∣η⟩,s,η)ρV(t,V∣η)]=0.

In standard mean-field reductions (e.g., Montbrió et al., 2015; Chen and Campbell, 2022), the function GV is quadratic in V, with two fixed points, one of which is stable, and the steady-state solution for ρV(t,V∣η) is the inverse of a quadratic i.e., a Lorentzian distribution. The Lorentzian Ansatz assumes that, at all times, such a distribution is preserved during the system’s dynamics. These conditions allow the derivation of closed solutions for a system of coupled spiking neurons. In our case, where GV has either one or three zeros, the steady state will correspond to delta functions at the fixed points, which are either one (stable) or three (of which two are stable). Therefore, a double-Lorentzian (or a piece-wise Lorenzian) could be a suitable form for ρV(t,V∣η). However, it is not clear under which conditions such an assumption would allow a solution to the continuity equation (Equation 20).

Nonetheless, in special cases, such as when the neuron population is acting close to synchrony, the neurons’ membrane potential distribution is faithfully described by a single Lorentzian peaked below or above the threshold (e.g., Figure 3a; notice that the orange distribution is bimodal, as a few neurons are idling in the subthreshold state). To proceed with the model derivation, we assume that the distribution of membrane potentials for neurons with excitability η is described by a single Lorentzian distribution (Figure 2e) on either side of the threshold V⋆

(23) ρ(V|η,t)=1πx(η,t)(V−y(η,t))2+x(η,t)2,

where x(η,t) and y(η,t) represent the width and the mean of the distribution, respectively. Thus, the validity of this approximation is limited to the cases when the neuronal population is distributed unimodally. Other cases, such as when a transition from sub- to suprathreshold activity of a subpopulation is described by the formation of two populations (e.g., Figure 2—figure supplement 1b), are not captured unless such transition is rapid enough. With these assumptions, we will derive the mean-field model equations in the next sections.

Firing rate definition

Request a detailed protocol

Previous mean-field derivations have adopted neuron models described by discontinuous quadratic equations, where the neuron is firing at a threshold Vth, after which the voltage is reset to Vreset. In these models, the threshold and reset are taken in the limit Vth=−Vreset→∞, and the quadratic expression for the membrane potential dynamics guarantees that a neuron takes a finite time to reach the threshold (Montbrió et al., 2015). Therefore, the firing rate can be defined as the flux at infinity, that is,

(24) r(t,η)=limVth→+∞ρV(t,Vth∣η)GV(Vth,⟨n∣η⟩,s,η)

In our continuous neuron model, the spiking and resetting are built into the dynamic equation and can be driven by several mechanisms, including slow potassium dynamics. Let us consider the case of a steady input current affecting the membrane potential, such that the nullcline function (here described by the piece-wise quadratic function) will have either one zero (stable fixed point) or three zeros (two stable and one unstable). In this case, we can describe the firing at the pitchfork bifurcation, where the system transitions from three fixed points to one. At the bifurcation point, a neuron with membrane potential values at −∞ will be first driven by the positive quadratic trend, until it crosses the V⋆ value, after which the negative parabolic trend will bring the neuron onto the fixed point. To have an operational definition of the firing rate and simplify our problem, we consider that once the V⋆ value is crossed the neuron membrane potential will reach the peak value at +∞, driven by the positive quadratic trend. Of course, in our model, this value is never reached in practice, because–for increasing values of the membrane potential–the negative parabola brings the membrane potential back with increasing speed. In this approximation, the firing rate definition reduces to (see Limitations of the model section).

(25) r(t,η)=limVth→∞1πx(t,η)[Vth−y(t,η)]2+x2(t,η)⋅[(R−(Vth−c−)2+I−)+Iext+Js(Esyn−Vth)]=R−πx(t,η).

Mean-field variables

Request a detailed protocol

We describe the mean-field variables as

(26) r(t)=R−πx(t)=R−π∫g(η)x(t,η)dη,=R−πx(η¯−iΔ,t),⟨V(t)⟩=∫∂ηg(η)∫∂VVρV(t,V|η)dV dη=∫∂ηg(η)y(t,η) dη=y(η¯−iΔ,t),⟨n(t)⟩=∫∂η∫∂V∫∂nnρ(t,V,n,η)dn dη dV=∫∂η⟨n|η⟩g(η)dη,

where we used the residue theorem to integrate out the parameter η (considering that there is only one pole of g(η) in the complex η-plane, we obtain that integrating out η corresponds to substituting η=η¯−iΔ).

Unlike the mean membrane potential ⟨V⟩ and the firing rate r, which can be explicitly derived from the continuity equation under the Lorentzian assumption, the expression for ⟨n(t)⟩ in Equation 26 is formal. In our mean-field model, the gating variable n is treated as a global population variable, evolving deterministically as a function of the average membrane potential. Therefore, ⟨n(t)⟩ corresponds to the collective gating variable assumed to be shared by all neurons, and is not computed by averaging distinct microscopic ni values.

Mean-field dynamics for the gating variable

Request a detailed protocol

Following (Chen and Campbell, 2022) we approximate the time derivative

(27) ⟨n⟩˙=∫∂η∫∂V∫∂nn∂∂tρ(t,V,n,η)dn dV dη=−∫∂η∫∂V∫∂nn(∂JV∂V+∂Jn∂n)dn dV dη.

The first term is integrated by parts, yielding

(28)  Term 1=∫∂η∫∂nn(∫−∞V⋆∂J−V∂V dV+∫V⋆Vmax∂J+V∂V dV)dn dη=∫∂η∫∂nn(∫−∞V⋆ dJ−V(t,V,n,s,η)+∫V⋆Vmax dJ+V(t,V,n,s,η))dn dη=∫∂η∫∂nn(J−V∣V⋆−J−V∣−∞+J+V∣Vmax−J+V∣V⋆)dn dη=0,

where J−V and J+V indicate the flux for V≤V⋆ and V≥V⋆ (preserved at V⋆, and governed by the positive and negative parabolas (Equation 21)), respectively, and where we assumed that the flux at V→−∞ and V→Vmax is zero (a safe assumption, as a neuron is pushed back toward V⋆ with (quadratic) infinite speed at these extreme values). The second term gives

(29)  Term2 =∫∂η∫∂V∫∂nn∂Jn∂n dn dV dη=∫∂η∫∂V∫∂nndJn dV dη=∫∂η∫∂V(nJn∣∂n−∫∂nJn dn)dV dη=−∫∂η∫∂V∫∂nJn dn dV dη=−∫∂η∫∂V∫∂nGn(⋅)ρ(t,V,n,η)dn dV dη=−⟨Gn(t,V,n)⟩.

Imposing ⟨n∞(V)⟩=n∞(⟨V⟩), we obtain the approximated result (see Limitations of the model section)

(30) ⟨n⟩˙=Gn(t,⟨V⟩,⟨n⟩)

Derivation of mean-field equations

Request a detailed protocol

Starting from Equations (21) and (5) the membrane potential dynamics for a coupled system of neurons is described on either side of V⋆ by the following equation:

(31) dVdt=AV2+BV+C,

with

(32) A=R,B=−2Rc−Jr,C=Rc2+I+JrE+η.

The parameters (R,c,I) were introduced in the previous sections as the curvature, center, and shift of the parabolas in the step-wise approximation. Here the indices (−, +) were dropped.

Substituting Equations (31) and (23) in the continuity Equation 20 we obtain a polynomial condition of the form

(33) k3V3+k2V2+k1V+k0=0,

Since this condition must be satisfied for all values of V, the solution for x(η,t) and y(η,t) is obtained by imposing kI≡0,∀I=0,1,2,3, leading to the mean-field equations

(34) x˙=2Axy+Bx,y˙=Ay2+By+C−Ax2.

Defining ω(η,t)=x(η,t)+iy(η,t), we can recast the above equations in the complex form

(35) ω˙=i(C−Aω2−iBω).

The mean-field equations are derived by integrating out η from the equation above. According to Equation 26, we find ω(η¯−iΔ,t)=x(t)+iy(t), that substituted in Equation 35 gives

(36) x˙={Δ+2R−(y−c−)x−Jrx; ∀y≤V⋆Δ+2R+(y−c+)x−Jrx; ∀y>V⋆,y˙={R−(y−c−)2+I−−R−x2+Jr(E−y)+η¯; ∀y≤V⋆R+(y−c+)2+I+−R+x2+Jr(E−y)+η¯; ∀y>V⋆,

which is a reduction of a system of infinite all-to-all coupled neurons described by a step-wise quadratic equation.

Neural mass model equations

Request a detailed protocol

From Equations (23) and (26), y(t)≡V¯(t) represents the average membrane potential of the population (from here on we will drop the 〈·〉 notation), while x(t) is an auxiliary variable, related to the firing rate as r(t)=R−x(t)/π (see previous section Mean-field variables). Finally, starting from Equation (36), we reverse the step-wise quadratic approximation and reintroduce the original current-based formulation (Equation 1) into the mean-field model, thereby recovering the full dynamics of the slow variables Δ[K+]int, [K]g. The mean-field approximation for a population of HH-type neurons consists of a five-dimensional system:

(37) drdt={R−Δπ+2R−(V¯−c−)r−Jr2; ∀V¯≤V⋆R−Δπ+2R+(V¯−c+)r−Jr2; ∀V¯>V⋆,dV¯dt={−1Cm(ICl+INa+IK+Ipump)−R−(π/R−)2r2+Jr(Esyn−V¯)+η¯; ∀V¯≤V⋆−1Cm(ICl+INa+IK+Ipump)−R+(π/R−)2r2+Jr(Esyn−V¯)+η¯; ∀V¯>V⋆,dndt=n∞(V¯)−nτn,dΔ[K+]intdt=−γωi(IK−2Ipump),d[K+]gdt=ϵ([K+]bath−[K+]ext}),

that we write in the compact form

(38) drdt=R−Δπ+2R±(V¯−c±)r−Jr2,dV¯dt=−1Cm(ICl+INa+IK+Ipump)−R±(π/R−)2r2+Jr(Esyn−V¯)+η¯,dndt=n∞(V¯)−nτn,dΔ[K+]intdt=−γωi(IK−2Ipump),d[K+]gdt=ϵ([K+]bath−[K+]ext}),

where the ± symbol refers to the cases V¯≤V⋆ (−) and V¯≤V⋆ (+).

Stimulation and coupling of neural mass models

Request a detailed protocol

In this work, the stimulation of the population is modeled by a common component Iinput(t) added to the membrane potential in Equation 38 as

(39) dV¯dt→dV¯dt+JIinput(Esyn−V¯).

Accordingly, N populations described by Equation 38 can be coupled together via their firing rate activity, so that the membrane potential equation at population P becomes:

(40) dVP¯dt→dVP¯dt+G∑QWPQrQ(Esyn−VP¯),

where WPQ is the structural weight proportional to the white matter fiber tract density between populations P and Q, and G is the global coupling tuning the impact of the macroscopic network structure over the local mesoscopic dynamics.

The dynamics of a single neuron and single mean-field population were simulated in Python. Spiking neural network simulations were performed using Brian2 package. Connectome-based simulations were run using The Virtual Brain platform (Sanz-Leon et al., 2015). In Supplementary file 1, we detail all the parameters and the values of the initial conditions used in the reported results.

Data availability

All code and data used in this study are publicly available on GitHub at https://github.com/grabuffo/ion-exchange-neural-mass, copy archived at Rabuffo, 2026.

References

  19. Book
    1. Coombes S
    2. Byrne Á
    (2019) Next generation neural mass models
    In: Corinto F, Torcini A, editors. Nonlinear Dynamics in Computational Neuroscience. Springer. pp. 1–16.
    https://doi.org/10.1007/978-3-319-71048-8_1
    1. Heinemann U
    2. Konnerth A
    3. Pumain R
    4. Wadman WJ
    (1986)
    Extracellular calcium and potassium concentration changes in chronic epileptic brain tissue
    Advances in Neurology 44:641–661.
  59. Book
    1. Heinemann U
    2. Arens J
    (1992)
    Production and calibration of ion-sensitive microelectrodes
    In: Kettenmann H, Grantyn R, editors. Practical Electrophysiological Methods. Wiley-Liss. pp. 206–212.
  73. Book
    1. Laing CR
    (2017) Phase oscillator network models of brain dynamics
    In: Moustafa AA, editors. Computational Models of Brain and Behavior. Wiley-Blackwell. pp. 505–517.
    https://doi.org/10.1002/9781119159193

Article and author information

Author details

  1. Giovanni Rabuffo

    Aix Marseille University, Institut de Neurosciences des Systèmes, UMR 1106, INSERM, Marseille, France
    Present address
    Center for Brain and Cognition, Universitat Pompeu Fabra, Barcelona, Spain
    Contribution
    Data curation, Software, Formal analysis, Investigation, Visualization, Writing – original draft, Project administration, Writing – review and editing
    Contributed equally with
    Abhirup Bandyopadhyay
    For correspondence
    giovanni.rabuffo@upf.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-3947-1662

  2. Abhirup Bandyopadhyay

    Aix Marseille University, Institut de Neurosciences des Systèmes, UMR 1106, INSERM, Marseille, France
    Contribution
    Conceptualization, Data curation, Software, Formal analysis, Investigation, Writing – original draft
    Contributed equally with
    Giovanni Rabuffo
    Competing interests
    No competing interests declared

  3. Carmela Calabrese

    Aix Marseille University, Institut de Neurosciences des Systèmes, UMR 1106, INSERM, Marseille, France
    Contribution
    Data curation, Software, Formal analysis, Investigation
    Competing interests
    No competing interests declared

  4. Kashyap Gudibanda

    Aix Marseille University, Institut de Neurosciences des Systèmes, UMR 1106, INSERM, Marseille, France
    Contribution
    Data curation, Software, Formal analysis, Investigation
    Competing interests
    No competing interests declared

  5. Damien Depannemaecker

    Aix Marseille University, Institut de Neurosciences des Systèmes, UMR 1106, INSERM, Marseille, France
    Contribution
    Data curation, Software, Formal analysis, Investigation
    Competing interests
    No competing interests declared

  6. Lavinia Mitiko Takarabe

    Aix Marseille University, Institut de Neurosciences des Systèmes, UMR 1106, INSERM, Marseille, France
    Contribution
    Data curation, Software, Formal analysis, Investigation
    Competing interests
    No competing interests declared

  7. Sourin Chatterjee

    Aix Marseille University, Institut de Neurosciences des Systèmes, UMR 1106, INSERM, Marseille, France
    Contribution
    Data curation, Software, Formal analysis, Investigation
    Competing interests
    No competing interests declared

  8. Maria Luisa Saggio

    Aix Marseille University, Institut de Neurosciences des Systèmes, UMR 1106, INSERM, Marseille, France
    Contribution
    Data curation, Software, Formal analysis, Investigation
    Competing interests
    No competing interests declared

  9. Mathieu Desroches

    1. MathNeuro, INRIA, University of Montpellier, Montpellier, France
    2. MCENteam, Basque Center for Applied Mathematics (BCAM), Bilbao, Spain
    Contribution
    Formal analysis, Investigation
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-9325-4207

  10. Anton Ivanov

    Aix Marseille University, Institut de Neurosciences des Systèmes, UMR 1106, INSERM, Marseille, France
    Contribution
    Data curation, Methodology, Writing – original draft
    Competing interests
    No competing interests declared

  11. Marja-Leena Linne

    Faculty of Medicine and Health Technology, Tampere University, Tampere, Finland
    Contribution
    Supervision, Writing – original draft
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-2577-7329

  12. Christophe Bernard

    Aix Marseille University, Institut de Neurosciences des Systèmes, UMR 1106, INSERM, Marseille, France
    Contribution
    Supervision, Writing – original draft
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-3014-1966

  13. Spase Petkoski

    Aix Marseille University, Institut de Neurosciences des Systèmes, UMR 1106, INSERM, Marseille, France
    Contribution
    Conceptualization, Data curation, Software, Formal analysis, Supervision, Validation, Investigation, Methodology, Writing – original draft, Project administration, Writing – review and editing
    Contributed equally with
    Viktor Jirsa
    For correspondence
    spase.petkoski@univ-amu.fr
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-4540-6293

  14. Viktor Jirsa

    Aix Marseille University, Institut de Neurosciences des Systèmes, UMR 1106, INSERM, Marseille, France
    Contribution
    Conceptualization, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Writing – original draft, Project administration, Writing – review and editing
    Contributed equally with
    Spase Petkoski
    For correspondence
    viktor.jirsa@univ-amu.fr
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8251-8860

Funding

Horizon 2020 Framework Programme

https://doi.org/10.3030/945539
  • Viktor Jirsa

Horizon 2020 Framework Programme

https://doi.org/10.3030/826421
  • Viktor Jirsa

Horizon 2020 Framework Programme

https://doi.org/10.3030/101137289
  • Viktor Jirsa

Horizon 2020 Framework Programme

https://doi.org/10.3030/101147319
  • Viktor Jirsa

Agence Nationale de la Recherche

https://doi.org/10.67599/anr-22-pesn-0012
  • Viktor Jirsa
  • Christophe Bernard

Agence Nationale de la Recherche

https://doi.org/10.67599/anr-24-rrii-0005
  • Viktor Jirsa

Agence Nationale de la Recherche (ANR-17-CE37-0001-01)

  • Christophe Bernard

Agence Nationale de la Recherche (ANR-20-NEUC-0005-02)

  • Christophe Bernard

HORIZON EUROPE Marie Sklodowska-Curie Actions

https://doi.org/10.3030/101199894
  • Giovanni Rabuffo

Eusko Jaurlaritza (BERC 2022-2025)

  • Mathieu Desroches

Agencia Estatal de Investigación (BCAM Severo Ochoa accreditation CEX2021-001142-S / MICIU / AEI / 10.13039/501100011033)

  • Mathieu Desroches

The funders had no role in study design, data collection, and interpretation, or the decision to submit the work for publication.

Acknowledgements

This research was supported by the European Union’s Horizon 2020 research and innovation program under Specific Grant Agreement 945539 (SGA3) Human Brain Project, No. 826421 Virtual Brain Cloud,

No. 101137289 Virtual Brain Twin Project and No. 101147319 EBRAINS 2.0, and government grants managed by the Agence Nationale de la Recherche references ANR-22-PESN-0012 and ANR-24-RRII-0005 (France 2030 program). CB received support from the Agence Nationale de la recherche projects ANR-17-CE37-0001-01 and ANR-20-NEUC-0005-02. GR is supported by the Marie Skłódowska–Curie Postdoctoral Fellowship (Project CAERUS) under the European Union’s Horizon Europe research and innovation programme, grant agreement No. 101199894. MD is supported by the Basque Government through the BERC 2022–2025 program and by the Ministry of Science and Innovation: BCAM Severo Ochoa accreditation CEX2021-001142-S/MICIU/AEI/10.13039/501100011033.

Ethics

In accordance with French law, ex vivo experiments do not require approval from an ethics committee. All procedures performed on live animals (deep sedation and euthanasia by decapitation prior to brain extraction) were carried out at the Centre for Scientific Functional Exploration (CEFOS) according to protocols approved under agreement No. J1305505, issued by the French Ministry of Agriculture, and in compliance with the ARRIVE guidelines and EU Directive 2010/63/EU. Experiments were performed on intact hippocampi taken from FVB NG mice of both sexes at postnatal days 5–7; in total, seven recordings were obtained from four mice. Every effort was made to minimize animal numbers and suffering.

Version history

  1. Preprint posted:
  2. Sent for peer review:
  3. Reviewed Preprint version 1:
  4. Reviewed Preprint version 2:
  5. Version of Record published:

Cite all versions

You can cite all versions using the DOI https://doi.org/10.7554/eLife.104249. This DOI represents all versions, and will always resolve to the latest one.

Copyright

© 2025, Rabuffo, Bandyopadhyay et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,955
    views
  • 141
    downloads
  • 8
    citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Citations by DOI

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Giovanni Rabuffo
  2. Abhirup Bandyopadhyay
  3. Carmela Calabrese
  4. Kashyap Gudibanda
  5. Damien Depannemaecker
  6. Lavinia Mitiko Takarabe
  7. Sourin Chatterjee
  8. Maria Luisa Saggio
  9. Mathieu Desroches
  10. Anton Ivanov
  11. Marja-Leena Linne
  12. Christophe Bernard
  13. Spase Petkoski
  14. Viktor Jirsa
(2026)
Biophysically inspired mean-field model of neuronal populations driven by ion-exchange mechanisms
eLife 14:RP104249.
https://doi.org/10.7554/eLife.104249.3

Share this article

https://doi.org/10.7554/eLife.104249

Categories and tags

Research organism