1 Introduction

Large-scale brain dynamics can be studied in silico with network models [1]. Local activity can be represented by neural mass models [2], which coupled together through synapses, time delays, and noise [3, 4] allow the emergence of whole brain activity that can be linked to empirical neuroimaging data [5]. 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 [6–8], for the analysis and identification of chaos in brain signals [9, 10], the study of epileptic seizure genesis and propagation [11–16], 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 [17–20]. Mean-field models demonstrated a predictive value for studying the mesoscopic dynamics of neuronal populations [21], providing statistical descriptions of neuronal networks [2, 17, 22–27], which can be used to address questions related to network-level mechanisms [10, 22]. 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 [28–33]. For example, statistical population measures, such as the firing rate, can be used to assess mesoscopic dynamics [1, 5, 28, 33–39].

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 [31]. The concentrations of Na⁺, K⁺, Ca²⁺, 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 [40, 41]. Changes in extracellular ion concentrations can induce fluctuations in the awake state [42, 43] as well as switch one brain state to another [40]. 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 [44–46], a result consistently reported in modeling studies [47–49]. Increases in [K⁺]_ext are tightly controlled by astrocytes [50], which can efficiently pump out [K⁺]_ext and distribute it via their syncytium to prevent hyperexcitability [51–55]. 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 [46, 56]. In such detailed models, the buffering action of astrocytes is represented by a parameter named [K⁺]_bath [46, 51, 52, 56, 57]. 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 meanfield 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 [58]). 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-like neurons operating near synchrony (Fig 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 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.

Fig. 1

Results

1.1 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 (Fig. 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 (37) were derived by approximating a locally homogeneous network of Hodgkin–Huxley (HH) type neurons operating near synchrony, in the thermodynamic limit of an infinitely large population (see Methods section).

The single neuron equations (15) were previously derived in [56] as a simplified version of HH neurons which includes three compartments: an intracellular space (ICS), and extracellular space (ECS) and an external bath (EB) in communication through ion-fluxes (Fig. 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 (Fig. 2.a). 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 (Fig. 2.b). The mean-field model inherits the single neuron biophysical parameters (Table1), 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., E_syn = 0mV) and low (e.g., E_syn = −80mV) 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 [23]. In our case, for any fixed value of the gating and potassium variables, the membrane potential equation resembles a cubic function (Fig. 2.c). Therefore, to derive the mean-field equations we assumed that: 1) the cubic-like profile can be described by a step-wise quadratic function corresponding to two parabolas with opposite curvature (Fig. 2.d); 2) The membrane potential distribution of HH type neurons is described at each time by a Lorentzian (Fig. 2.e), which is a good approximation in situations where neurons are operating close to synchrony (e.g., Fig. 3.a), but is not a suitable description for other dynamical states better described by a bimodal distribution (e.g., Suppl. Fig. S1.a, orange distribution). Using these key approximations, we derived the primary result of this work: a closed form for the neural mass equation (Eq.37) expressed in terms of mean-field variables.

Fig. 2

Fig. 3

1.2 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 Hodgkin–Huxley type single neuron dynamics (Eq. (15)) driven by an ion-exchange mechanism, has revealed that the parameter [K⁺]_bath has a significant impact on the dynamics (Fig. 2.a). Thus, for different [K⁺]_bath we simulated a network of N = 3000 coupled neurons described by Eq. (15), coupled as in Eq.(9) with synaptic strength J = 1 and heterogeneous excitability distributed with mean and width Δ = 1. Using the mean-field equations Eq.(37), we show that the mean-field model membrane potential V_MF matches qualitatively the average membrane potential V_pop of the neuronal population (Fig. 3.b). We stress that V_MF was obtained by solving 5 coupled equations Eq.(37), while V_pop was obtained by averaging the solution of 4×3000 equations (15) coupled as in Eq.(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 strongly synchronized regimes, including spike trains, tonic spiking, bursting, seizure-like events, status epilepticus-like events, and depolarization block (such as in Fig. 3.b). Notice that in such strongly synchronized regimes the population average displays dynamics qualitatively similar to the single neuron simulation (Fig. 2.b). 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 Fig. 3.b than in the single neuron Fig. 2.a). A frequency shift is also present in the mean-field model compared to the population dynamics. However, the definition of an effective can adjust this frequency shift (Fig. 3.c). In Fig. 3.c we report the 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 (Fig. 3.b).

1.3 Comparison with in vitro experiments

Performing in vitro experiments we showed that the slow fluctuations of external potassium concentration [K⁺]_ext drive the local field potential bursting (Fig. 4.a). Potassium fluctuations are rather slow, with a period of approximately 4 minutes. The biophysical neural mass model (37) is capable of simulating such slow fluctuations by modifying the parameters (e.g, C_m = 16, τ_n = 8, ϵ = 0.0001, γ = 0.00025 ; Fig. 4.b). In Figure 4.c we show the average membrane potential and external potassium for a simulation of N = 3000 coupled HH-like neurons showing a similar behavior, although the parameters were modified to simulate shorter fluctuations for computational efficiency.

Fig. 4

Changing the parameters of the model can also reveal more complex patterns not observed in the single-neuron analysis e.g., in the form of isolated bursts in the up-state (Fig. 4.e; see next section). Such an activity pattern is qualitatively similar to patterns observed in vitro (Fig. 4.d).

1.4 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 [56]. 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 bifurcation analysis to the neural mass model (37), treating the slow variables (Δ[K⁺]_int and [K⁺]_g) as parameters and exploring the bifurcation sets of the fast subsystem by varying these parameters (see Fig. 5.a). 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 (grey area).

Fig. 5

The qualitative dynamics, illustrated in Fig. 3.b, 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 Fig. 5.b), 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 [56], 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 (REF). Similarly, in the neural mass model, seizure-like events (Fig. 3.b 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 (Fig. 5.b, right versus left panels), with new complex activity regimes, which is also supported by the simulation of new behaviors as in Fig. 4.e.

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.

1.5 Coupled Neural Masses

Finally, we show that the presented mean-field model can be used to perform large-scale network simulations of brain activity e.g., 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 (Fig. 6.a). Each population (network node) is described by the mean-field model derived in this article (Eq. (37)). Nodes A, B, C, E, 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 (Fig. 6.b) and G = 100 (Fig. 6.c). 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.

Fig. 6

1.6 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 (37), 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. 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 meanfield 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 strongly synchronous regimes (Fig. 3.a). In fact, it is not guaranteed that such a population of HH-type neurons allows for a closed formulation of the ther-modynamic 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; Fig. 4.b, 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; Fig. 4.b, 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, which might unveil new dynamical regimes and extend the validity of the thermodynamic limit beyond strongly synchronous regimes. 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 (19) and the assumptions that the functions (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 (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 (25), however, more refined definitions (e.g., accounting for a heterogeneous firing threshold) can be considered [59]. 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 (Fig. 5) and sensitivity to external coupling (Fig. 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.

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 [60], 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 [61]. Hence, an appropriate way to model recorded brain activities is to write mathematical descriptions for large neuronal populations [5], 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 [62] and mean ensemble frequency [63] for coupled oscillators, or more specifically the firing rate in the case of neuronal population [64]. 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 Hodgkin–Huxley type neurons which are characterized by ion-exchange mechanism across the cellular space. The intermediate approximation into the step-wise quadratic function allowed us to apply the analytical formalism to the complex Hodgkin–Huxley 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 meanfield 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., Fig. 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 [65]. This classification was supported by a model-based bifurcation analysis using the “Epileptor” model as a neural mass model [28]. 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 invitro. 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 [66], 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 [28, 67, 68].

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 [69, 70] or quadratic integrate and fire neurons [23, 71], and rate models [20]. 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 [72, 73] and in vivo under several conditions such as quiet waking, anesthesia, slow-wave sleep and during perceptual task across several species [74–78]. 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 [8], as observed in large-scale brain recordings [79].

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 (Fig. 5.d). 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 [80, 81]. 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 [5], did not take into consideration the extracellular space. Our mean-field model is a first step towards the integration of biophysical processes that may play a key role in controlling network behavior, as shown at the spiking network level [82].

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 [40]. Although a causal relationship is not clearly established, seizures and spreading depression, which is assumed to underlie certain forms of migraine [83, 84], are associated with large (>6 mM) and very large (>12 mM) concentrations of K⁺ [83, 85]. The invariant increase of K⁺ concentration during seizure [86], provokes a saturation of potassium buffering by the astrocytes. In fact, the extracellular concentration of K⁺ is tightly controlled by astrocytes [51, 52, 54]. Most large-scale simulations do not integrate glial cells, which make up half of the brain cells [51, 54]. Their functions are altered in most, if not all, brain disorders, in particular, epilepsy [53, 55]. Our approach indirectly accounts for the astrocytic control of extracellular K⁺. Other ion species also vary during arousal/sleep and seizures, in particular, Ca²⁺ [40, 87–89]. A decrease in Ca²⁺ will decrease neurotransmission and thus change cell-to-cell communication [87–89]. Future studies are needed to integrate Ca²⁺ 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 FigS1 will be undertaken in future studies.

To date, the Epileptor model [28] is a gold standard for inferring the location of the epileptogenic zone, with direct applications in clinics [90, 91]. 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.

2 Materials and Methods

2.1 Immature hippocampus in vitro preparation

All protocols were designed and approved according to INSERM and international guidelines for experimental animal care and use. Experiments were performed on intact hippocampi taken from FVB NG mice between postnatal (P) days 5 and 7 (P0 was the day of birth). The immature mice were decapitated rapidly, and the brain was extracted from the skull and transferred to oxygenated (95% O₂ / 5% CO₂) ice cold (4 °C) artificial cerebrospinal fluid (aCSF) containing (in mM): NaCl 126; KCl 3.5; CaCl₂ 2.0; MgCl₂ 1.3; NaHCO₃ 25; NaHPO₄ 1.2; glucose 10 (pH = 7.3). After at least 2 hours 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 Mg²⁺ containing aCSF solution, the media was switched to one without added Mg²⁺ ion. In this condition, the extracellular concentration of Mg²⁺ may be influenced by other constituents of the aCSF, possibly near 0.08 mM [92]. Therefore, we use the term low-Mg²⁺ aCSF rather than zero-Mg²⁺ aCSF.

2.2 Electrical activity monitoring

Extracellular glass electrode filled with low Mg²⁺ aCSF was placed into the mid CA1 region of the ventral hippocampus. Local field potentials (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).

2.3 External Potassium concentration Measurements

Potassium-selective microelectrodes were prepared using the method described by Heinemann and Arens [93]. 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 [93]) 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.

2.4 Biophysical neuron model

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 [46, 56, 57, 94]. The mechanisms of ion exchange in the intracellular and extracellular space of the neuronal membrane are represented schematically in Fig. 1.

This biophysical interaction and ion-exchange mechanism across the membrane of a neuronal cell can be described as a Hodgkin–Huxley type dynamical process, represented by the following dynamical system as described in [56].

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 (Fig. 2.b). In Eq. (15) the input currents due to different ionic substances and pump are represented as follows [46, 56]:

The gating functions for the conductance are modeled as

Notice that, compared to [56], these equations were reparametrized to account for changes in timescales (see also 1). In this model the concentration of chloride ions inside and outside of the membrane is invariant. i.e., [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 [56]:

The biophysically relevant values of the parameters could be obtained from several previous studies, both from in vivo and in vitro experiments [46]. 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 (Fig. 2.a).

Table 1

2.5 Network of coupled biophysical neurons

We aim to develop a mean-field model for a heterogeneous population of the all-to-all coupled biophysical neurons described by Eq. (15). 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 Hodgkin–Huxley (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

where g_syn,i is the synaptic conductance, V_i is the postsynaptic potential and E_syn is the synaptic reversal potential at which the direction of the net current flow reverses. In the present work, we fixed E_syn = 0mV, corresponding to excitatory neurons. Lower reversal potential values (e.g., E_syn = −80mV) model inhibitory neurons. Typically, at the arrival time of the k-th 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 g_syn,i ∝ Js_i(t) where J is a constant conductance value and the adimensional dynamical variable s_i(t) is the solution of the following equation (see for example [26]):

where s_syn = 1 is the synaptic strength, W_ij 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

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

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

where we fix the capacitance to C_m = 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

2.6 Continuity equation

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

where 𝒥 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., Fig. 2.b), 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. 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

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

where the flux in the V and n directions is

with

and

Following [95], we obtain a modified version of the continuity equation by integrating the continuity equation 13 with respect to n. Using the normalization condition on the marginal density of n, the first term gives

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

where we assumed that the flux along n vanishes on the boundary ∂n. Finally, putting together these factors and assuming the moment closure condition (see Limitations of the model section)

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

2.7 Stepwise quadratic approximation

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

corresponding to two parabolas with opposite curvature (R₋ > 0, R₊ < 0), centered at c₋ and c₊ and shifted by I₊ and I₋, respectively (Fig. 2.d). 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 (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).

2.8 Steady-state solution and Lorentzian Ansatz

The steady-state solution for the (20) corresponds to

In standard mean-field reductions (e.g., [23, 95]), the function G^V 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 G^V 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 (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., Fig.3.a; 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 (Fig. 2.e) on either side of the threshold V ^⋆

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., Suppl. Fig.2.b), are not captured unless such transition is rapid enough. With these assumptions, we will derive the mean-field model equations in the next sections.

2.9 Firing rate definition

Previous mean-field derivations have adopted neuron models described by discontinuous quadratic equations, where the neuron is firing at a threshold V_th, after which the voltage is reset to V_reset. In these models, the threshold and reset are taken in the limit V_th = −V_reset → ∞, and the quadratic expression for the membrane potential dynamics guarantees that a neuron takes a finite time to reach the threshold [23]. Therefore, the firing rate can be defined as the flux at infinity i.e.,

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)

2.10 Mean-field variables

We describe the mean field variables as

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.

2.11 Mean-field dynamics for the gating variable

Following [95] we approximate the time derivative

The first term is integrated by parts, yielding

where 𝒥₋ and 𝒥₊ indicate the flux for V ≤ V^⋆ and V ≥ V^⋆ (preserved at V^⋆, and governed by the positive and negative parabolas (21)), respectively, and where we assumed that the flux at V → −∞ and V → V_max 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

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

2.12 Derivation of mean-field equations

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:

with

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

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

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

The mean-field equations are derived by integrating out η from the equation above. According to (??), we find , that substituted in (35) gives

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

2.13 Neural Mass Model equations

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 replace the step-wise quadratic approximation with the original functional form (15) in the mean-field model, thus reintegrating the dynamics of the other variables Δ[K⁺]_int, [K]_g. The mean-field approximation for a population of HH-type neurons consists of a five-dimensional system:

2.14 Stimulation and Coupling of Neural Mass Models

In this work, the stimulation of the population is modeled by a common component I_input(t) added to the membrane potential in Eq. 37 as:

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

where W_{P Q} 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 [5]. In Table S2, we detail all the parameters and the values of the initial conditions used in the reported results.

Supplementary information

Fig. S1

Table S2

Data and materials availability

All data, code, and materials used in the analyses will be made available at a Github repository depending on the outcome of the editorial process.

Additional information

Funding

This research was supported by the European Union’s Horizon 2020 research and innovation program under grant agreement 945539 (SGA3) Human Brain Project and by grant agreement No. 826421 Virtual Brain Cloud. CB received support from the Agence Nationale de la recherche projects ANR-17-CE37-0001-01 and ANR-20-NEUC-0005-02.

Author contributions

G.R., A.B., S.P., and V.K.J. designed research; G.R., A.B., C.C., K.G., D.D., L.MT, ML.S., M.D., S.P., and V.K.J. performed research; G.R., A.B., C.C., K.G., D.D., L.MT, M.D. simulated the models and analyzed the results; G.R., L.MT designed the figures. A.I. performed the experiments. G.R., A.B., D.D., ML. S., A.I., C.B., ML.L., S.P., and V.K.J. wrote the paper. *G.R. and *A.B. contributed equally to this work as first authors. +S.P. and +V.K.J. contributed equally to this work as the last authors. The authors declare no competing financial interests. Correspondence should be addressed to Giovanni Rabuffo at giovanni.rabuffo@univamu.fr or Spase Petkoski at spase.petkoski@univ-amu.fr or Viktor K. Jirsa at viktor.jirsa@univ-amu.fr.