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.

Single Hodgkin-Huxley type neuron model:

(a) Different patterns of electrophysiological activities previously identified in [56] 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 sub-system 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 distri bution 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 is 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 Section Steady-state solution and Lorentzian Ansatz

Mean-field model versus Neuronal network in strongly synchronized 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 in strongly synchronized regimes (parameters in Table S2). We compare the mean membrane potential Vpop and external potassium [K+]ext of such population (in blue) with the results obtained using the meanfield 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 , represented in panel (c) (the dashed line represents the identity).

Comparison of numerical results and in vitro experiments:

(a) Experimental co-recording of extracellular local field potential (blue traces) and extracellular K+ concentration (green trace). (b) The mean-field model qualitatively reproduces the experimentally observed activity showing a modulation of the amplitude of the fast oscillations during a slow oscillation cycle of the potassium variables. (c) The mean membrane potential and external potassium concentration of a population of connected HH type neurons. The matching of experimental timescales is possible by appropriate re-tuning of the neuron parameters Cm, τn, ϵ and γ (parameters in Table S2). (d) In vitro experiments show that the LFP (blue trace) is modulated by the extracellular concentration of potassium ions (green trace), fluctuating at a slow timescale. The external [K+]ext concentration is relatively stable for several minutes (after the first black arrow that points on the onset of the low Mg2+ administration) before spontaneously transitioning towards a second similarly stable state (occurring at the double-arrow) with higher potassium concentration. From this state, after a few minutes, the ion concentration starts an oscillation that induces seizure-like events evidenced by the LFP recording. (e) The mean-field model gives rise to an emergent behavior not present in the single-neuron model. In this example, the activity is characterized by a modulation of the bursting profile presenting isolated bursts in the up state (parameters in Table S2).

Fast Subsystem Bifurcation Diagram in Two Parameters:

(a) The bifurcation diagram shows the behavior of the fast subsystem in terms of the slow subsystem variables Δ[K+]int and [K+]g, which were fixed during our analysis. Each curve represents parameter values for which specific bifurcations occur, obtained through continuation methods. Saddle Node (SN) bifurcations are shown in black, Saddle Homoclinic (SH) bifurcations in green, Fold Limit Cycle (FLC) bifurcations in red, and Hopf bifurcations in blue. The gray region highlights where oscillations can emerge. The yellow star indicates where SH and SN2 bifurcation sets intersect, a point also identified in the generic model, single neuron model, and epileptor model. (b) 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, and the star marks where SH and SN2 bifurcation sets intersect. 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.

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, 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 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.

List of parameters and their values used for the single-neuron simulation.

Quadratic approximation of nullcline geometry:

(a) Neuronal network simulations of all-to-all coupled HH-like neurons can display bimodal distribution of the membrane potential when a subpopulation transitions from sub-to supra-threshold activity. (b) By changing the values of the parabola coefficients (c−, R−, c+, R+), which describe the nullcline geometry of the model, the results of the mean-field simulation change qualitatively. We show this parameters’ effect by plotting the power spectral density for several values of these parameters (here we used [K+]bath = 15.5). We can observe that for all the parameters, frequency shifts and sudden state transitions occur.

List of parameters’ values used for each simulation in the manuscript.

If not present, the parameter is not applicable for that case, or different values have been used, as specified in the corresponding caption.