Development of Mechanistic Neural Mass (mNM) Models that Link Physiology to Mean-Field Dynamics

Brain rhythms emerge from the mean-field activity of networks of neurons. There have been many efforts to build mathematical and computational embodiments in the form of discrete cell-group activities—termed neural masses—to understand in particular the origins of evoked potentials, intrinsic patterns of activities such as theta, regulation of sleep, Parkinson’s disease related dynamics, and mimic seizure dynamics. As originally utilized, standard neural masses convert input through a sigmoidal function to a firing rate, and firing rate through a synaptic alpha function to other masses. Here we define a process to build mechanistic neural masses (mNMs) as mean-field models of microscopic membrane-type (Hodgkin Huxley type) models of different neuron types that duplicate the stability, firing rate, and associated bifurcations as function of relevant slow variables - such as extracellular potassium - and synaptic current; and whose output is both firing rate and impact on the slow variables - such as transmembrane potassium flux. Small networks composed of just excitatory and inhibitory mNMs demonstrate expected dynamical states including firing, runaway excitation and depolarization block, and these transitions change in biologically observed ways with changes in extracellular potassium and excitatory-inhibitory balance.

1 Introduction

Neural masses are mean field models of neural cell group activities. Neural mass models were first introduced by Walter Freeman in the 1970s (Freeman, 1972a; Freeman, 1975) as a model to relate the activity, represented by the field potentials or EEG from one brain region to the EEG in another region, without the intricate details of neuron-level connectivity and action-potential dynamics and timing.

In the past, networks of NMs (neural mass models NMMs) have been constructed and used to understand the origins of evoked potentials (Jansen et al., 1993; Jansen and Rit, 1995), intrinsic patterns of activities such as theta (Segneri et al., 2020), regulation of sleep (Diniz Behn and Booth, 2010; Fleshner et al., 2011; Bhattacharya, 2013; Costa et al., 2016; Schellenberger Costa et al., 2016), Parkinson’s disease related dynamics (Liu et al., 2016; Liu et al., 2017; Basu et al., 2018; Liu et al., 2021), and instabilities such as seizure dynamics (Wendling et al., 2002; Wendling, 2008; Aarabi and He, 2014; López-Cuevas et al., 2015; Kameneva et al., 2017; Köksal Ersöz et al., 2020). But their firing dynamics are disconnected from real physiological parameters and instead are assumed to be sigmoidal functions of their summed input.

It is our objective here to build discrete mean-field model elements - neural masses (NM) - whose dynamics, detailed parameters, and interactions with each other and their environment can be linked quantitatively across scales directly to membrane-level mechanics of single neurons. Such mechanistic neural masses can then be used to build network models that interact with the environment and have measures that can quantitatively be related to physiological measures. Our interest here is especially to embody some of the underlying physiological mechanics that are hypothesized to be involved in epileptic dynamics.

Epilepsy - the occurrence of recurrent spontaneous seizures - is a major human health concern, leading to significant deficit in quality of life, higher mortality, and significant economic impact. In developed countries, epilepsy rates are close to 1% of the population. There are pharmacological interventions, but these provide successful abatement of seizures in only as few as 2/3 of patients, and often have significant side effects (Löscher and Schmidt, 2011; Svalheim et al., 2015; Meador and Loring, 2016; Chen et al., 2017).

At the network level, seizure susceptibility is often described as resulting from an excitatory-inhibitory imbalance. At the cellular level, this can mechanistically result from transient inactivation of inhibitory neurons leading to runaway excitation (McCormick and Contreras, 2001).

Pre- or post-seizure side effects or related neurological phenomena include postictal generalized EEG suppression (PGES) (Somjen, 2004; Poh et al., 2012; Rajakulendran and Nashef, 2015), post-ictal amnesia (Corcoran and Thompson, 1992; Butler et al., 2007; Lanzone et al., 2018), and sudden unexplained death in epilepsy (SUDEP) (Ryvlin et al., 2013; Rajakulendran and Nashef, 2015; Noebels, 2019). Mortality in persons with epilepsy are nearly two times that of the population at large (Thurman et al., 2017). In 2015, Aiba and Nobels demonstrated that spreading depolarization (SD) invading the brainstem could cause autonomic shutdown, and therefore might be one of the mechanisms of SUDEP (Aiba and Noebels, 2015), and once looked for, spontaneous seizure-associated SD has been observed in animal models of epilepsy (Ssentongo et al., 2017; Bahari et al., 2018; Loonen et al., 2019).

Leão discovered spreading depolarization - which he denoted spreading depression - in 1944 in the context of acute induced seizure models (Leao, 1944; Leo, 1944; Leo and Morison, 1945). SD is readily-associated with a range of other neurological disorders including migraine, stroke, sub-arachnoid hemorrhage, and traumatic brain injury (Pietrobon and Moskowitz, 2014; Brennan and Pietrobon, 2018; Tolner et al., 2019).

Although many mechanisms can contribute to SD, the fundamental component is that elevation of the extracelluar potassium concentration leads individual neurons into depolarization block (DB) - a state in which the potassium channels are substantially open, but the potassium current is not enough to repolarize the membrane potential and restore the channels to the normally closed state. In this state, the potassium flux can further increase the extracellular potassium concentration and, through diffusion, induce nearby neurons into DB (Somjen, 2004).

These mechanics are well known and captured in membrane and compartment models of single neurons (Kager et al., 2000; Wei et al., 2014).

In this work we define a process to develop mechanistic neural mass (mNM) elements derived from single-cell membrane dynamics. Functionally, these mass elements should serve as plug-in replacements for commonly-used elements in existing neural mass models. These differ in that they parametrize state transitions of the single-cell models to include depolarization block, their sensitivity to slow variables such as extracellular potassium and accommodation state, and their coupling back to those slow variables. We find that small excitatotry-inhibitory networks built from those networks embody dynamical transitions that underlie seizure and spreading depolarization dynamics.

It is our expectation that network models constructed with the resulting mechanistic neural masses (mNMs) we develop can reproduce normal dynamics as well as the pathological dynamics of seizures and SD in a way that can be linked to measurable biological features such as extracellular potassium, make quantitative predictions about their effects, and meaningfully give insight as to why they have such effects. As an example, we aim to give insight as to why a modest change in excitatory-inhibitory balance yields a network that is not continually seizing or undergoing SD but is significantly more susceptible to these instabilities.

As a primary step, we outline the procedure for deriving mechanistic neural masses and demonstrate it for two canonical neuron types: a non-accommodating mammalian inhibitory neuron, based on the Wang-Buzsáki (WB) neuron (Wang and Buzsáki, 1996), and an excitatory one with accommodation based on the Pinsky-Rinzel (PR) pyramidal neuron (Pinsky and Rinzel, 1994).

We then demonstrate that a two-mass inhibitory-excitatory network built from these elements, illustrated in Figure 1 already express the range of normal balanced firing, runaway excitation and depolarization block that form the fundamentals of normal and pathological dynamics for transitions into seizure and SD, and how changes in extracellular potassium modifies these transitions. Furthermore, small changes in static excitatory-inhibitory balance, as might occur from interneuron damage, significantly change the susceptibility of this network to such transitions - explicilty making the network more sensitive to small changes in extracellular potassium. We leave for future work larger networks and coupling these mNMs to the extracellular space and its ion concentrations and to glial networks.

FIGURE 1

FIGURE 1. Schematic of coupled masses: The figure shows a couplet of an excitatory (blue) and an inhibitory (pink) mass. The two masses represent mean dynamics of ensemble of respective neuron types in terms of mean firing rates (⟨FR⟩), mean potassium ion flux (⟨K_Flux⟩) and mean membrane potential (⟨V_m⟩). The means of these dynamical quantities for the excitatory(E) and inhibitory (I) masses are estimated from the canonical excitatory and the inhibitory neuron models, respectively [K]_o is the potassium concentration of the common extracellular bath of the masses. The masses couple to each other through their firing with specific coupling strengths and interact with extracellular potassium ion bath. External inputs are represented by blue and green arrows.

2 Methods

Neural masses (Freeman, 1972a; Freeman, 1972b; Wilson and Cowan, 1972; Lopes da Silva et al., 1974; Freeman, 1975; Jansen and Rit, 1995; Cona et al., 2014) are simplifications of mean-field dynamics that are thought to represent the input-output dynamics of a sub-population of neurons. In NMMs, the input to one NM is a sum over the output from the other masses. Classically, the output of a NM is its firing rate (FR), and its contribution to input on downstream masses is the FR convolved in time with a synaptic response function, described by an alpha function (Rall et al., 1967), to yield an instantaneous input current. Classically, the FR is taken as a sigmoidal function of the total input current.

Neural mass dynamics are supposed to represent the activity of groups of hundreds to millions of neurons (Freeman, 1972a). While this could be done by integrating a like number of membrane-type models, it is quite computationally intensive. Alternately, one could replace the mass with a small number - or even one - membrane-type model, but the resulting output would suffer from the precise action potential timing of the elements which would dominate the dynamics.

Critical in the process of creating neural masses is that the dynamics are significantly simplified from the detailed firing/action potential dynamics of Hodgkin-Huxley type membrane dynamics. The simplifications achieved with neural mass dynamics are very low computational load, and loss of action potential timing sensitivity by conversion to firing rates.

To achieve this mean field approach from Hodgkin-Huxley type single neuron models - to create mechanistic neural masses - we utilize a time-scale separation approach in which the fast dynamics associated with the action potentials are separated from slower dynamics by treating those other dynamics as quasi-static. We then parameterize the dynamics of that fast model, explicitly fitting as output the firing rate, mean membrane potential, and potassium flux, as a function of those quasi-static inputs. These outputs are then used in the evolution dynamics of the slow input variables, as well as the coupling of the mass to its environment and to other masses.

2.1 Procedure to Obtain Mechanistic Neural Mass

A key idea behind the development of mNMs, is that the time-averages of the fast dynamics of the neuron models can be expressed as simple functions of the total current and quasi-static (slow) state variables. Hence, instead of looking at detailed firing response of neurons in terms of instantaneous changes in membrane potentials, potassium fluxes and firing rates, we look steady-state response of the neuron models and use them to characterize the corresponding neural masses. These steady state responses are then parametrized (fitted) with simple analytic functions of the net input current and the slow state variables. These analytical functions fully characterize our neural mass outputs in terms of mean firing rate, mean membrane potential, and mean potassium flux.

We are interested in dynamics of the masses in the bifurcation space of injected current and changes in the transmembrane potassium Nernst potential (ν_K) from its nominal value. If the nominal and actual Nernst potentials are ${\widehat{\nu }}_K$ and ν_K, respectively, the change in it is expressed as Δν_K.

$\begin{array}{l}\hfill \mathrm{\Delta }{\nu }_K={\nu }_K-\widehat{{\nu }_K}\\ \hfill =\frac{RT}{F}\left[\ln \left(\frac{{\left[K\right]}_o}{{\left[K\right]}_i}\right)-\ln \left(\frac{{\left[\widehat{K}\right]}_o}{{\left[\widehat{K}\right]}_i}\right)\right]\\ \hfill =\frac{RT}{F}\left[\ln \left(\frac{{\left[K\right]}_o}{{\left[\widehat{K}\right]}_o}\right)-\ln \left(\frac{{\left[K\right]}_i}{{\left[\widehat{K}\right]}_i}\right)\right]\\ \hfill \approx \frac{RT}{F}\ln \left(\frac{{\left[K\right]}_o}{{\left[\widehat{K}\right]}_o}\right)\end{array}(1)$

Here $\frac{RT}{F}=\frac{k_{b}T}{e}\approx 26.64$. [.]_o and [.]_i represent extracellular and intracellular concentrations, respectively. The nominal values for Nernst potentials for excitatory and inhibitory masses are −75 mV and −90 mV, respectively. For these ${\widehat{\nu }}_K$s and chosen [K]_i of 140 mM, this gives nominal extracellular concentrations as 4.77 mM and 8.38 mM, respectively. Because, the intracellular volume fraction is much larger than the extracellular volume fraction, and the intracellular potassium concentration is much larger than the extracellular concentration, the fractional changes in intracellular concentration is generally negligible. Hence, Δν_K is dominantly due to fractional changes in extracellular potassium concentration, as expressed by last line of Eq. 1.

Given a Hodgkin-Huxley type single cell neuron model (i.e., the Single-compartment Excitatory Accommodating Neuron (SEAN) and Wang- Buzsáki (WB) models, described next), the procedure for developing the corresponding mNMs is:

Identify the dynamical regimes of interest. We obtain the dynamics of the neuron model at different values of injected currents (I_inj) and ν_K to obtain membrane potential traces. Following this, we find the corresponding mean firing rate (⟨FR⟩), and mean membrane potential (⟨V_m⟩) and mean potassium flux (⟨K_flux⟩) for each values of I_inj and ν_K by averaging the dynamics. As shown in Figure 2 for both the WB and the SEAN models for discrete choices of I_inj and ν_K values, this subdivides (I_inj, ν_K) parameter space into regions of not-firing (NF), firing (F), and depolarization block (DB).

FIGURE 2

FIGURE 2. Dynamics of WB and SEAN models: Shown for the WB (A) and SEAN (B) models are the mean firing rate ⟨FR⟩ (top row); mean membrane potential ⟨V_m⟩ (middle row), and mean potassium flux ⟨K_FLux⟩ (bottom row). For the red portions of the colour maps in the first row (⟨FR⟩ = 0), the models are either not firing (NF) or in DB. The difference in firing (F) and DB regions can be understood in terms of ⟨V_m⟩, and ⟨K_FLux⟩ colour maps. While the NF region shows very low negative values for ⟨V_m⟩, the DB region has them at higher negative values. Similarly, in NF region, the ⟨K_FLux⟩ is either negative or zero, in the DB region it is finite and always positive. The nominal values of the Nernst potentials correspond to Δν_K = 0 for both the models, and these are ν_K = −90 mV for the WB model, and ν_K = −75 mV for the SEAN model.

The depolarization block region that occurs at higher values of I_inj and ν_K is characteristically different from not-firing region that occurs at low I_inj values. While the NF region has either negative or zero ⟨K_flux⟩, the DB region has finite and positive ⟨K_flux⟩. Similarly, the ⟨V_m⟩ shows very negative values in the NF region, whereas it shows higher values in the DB region.

Identification of bifurcation boundaries. After identification of the three dynamical regimes of interest, we then analytically identify the bifurcation boundaries.

For the WB and SEAN models, we identify the F and DB thresholds in I_inj and ν_K state space. For this, we find stability of the fixed points of the models using eigenvalues of the model Jacobian. In practice the fixed points in these models are more easily solved as functions of V_m. The point where the largest eigenvalue crosses zero marks the transitions between stable and unstable regions. The two boundaries of the unstable triangular region are shown as blue (ISS_F) and red (ISS_DB) lines in top row of Figure 3 for both the models. The boundaries are the steady state currents at F onset and DB onset. We can observe from this figure that for each ν_K (also each row in the heat map in Figure 2), there are two values of currents that mark the firing (unstable) region. At the lower of these current values the ⟨FR⟩ is zero and at the higher values the ⟨FR⟩ is at its maximum value, for that ν_K.

FIGURE 3

FIGURE 3. Parametrization of the neural masses for the WB (A) and SEAN (B) models: In the top row are shown the boundaries of the firing (F) region in terms of steady state current (ISS) and change in potassium Nernst potential $\mathrm{\Delta }{\nu }_K$. The blue and red boundaries (ISSF and ISSDB, respectively) are obtained from the steady state analysis of the WB model (refer to the supplementary). We parametrize the maximal firing rate, and related outputs along the green line (ISSDB − Δ). Shown in the subsequent rows are the mean firing rate <FR>, mean transmembrane potential <V_m>, and mean potassium flux <K_flux>, as a function $\mathrm{\Delta }{\nu }_K$ along the boundary at (ISSDB − Δ) (closed symbols), along with the polynomial fits (solid lines), used to parameterize the masses. The nominal values of the Nernst potentials correspond to $\mathrm{\Delta }{\nu }_K=0$ for both the models. The nominal values are$\mathrm{\Delta }{\nu }_K=-90$ mV for the WB model, and $\mathrm{\Delta }{\nu }_K=75$ mV for the SEAN model.

Note that while the dynamics in the lower triangular region bounded by the ISS_F and ISS_DB lines is characterized by unstable fixed point dynamics, the upper triangular regions (at high Δν_K are characterized by two stable fixed points - one that is NF and the other in DB.

Find mean dynamical quantities near bifurcation boundaries. After having obtained the thresholds, we find the mean of dynamical variables (⟨FR⟩, ⟨V_m⟩, and ⟨K_flux⟩) near the thresholds (ISS_F and ISS_DB). We use these means to parameterize these quantities within the enclosed firing region. Outside these regions the means are analytically defined by the stable fixed point solutions.

The mean firing rate is a smooth, nearly square-root shaped function of ISS, except near the DB boundary (refer Supplementary Figure S1), where it increases sharply and the models produce incomplete action potentials. Therefore, in practice, we compute these averages on a line that is slightly shifted to (ISS_DB(ν_K) − Δ), as shown in green in Figure 3. Moreover, by using this left shifted boundary, we avoid taking very noisy and high firing rates with incomplete action potentials that likely do not produce synaptic transmission to postsynaptic cells. Averages quantities along the green boundary, for a range of changes in ν_K, for these models are shown in next three rows of Figure 3. We treat the ⟨FR⟩ and other mean dynamical quantities at ISS_DB − Δ as their “maximum” values in the firing region before the model goes into DB.

Neural Mass parametrization: Finding functional fits to time-averaged dynamical quantities. After obtaining the fairly smooth boundaries of unstable region and mean dynamical quantities as a function of Δν_K, we find and polynomial functions that parameterize the bifurcation boundaries ISS_F(ν_k) and ISS_DB, as well as ⟨FR⟩, ⟨V_m⟩, and ⟨K_flux (ν_k)⟩ along these boundaries. So, for any given ν_k, these analytical functions would give the F, DB thresholds, and “maximum” values of the dynamical quantities. These functions are then used to parametrize the dynamical quantities between the boundaries, which appear with nearly square-root behavior. These resulting equations form our mNM behavior.

Check Fidelity of neural mass parametrization. Within the firing region, defined by the parametrized boundaries $IS{S}_F(\mathrm{\Delta }{\nu }_K)$ and $IS{S}_{DB}(\mathrm{\Delta }{\nu }_K)$, we find that the neural mass outputs $FR$, $V_m$, and $K_{flux}$, are well approximated by square root functions of the normalized injected current times the parametrized maximal value of the quantity found along the DB boundary. To check fidelity of this parameterization, we then compare the neural mass outputs with the means of the same quantities obtained from the ordinary differential Equation (ODE) simulations of the respective neuron models to ascertain if the mNM outputs are in the reasonable limits of the error.

2.2 The Wang-Buzsáki and Single-Compartment Excitatory Accommodating Neuron Neuron Models Used

In the next two sections we give the ODE equations and details of the WB and SEAN neuron models that we use to obtain the I and E neural masses.

2.2.1 Wang-Buzsa´ki Model for Inhibitory Neuron

The model equations of the Wang-Buzsáki (Wang and Buzsáki, 1996) neuron are as follows:

$\begin{array}{l}\hfill \frac{d{V}_m}{dt}=\frac{1}{C_m}\left(I_{\mathit{inj}}-G_K\left(V_m-{\nu }_K\right)-G_{Na}\left(V_m-{\nu }_{Na}\right)-G_L\left(V_m-{\nu }_l\right)-I_{\mathit{pump}}\right)\\ \hfill \frac{dn}{dt}=\psi \left({\alpha }_n\left(V_m\right)\left(1-n\right)-n{\beta }_n\left(V_m\right)\right)\\ \hfill \frac{dh}{dt}=\psi \left({\alpha }_h\left(V_m\right)\left(1-h\right)-h{\beta }_h\left(V_m\right)\right)\end{array}(2)$

where,

$\begin{array}{ll}\hfill G_K& =g_k{n}^4\hfill \\ \hfill G_{Na}& =g_{Na}{m}_{\infty }^{3}h\hfill \\ \hfill G_L& =g_L\hfill \\ \hfill I_{pump}& =\frac{\rho }{\left(1+\exp (\frac{{[\widehat{N}a]}_i-{[Na]}_i}{3})\right)\left(1+\exp ({[\widehat{K}]}_o-{[K]}_o)\right)}\hfill \end{array}(3)$

The expression for the pump current and the value of pump strength ρ = 1.25 is taken from (Wei et al., 2014). Nominal ion concentrations for intracellular sodium and extracellular potassium are ${\left[\widehat{N}a\right]}_i=25\mathrm{m}\mathrm{M}$, ${\left[\widehat{K}\right]}_o=4.77\mathrm{m}\mathrm{M}$, and the nominal Nernst potentials for these ions are ${\widehat{\nu }}_K=-90\mathrm{m}\mathrm{V},{\widehat{\nu }}_{Na}=55\mathrm{m}\mathrm{V}$. The Nernst potential for sodium ion is assumed to be constant, i.e., ${\nu }_{Na}={\widehat{\nu }}_{Na}$. Non-hatted versions of concentrations incorporate the changes in Nernst potentials. In this work, we assume that the intracellular concentration remain constant at all times, hence in the pump current (I_pump), the changes in ν_K is fully accounted by [K]_o.

The activation functions for the ion channels are given by the following equations:

$\begin{array}{l}\hfill {\alpha }_m\left(V_m\right)=\frac{a_{1_W}\left(V_m+V_{1_W}\right)}{\left(1-\exp \left(-\left(V_m+V_{1_W}\right)/b_{1_W}\right)\right)}\\ \hfill {\beta }_m\left(V_m\right)=a_{2_W}\exp \left(-\left(V_m+V_{2_W}\right)/b_{2_W}\right)\\ \hfill {\alpha }_n\left(V_m\right)=\frac{a_{3_W}\left(V_m+V_{3_W}\right)}{\left(1-\exp \left(-\left(V_m+V_{3_W}\right)/b_{3_W}\right)\right)}\\ \hfill {\beta }_n\left(V_m\right)=a_{4_W}\exp \left(-\left(V_m+V_{4_W}\right)/b_{4_W}\right)\\ \hfill {\alpha }_h\left(V_m\right)=a_{5_W}\exp \left(-\left(V_m+V_{5_W}\right)/b_{5_W}\right)\\ \hfill {\beta }_h\left(V_m\right)=\frac{1}{\left(1+\exp \left(-\left(V_m+V_{6_W}\right)/b_{6_W}\right)\right)}\end{array}(4)$

The values of the model parameters are given in the first column of Supplementary Table S1. The mean outputs of the model obtained from the ODE simulations are depicted as colour maps in Figure 2A. The ⟨K_Flux⟩ current has contributions from active potassium current (G_K(V_m − ν_K)) and the pump current.

The fixed point of the system of WB model equations for a given membrane potential are given by (V_m, n_∞(V_m), m_∞(V_m), h_∞(V_m)), and its stability is computed from the Jacobian of the dynamics at that point. At the fixed point, the input current I_SS for a given V_m and ν_K is given by:

$I_{SS}=g_K{n}_{\infty }^4\left(V_m\right)\left(V_m-{\nu }_K\right)+g_{Na}{m}_{\infty }^3\left(V_m\right)h_{\infty }\left(V_m\right)\left(V_m-V_{Na}\right)+g_L\left(V_m-V_l\right)+I_{\mathit{pump}}(5)$

The boundaries of the unstable region in terms of the Firing boundary (ISS_F) and Depolarization block boundary (ISS_DB) are calculated using the above expression and depicted in the left column of Figure 3.

2.2.2 Single-Compartment Excitatory Accommodating Neuron Model

We obtain the SEAN model from the Pinksy-Rinzel model (Pinsky and Rinzel, 1994) by converting it to a single-compartment model and retaining major active currents along with the accommodation (I_K−AHP) and the pump current (I_pump).

The equations of the SEAN model are as follows:

$\begin{array}{l}\hfill \frac{d{V}_m}{dt}=\frac{1}{C_m}\left(I_{\mathit{inj}}-G_K\left(V_m-{\nu }_K\right)-G_{Na}\left(V_m-{\nu }_{Na}\right)-G_L\left(V_m-{\nu }_l\right)-I_{\mathit{pump}}-I_{K-AHP}\right)\\ \hfill \frac{dn}{dt}={\alpha }_n\left(V_m\right)\left(1-n\right)-n{\beta }_n\left(V_m\right)\\ \hfill \frac{dh}{dt}={\alpha }_h\left(V_m\right)\left(1-h\right)-h{\beta }_h\left(V_m\right)\end{array}(6)$

where,

$\begin{array}{ll}\hfill G_K& =g_k{n}^4\hfill \\ \hfill G_{Na}& =g_{Na}{m}_{\infty }^{3}h\hfill \\ \hfill G_L& =g_L\hfill \\ \hfill I_{pump}& =\frac{\rho }{\left(1+\exp (\frac{{[\widehat{N}a]}_i-{[Na]}_i}{3})\right)\left(1+\exp ({[\widehat{K}]}_o-{[K]}_o)\right)}\hfill \\ \hfill I_{K-AHP}& =G_{K-AHP}{q}_{\infty }(Ca)(V_m-{\nu }_K)\hfill \end{array}(7)$

where ρ = 1.25 is the pump strength. ${[\widehat{N}a]}_i=25$ mM, ${\left[\widehat{K}\right]}_o=8.38$ mM stand for intracellular sodium and extracellular potassium concentrations for nominal Nernst potentials for these ions (ν_K = −75 mV, ν_Na = 60 mV). Non-hatted versions of these quantities incorporate the changes in Nernst potentials.

The activation variables of different ions as a function V_m are:

$\begin{array}{l}\hfill {\alpha }_m\left({\mathit{V}}_{\mathit{m}}\right)=\frac{-a_{1_S}\left(V_{1_S}+V_m\right)}{\left(\exp \left(-\left(V_{1_S}+V_m\right)/b_{1_S}\right)-1.0\right)}\\ \hfill {\beta }_m\left({\mathit{V}}_{\mathit{m}}\right)=\frac{a_{2_S}\left(V_m+V_{2_S}\right)}{\left(\exp \left(\left(V_m+V_{2_S}\right)/b_{2_S}\right)-1.0\right)}\\ \hfill {\alpha }_n\left({\mathit{V}}_{\mathit{m}}\right)=\frac{-a_{3_S}\left(V_{3_S}+V_m\right)}{\left(\exp \left(-\left(V_{3_S}+v\right)/b_{3_S}\right)-1.0\right)}\\ \hfill {\beta }_n\left({\mathit{V}}_{\mathit{m}}\right)=a_{4_S}\exp \left(-\left(V_{4_S}+V_m\right)/b_{4_S}\right)\\ \hfill {\alpha }_h\left({\mathit{V}}_{\mathit{m}}\right)=a_{5_S}\exp \left(-\left(V_{5_S}+V_m\right)/b_{5_S}\right)\\ \hfill {\beta }_h\left({\mathit{V}}_{\mathit{m}}\right)=\frac{a_{6_S}}{\left(1.0+\exp \left(-\left(V_{6_S}+V_m\right)/b_{6_S}\right)\right)}\end{array}(8)$

The potassium after-hyper-polarization current I_K−AHP is the current responsible for accommodation in this model and calcium activation variable q is assumed to be at its steady value:

$q_{\infty }\left(Ca\right)=\frac{{\alpha }_q}{\left({\alpha }_q+{\beta }_q\right)}(9)$

where α_q = min (0.00002∗Ca, 0.01) and β_q = 0.001.

As in (Pinsky and Rinzel, 1994), the calcium dynamics themselves evolve relatively slowly according to:

$\frac{dCa}{dt}=-0.13I_{Ca}-0.075Ca(10)$

We then make the follow approximation to separate the effects of the I_K−AHP into fast dynamical (I_K−AHP,fast) effects computed at the nominal state and slow modulatory effects that depend on calcium (I_K−AHP,slow (Ca)):

$\begin{array}{l}\hfill I_{K-AHP}=g_{K-AHP}\left(\widehat{q}+q_{\infty }\left(Ca\right)-\widehat{q}\right)\left(V_m-{\nu }_K\right)\\ \hfill =g_{K-AHP}\widehat{q}\left(V_m-{\nu }_K\right)+g_{K-AHP}\left(q_{\infty }\left(Ca\right)-\widehat{q}\right)\left(V_m-{\nu }_K\right)\\ \hfill \approx g_{K-AHP}\widehat{q}\left(V_m-{\nu }_K\right)+\langle g_{K-AHP}\left(q_{\infty }\left(Ca\right)-\widehat{q}\right)\left(V_m-{\nu }_K\right)\rangle \\ \hfill \approx g_{K-AHP}\widehat{q}\left(V_m-{\nu }_K\right)+g_{K-AHP}\left(q_{\infty }\left(Ca\right)-\widehat{q}\right)\left(\langle V_m\rangle -{\nu }_K\right)\\ \hfill =I_{K-AHP,fast}+I_{K-AHP,slow}\left(Ca\right)\end{array}(11)$

In practice, the mNM parametrization is computed just using the fast dynamics I_K−AHP,fast, where we have chosen $\widehat{q}=q_{\infty }(Ca=0.2)$. The slow dynamics appear as an additive correction to the input current I_inj,eff = I_inj + I_K−AHP,slow.

The values of the model parameters are given in the second column of Supplementary Table S1. The mean outputs of the model obtained from the ODE simulations are depicted in Figure 2B. The ⟨K_Flux⟩ current has contributions from the active potassium current (G_K(V_m − ν_K)), the pump current, and the accommodation current.

The model fixed points are determined similarly as the previous model. The steady state value of the current at given (V_m and ν_K) are determined as:

$I_{SS}=g_K{n}_{\infty }\left(V_m\right)\left(V_m-{\nu }_K\right)+g_{Na}{m}_{\infty }^3\left(V_m\right)h_{\infty }\left(V_m\right)\left(V_m-V_{Na}\right)+g_L\left(V_m-V_l\right)+I_{\mathit{pump}}+I_{K-AHP}(12)$

The firing onset (ISS_F) and DB (ISS_DB) boundaries in terms of steady state current are depicted in Figure 3.

2.3 Mechanistic Neural Mass Parametrization

We use the functional fits to the firing and DB thresholds, and the mean of dynamical variables at (adjusted) DB threshold (ISS_DB − Δ) to parameterize the mNMs. For both the neuron models, we use polynomial functions up to 3rd order to fit these quantities as a function of ν_K (x). For instance, for functional fit to ISS_F of WB model, we obtain the best linear fit (y = ax + b) with a = 0.014, b = 1.746 as parameters. All the functional best fits for both the models are shown in Supplementary Figures S2 and S3. The order of polynomial was simply chosen so as to fit the profile of the simulated data (see Supplementary Figures S2 and S3), and python’s Scipy (docs.scipy.org) module was used to obtain the best fit. These functions are used in ascertaining the limits of the three dynamical quantities, for a given ν_K.

As a next step we separately identify another set of simple functions that would approximate the values of the dynamical variables of the mass, within the predetermined limit, as a function of the injected current. These functions (example $y=M\sqrt{x}+C$) use the input current (injected current) as independent variable (x) and the predetermined limits of dynamical variables (M) to yield mean dynamical output (y) in terms of mean firing rate, mean membrane potential, or, mean potassium flux.

In summary, for any given combination of (I_inj, ν_K), we first determine the threshold currents, and maximum values of mean dynamical variables from the first set of analytical functions. Next, we use the I_inj, as the independent variable in the next set of functions and find the dynamical variables corresponding to the mass. The choice of these functions was guided by the profiles of average firing rate ⟨FR⟩, membrane potential ⟨V_m⟩, and potassium flux ⟨K_flux⟩ of individual neuron models as function of I_inj. These are shown for three different ν_K values for both the neuron models in Supplementary Figure S1.

After the neural masses are characterized, we check their fidelity by comparing the mass outputs with that of mean outputs from the ODE simulations of the neuron models across the ranges of I_inj and ν_K.

2.3.1 Determination of Thresholds for a Given Change in the Nernst Potential

The thresholds current values and maximum of the means of FR, V_m, and K_Flux within the firing range are determined using simple linear and polynomial functions of ν_K of the two masses. Hence, firstly, for a given change in the potassium Nernst potential Δν_K, the corresponding (new) Nernst potentials of the two masses are:

${\nu }_{K_I}={\widehat{\nu }}_{K_I}+\mathrm{\Delta }{\nu }_K(\mathrm{13a})$

${\nu }_{K_E}={\widehat{\nu }}_{K_E}+\mathrm{\Delta }{\nu }_K(\mathrm{13b})$

The hatted versions of the Nernst potentials (${\widehat{\nu }}_{K_I}$ and ${\widehat{\nu }}_{K_E}$) are the nominal Nernst potential values of the two masses. Note that amount of change in Nernst potential is assumed to be caused only due to change in extracellular potassium concentrations for the two masses, and these masses share a common extracellular bath (see Figure 1). Hence Δν_K is same for both the masses.

The functions that determine the firing threshold $(i_{I_{\mathit{th1}}})$, depolarization threshold $(i_{I_{\mathit{th2}}})$, maximum $⟨\text{FR}⟩(M_1^I)$, maximum $⟨V_m⟩(M_2^I)$, and maximum $⟨K_{\mathit{Flux}}⟩(M_3^I)$, as functions of ${\nu }_{K_I}$ for the I mass are:

$i_{I_{\mathit{th1}}}=a_{I_{\mathit{FRT}}}{\nu }_{K_I}+b_{I_{\mathit{FRT}}}(\mathrm{14a})$

$i_{I_{\mathit{th2}}}=a_{I_{\mathit{DBT}}}{\nu }_{K_I}+b_{I_{\mathit{DBT}}}(\mathrm{14b})$

$M_1^I=a_{I_{FR}}{\nu }_{K_I}^3+b_{I_{FR}}{\nu }_{K_I}^2+c_{I_{FR}}{\nu }_{K_I}+d_{I_{FR}}(\mathrm{14c})$

$M_2^I=a_{I_{V_m}}{\nu }_{K_I}^2+b_{I_{V_m}}{\nu }_{K_I}+c_{I_{V_m}}(\mathrm{14d})$

$M_3^I=a_{I_{KF}}{\nu }_{K_I}+b_{I_{KF}}(14e)$

Similarly, the firing threshold $(i_{E_{\mathit{th1}}})$, depolarization threshold $(i_{E_{\mathit{th2}}})$, maximum $⟨V_m⟩(M_1^E)$, maximum $⟨V_m⟩(M_2^E)$, and maximum $⟨K_{\mathit{Flux}}⟩(M_3^E)$ functions (of ${\nu }_{K_E}$) that determine the threshold values for the E mass are:

$i_{E_{\mathit{th1}}}=a_{E_{\mathit{FRT}}}{\nu }_{K_E}+b_{E_{\mathit{FRT}}}(\mathrm{15a})$

$i_{E_{\mathit{th2}}}=a_{E_{\mathit{DBT}}}{\nu }_{K_E}+b_{E_{\mathit{DBT}}}(\mathrm{15b})$

$M_1^E=a_{E_{FR}}{\nu }_{K_E}^2+b_{E_{FR}}{\nu }_{K_E}+c_{E_{FR}}(\mathrm{15c})$

$M_2^E=a_{E_{V_m}}{\nu }_{K_E}^2+b_{E_{V_m}}{\nu }_{K_E}+c_{E_{V_m}}(\mathrm{15d})$

$M_3^E=a_{E_{KF}}{\nu }_{K_E}^3+b_{E_{KF}}{\nu }_{K_E}^2+c_{E_{KF}}(15e)$

The values of the parameters appearing in the above equations Eqs. 14 and 15 are presented in Supplementary Table S2. The functions are plotted against actual thresholds for a range of ν_K values and shown in Supplementary Figures S1 and S2, for the I and the E mass, respectively, in the SI.

2.3.2 Neural Mass Functionals

The parameters of the above linear and polynomial functions can be used to determine these thresholds for any arbitrary Δν_K. These thresholds are further used in another set of functions to determine the mean dynamical outputs of the mass for any given injected current. As these functions fully characterize the neural mass with respect to the bifurcation parameters (I_inj and ν_K), we address them as the neural mass functionals. We present these functionals below for the I and the E masses.

The mean firing rate function (FR_I), mean membrane potential function $(V_{m_I})$, mean potassium flux functions (KF_I) for the I mass are:

$F{R}_I=M_1^I\sqrt{Q_I}(\mathrm{16a})$

$V_{m_I}=-M_2^I\left(\sqrt{Q_I}\right)+C1(\mathrm{16b})$

$K{F}_I=M_3^I\left(\sqrt{Q_I}\right)(\mathrm{16c})$

where, C1 = −60, $M_1^I$, $M_2^I$, and $M_3^I$ are the maximum of the respective quantities as determined for a given ${\nu }_{K_I}$ from Eq. 14. Q_I is the normalized total current that the mass receives and is expressed as:

$Q_I=\frac{X_{I_{\mathit{tot}}}-i_{I_{\mathit{th1}}}}{i_{I_{\mathit{th2}}}-i_{I_{\mathit{th1}}}}(17)$

Notice that, for the isolated mass this total current $(X_{I_{\mathit{tot}}})$ is just the injected current $(I_{I_{\mathit{inj}}})$. Alternatively, when this mass is coupled to other mass (es), the total current is the sum of injected current and the current due to finite ⟨FR⟩ of other mass (es). Hence, the expression for the total current for the mass (U), that couples to another mass (V) is:

$X_{U_{\mathit{tot}}}=I_{U_{\mathit{inj}}}+I_{U_{\mathit{pump}}}\left(\mathrm{\Delta }{\nu }_K\right)+g_{VU}F{R}_V(18)$

The mean firing rate function (FR_E), mean membrane potential function $(V_{m_E})$, mean potassium flux functions (KF_E) for the E mass are:

$F{R}_E=M_1^E\sqrt{Q_E}(\mathrm{19a})$

$V_{m_E}=-M_2^E\left(Q_E^{0.62}\right)+C1(\mathrm{19b})$

$K{F}_E=M_3^E\left(Q_E^{0.28}\right)+C2(\mathrm{19c})$

where, C1 = −65 and C2 = 0.1 and normalized total current (Q_E) is expressed as:

$Q_E=\frac{X_{E_{\mathit{tot}}}-i_{E_{\mathit{th1}}}}{i_{E_{\mathit{th2}}}-i_{E_{\mathit{th1}}}}(20)$

Here again, the total current $X_{E_{\mathit{tot}}}$ can have up to three contributions as expressed in Eq. 18.

2.4 Synaptic Coupling of Masses

So far, we have presented how we develop the neural masses of both kinds (E and I) from their respective neuron models. Next, we would like to understand how does the coupled E − I model, in which the two masses interact with each other through synaptic coupling, behaves under the conditions of changing injected currents and changes in the ν_K of both the masses.

Mass coupling can also be understood from an averaging approach from classic Hodgkin-Huxley (HH) type compartment models. We write this out both for completeness, and to illustrate a slight change from existing NMMs.

In the HH like models originally input is only in the form of current. In neuronal networks, this is mediated by synaptic currents, whose impulse-response function in response to neurotransmitter release at time t_ap has been characterized by alpha functions α(t) (Rall et al., 1967; Brown and Johnston, 1983; Traub and Miles, 1991).

$\begin{array}{l}\hfill I\left(t-t_{ap}\right)=G\alpha \left(t\right)\left(V_m-{\nu }_s\right)\\ \hfill \alpha \left(t\right)=t{e}^{-t/{\tau }_s}\end{array}(21)$

Note that the time constant τ_s is particular to a type of neurotransmitter, and we will in the future denote as a generalization α_x(t) to have time constant τ_x. Note that the term V_m − ν_s is the difference between membrane potential V_m and the Nernst-potential of carrier ions through the channel of type s. Commonly used time constants include τ_p = 10 ms for pyramidal neurons, τ_I,f = 2 ms for fast (γ-aminobutyric acid GABA_A) inhibitory, and τ_I,s = 20 ms for slow (GABA_B) inhibitory interneurons (Traub and Miles, 1991).

For a single neuron, denoted as neuron i, with firing times at its soma of T_k, its synapses onto neurons will occur at axonal terminals, and there will be a delay time D mostly proportional to the distance such that action potentials will be seen at the synapse at times T_i,k + D_i,b. Critically, delay time D_i,b is specific to the connection distance between neurons i and b.

We note that a particular neuron type, such as pyramidal neuron, will only release a single neurotransmitter type. Then the time course of postsynaptic currents from firing of that neuron to postsynaptic neuron b, will have the form:

$I_b\left(t\right)=G_{i,b}\left(V_{m,b}-{\nu }_b\right)\sum _k{\alpha }_i\left(t-T_{i,k}-D_{i,b}\right)(22)$

It is now convenient to introduce a pre-synaptic centric conductance function h_a,i(t) for the ith single cell from population a, assumed to be all of the same cell type:

$h_{a,i}\left(t\right)\equiv \sum _k{\alpha }_{a,i}\left(t-T_{a,k}\right)(23)$

In this form, any postsynaptic current in neuron b from action potentials in i will be characterized as,

$I_b\left(t\right)=G_{a,b}\left(V_{m,b}-{\nu }_b\right)h_{a,i}\left(t-D_{a,b}\right)(24)$

In our neural mass formalism, we now average over a population a, composed of neurons i, all with the same neurotrasmitter type, and assumed to have common time delay D_a,b. As such α_a,i = α_a. As a result we arrive at an equivalent averaged conductance function h_a(t),

$h_a\left(t\right)={⟨\sum _{i,k}{\alpha }_{a,i}\left(t-T_{i,k}\right)⟩}_{i,k}(\mathrm{25a})$

$={⟨\sum _{i,k}\int d\tau {\alpha }_a\left(t-\tau \right)\delta \left(\tau -T_{i,k}\right)⟩}_{i,k}(\mathrm{25b})$

$=\int d\tau {\alpha }_a\left(t-\tau \right){⟨\sum _{i,k}\delta \left(\tau -T_{i,k}\right)⟩}_{i,k}(\mathrm{25c})$

$=\int d\tau {\alpha }_a\left(t-\tau \right)F_a\left(\tau \right)(\mathrm{25d})$

Here we’ve defined the firing rate for the neuron group defined by a as F_a(τ).

It can be shown that h_a(t) can be computed without the time convolution by instead integrating the second order ordinary differential equation:

${\ddot{h}}_a=\frac{F_a-h_a}{{\tau }_a^2}-\frac{2{\dot{h}}_a}{{\tau }_a}(26)$

Note that under this parametrization, under steady state conditions F_a(t) = F_a and ${\dot{h}}_a=0$ then h_a = F_a.

For computational efficiency, for masses that are coupled into networks with both near and far connections, we only compute the synaptic dynamics once per mass, and accommodate the delays in their connections. Explicitly, for NM for population b, whose presynaptic masses are denoted as members a, will have a total input current

$I_{b,in}\left(t\right)=I_{b,pump}+\sum _a{G}_{a,b}\left(⟨V_m⟩-{\nu }_b\right)h_a\left(t-D_{a,b}\right)(27)$

and idealized instantaneous firing rate defined by our parametrization f_b (I_b,in(t), θ_b(t)).

Note that because synaptic current is slow with respect to the fast dynamics, and appears after averaging across the post-synaptic network, the term V_m in Eq. 21 is replaced with the NM state ⟨V_m⟩. This is commonly further replaced with a nominal value (as in (Deschle et al., 2021)) ${\widehat{V}}_{m,b}$. This may be rationalized by assumption that the synapses, which typically appear distant from the neural Soma modeled, are not well-described by the mean Soma membrane potential. In the latter case the term $G_{a,b}({\widehat{V}}_{m,b}-{\nu }_b)$ is constant, and replaced with an overall coupling constant g_a,b.

$I_{b,in}\left(t\right)=I_{b,pump}+\sum _a{g}_{a,b}{h}_a\left(t-D_{a,b}\right)(28)$

3 Results

As described, we have developed mNMs as parameterizations of single-neuron Hodgkin-Huxley style neurons for both inhibitory and excitatory neuron types, based on the WB and SEAN models. As results, we first demonstrate that the individual mNMs parametrizations reproduce the outputs of the demonstrated models.

We then investigate the dynamics of the simple 2-element excitatory/inhibitory network illustrated in Figure 1 and demonstrate that as a function of both input and extracellular potassium, it reveals dynamics expected to underlie seizure and spreading depression dynamics.

3.1 Mechanistic Neural Mass Fidelity

To check the fidelity of the mechanistic neural mass parametrization, we compare the outputs obtained from the I and E neural functionals with that of mean outputs from the WB and SEAN neuron ODE models. For both the neurons, we begin by choosing a range of injected currents for which the neurons fire and obtain the mean quantities (⟨FR⟩_ODE, ${⟨V_m⟩}_{\mathit{ODE}}$, and ${⟨K_{\mathit{Flux}}⟩}_{\mathit{ODE}}$) at three different ν_K values. We also find out these three quantities (⟨FR⟩_mass, ${⟨V_m⟩}_{\mathit{mass}}$, and ${⟨K_{\mathit{Flux}}⟩}_{\mathit{mass}}$), using the neural mass functionals for the same range of injected currents and for the same ν_K values. We plot these mean quantities, one from solving the ODE and other from the mass functionals against one another as in Figure 4. The green line depicts y = x line in all the six subplots. For a perfect match between the ODE and mass results quantities, the points would lie along this line. We observe a good match between the results for both kinds of masses for ⟨FR⟩ and ⟨V_m⟩ with error of less than 10% for most of the firing range of the ODE models. There are a few exceptions, however, at the onset of firing and at the DB onset, where the output of the ODE model is not exactly approximated by the mass functionals.

FIGURE 4

FIGURE 4. Fidelity of mass parametrization: The figures show a comparison of mean dynamical variables obtained from ODE simulations with those obtained from neural mass parametrization at three different values of ν_K for the I mass (A) and the E mass (B), respectively. The green line in all the plots depicts a 1:1 line for the ODE and neural mass outputs. The colours of the dots in the figure are indicative of different potassium Nernst potential values, as depicted in the middle subplots. The points parallel the x axis in the subplots at firing onset and at DB onset, depict areas where the mass model is not exactly approximating the mean outputs from the neuron models.

3.2 Dynamics of the Coupled Model: Ordinary Differential Equation Formulation

Using the methodology as described in Coupled Model section in Methods, we couple the two kinds of masses to each other, so as to study the dynamical behaviour of the E − I couplet. Here we investigate the system behaviour in terms of the E and I mass firing rates variations in potassium ion concentration, injected current and coupling strength between the masses. Hence, the mean FR of each of the masses is a function of total current expressed as combination of injected current and the mean FR of the other masses that it couples to. The ODEs of the E − I coupled model are expressed in terms of FR_E and FR_I as:

$\dot{F{R}_E}=\left({\mathrm{\Phi }}_E\left(I_E\left(t\right)\right)-F{R}_E\right)/{\lambda }_E(29)$

$\dot{F{R}_I}=\left({\mathrm{\Phi }}_I\left(I_I\left(t\right)\right)-F{R}_I\right)/{\lambda }_I(30)$

This form, which is a discretized version of the Wilson-Cowan model (1972) (Wilson and Cowan, 1972), was adopted to acknowledge that firing rates - which are an observation of the system - cannot change instantaneously. λ_E and λ_I are the time constants of the respective masses. The functions Φ are FR functions of the total currents to these masses, as given in Eqs. 16a and 19a. I_E(t) and I_I(t) are the total contributions to the E and I masses, respectively, expressed as:

$I_E\left(t\right)=I_{E,inj}+I_{E,pump}-g_{I,E}{h}_I\left(t\right)(31)$

$I_I\left(t\right)=I_{I,inj}+I_{I,pump}+g_{E,I}{h}_E\left(t\right)(32)$

where the currents I_x,inj denote the input currents from outside the network, I_x,pump denotes the pump currents, and the coupling constants g_xy and synaptic functions h_x(t), and the differential equations that govern them, follow the definitions in Section 2.4. Note that the sign of these synaptic inputs are written out explicitly, with inhibitory input appearing as a negative current.

The values of the parameters used for the simulations are, in ms, λ_E = λ_I = 1, τ_E = 10, τ_I = 20, which corresponds to slow (GabaB-type) synapses. Nominal coupling constants used here are ${\widehat{g}}_{EI}=0.2$ and ${\widehat{g}}_{IE}=0.04$. Excitatory to inhibitory coupling was kept at the nominal value, and inhibitory to excitatory coupling was varied with respect to the nominal value. Additionally, no external input was applied to the inhibitory mass I_I,inj = 0.

Some of the fundamental dynamics of this E − I network, under external injected current input I_E,inj, are shown under four different state conditions in Figure 5. Here the external input is a series of current pulses with a 250 ms ramp input, and 0.75 duty cycle, with increasing amplitudes (lower traces). The upper left panel corresponds to nominal ν_K and g_IE. At very low pulse height, this small network oscillates. We expect that addition of fast (GABA-A) inhibition could be used to tune that behavior.

FIGURE 5

FIGURE 5. Dynamics of the coupled model: The figures show the dynamical response of the coupled masses in terms of their mean firing rates (FR_E and FR_I) for a square wave stimulus (shown in black in the bottom most subplots) for Δν_k = 0 (left) and Δν_k = 15 (right). The top subplots are for nominal I to E coupling strength $g_{IE}={\widehat{g}}_{IE}$ and the middle subplots show the dynamics when $g_{IE}=\frac{{\widehat{g}}_{IE}}{2}$.

At moderate pulse height, we observe balanced firing of both the excitatory and inhibitory masses. As the pulse height increases, the inhibitory mass is driven into DB (at approximately time t = 12 s). At this point the excitatory mass fires at an even higher rate, and is equivalent to the Run-away excitation (RAE). At even higher pulse height, the excitatory mass is driven into depolarization block and both their firing rates go to zero during the pulse. We denote this network state as DB. Note that upon pulse initiation, because the pulses have a sloped front end, both masses initially fire before reaching DB.

In the left, middle panel, we have decreased g_IE by a factor of two. The result is that the transition to RAE occurs at a lower pulse input height.

In the right column we show the results for slightly elevated extracellular potassium (Δν_K = 15 mV, or $\sim 75\%$ increase). In this case the transitions to both RAE and DB occur at lower input currents, are further compressed with decreased g_IE.

The network dynamics are hysteretic with respect to injected current. This is apparent from the detailed response of the model to a symmetric ramp current function as shown in Figure 6. Here the output is not the same for increasing input as for decreasing input. A first indication of this is that the ramp points where the I mass goes in and out of DB occur at different values of $I_{E_{\mathit{inj}}}$. This is better observed in Figure 6B where the firing rates are plotted as a function of $I_{E_{\mathit{inj}}}$, and the sign of ${\dot{I}}_{E_{\mathit{inj}}}$ is denoted with solid (positive) or dashed (negative) line types. Therefore, in state space, there is bistable region in which the network can either be in stable firing (F) or RAE. We denote this region as BS.

FIGURE 6

FIGURE 6. Dynamics of the coupled model to ramp stimulus: The figures show the dynamical response of the coupled masses in terms of their mean firing rates (FR_E and FR_I) for a ramp stimulus (shown in black in part (A)) for Δν_k = 0 and $g_{EI}={\widehat{g}}_{EI}/2$. As the injected current to the E mass increases the firing rates increases proportionally, until the I mass goes into DB, which results in unbalanced firing of E mass or run-away excitation. As the current is further increased the E mass goes into DB, at around $I_{E_{\mathit{inj}}}=81\mu A/c{m}^2$. This is also accompanied by a simultaneous spike in FR_I, owing to momentary drop in current (from its DB threshold) it received from E mass firing. This non-zero FR_I is not sustained any further as the E mass goes into DB due to very high injected current. This is repeated in reverse as the current goes down from its maximum value in the second half of the stimulus. However, the point at which the I mass recovers from DB (during second half of the current), and at which it goes into DB (in the first half) are not exactly same. These hysteretic dynamics result in bi-stable (BS) region as is shown with the other three regions (F, RAE, and DB) in part (B) as a function of injected current.

3.3 Analytical Boundaries of the IE Network Dynamics

As explored in the previous section, the coupling of the E mass and I mass models can result in four dynamical states as defined not-firing NF, stable firing F, RAE, and DB, and potentially can have bistable state-space regions that can have either. In this section, we delineate analytically the boundaries between these behaviors in parameter-space composed of Δν_K, g_IE, and I_E,inj.

In the NF state, the masses do not fire due to insufficient total current for E mass firing. In the DB state the E mass does not fire because of very high total current. In the RAE state the E mass shows uncontrolled (pathological) firing due to failure of the I mass to inhibit its firing. This in turn is because the E mass’s firing rate is high enough to put the I mass into DB. In the BS region, the E mass shows two stable firing patterns for the same injected current depending on whether the I mass is in DB or not - and this depends on history. When the I mass is in the F state and injected current is on, the E mass shows one FR (lower), while when I mass is in the DB state, and the same injected current is still on, the E mass shows another FR (higher). This property marks the bi-stable region in the $I_{E_{\mathit{inj}}}$ - Δν_K bifurcation space. Also, at the boundaries of the bi-stable region, the FRs of the E mass are the same. The F state corresponds to the E mass normal firing.

The analytical boundaries of these regions are shown in Figure 7. Note that in absence of coupling to the I mass, the blue and red lines in the figure enclose the unstable fixed points region of the E mass. On the left of the blue boundary, the state is NF, the state beyond the red boundary is DB. When the coupling between E and I mass is on, the orange region limited by the blue line (firing boundary) on the left and orange line (bistability boundary) on the right is the firing region (F state). The yellow region between the orange and the green line (RAE boundary) is the bistable region (BS state). The red region between the green line and the red line (DB boundary) is the RAE region.

FIGURE 7

FIGURE 7. Analytical bounds of dynamical regimes of the coupled model: The figures (A–D) show the analytical boundaries for not firing, firing, bi-stability, run-away excitation and depolarization block regions as a function of couplings g_IE and Δν_K. The boundaries and regions enclosed by them for a range of values of Δν_K at two fixed couplings g_IE are shown in (A,C). (B,C) show the firing rates of the two masses as a function of ramp current at the two Δν_K values 0 and 15 mV (marked in grey line in parts (A,C)). The solid lines are firing rates for the increasing current and the dotted lines show the firing rate for the decreasing current.

Steps to obtain the Firing onset (FO), Bi-stability (BS), Run-away excitation (RAE) and Depolarization block (DB) boundaries are calculated as follows, where the Nernst potentials are adjusted using Eq. (13).

1) The FO boundary is obtained using the parameters to the functional fits as in the equation of Eq. 15a:

$i_{E_{\mathit{th1}}}=a_{E_{\mathit{FRT}}}{\nu }_{K_E}+b_{E_{\mathit{FRT}}}$

2) The RAE boundary is characterized by a point where I mass goes into depolarization, after reaching its maximum firing rate and as a result the E mass fires in uninhibited manner. Hence, we first find the total current to the I mass (Q_I) that results in its maximum FR using Eq. 16a.

$Q_I={\left(\frac{F{R}_I}{M_1^I}\right)}^2(33)$

Using the obtained Q_I, and Eq. 18, we find FR_E for the above Q_I as below. Notice that there is no injected current to the I mass (0 in the equation below).

$F{R}_E=\left(1/g_{EI}\right)\left(Q_I-0+I_{I_{\mathit{pump}}}\right)(34)$

Using the above obtained FR_E, we obtain the total current to the E mass, i.e., Q_E using Eq. 19a,

$Q_E={\left(\frac{F{R}_E}{M_1^E}\right)}^2(35)$

Now since we know the total current, we can again use equation. 18, to find injected current which is our I_RAE, as follows.

$I_{\mathit{RAE}}=Q_E+I_{E_{\mathit{pump}}}+g_{IE}{M}_1^I(36)$

3) The bi-stability boundary is obtained by invoking the condition: FR_E(BS) = FR_E(RAE), which implies that Q_BS = Q_RAE. Using Eq. 18, we arrive at:

$I_{BS}-g_{IE}F{R}_{I_{BS}}-I_{E_{\mathit{pump}}}=I_{\mathit{RAE}}-g_{IE}{M}_1^I-I_{E_{\mathit{pump}}}(37)$

Since FR_I at BS zero, the second term on the left hand side of the above equation goes to zero.

Hence, we arrive at:

$I_{BS}=-g_{IE}{M}_1^I+I_{\mathit{RAE}}(38)$

4) The DB boundary is using the parameters of the fit to the DB threshold as in Eq. 15b.

$i_{E_{\mathit{th2}}}=a_{E_{\mathit{DBT}}}{\nu }_{K_E}+b_{E_{\mathit{DBT}}}(39)$

These four boundaries separate the Δν_K vs. $I_{E_{\mathit{inj}}}$ state space into five regions as shown in Figure 7.

In Figure 7A is shown this state space for the default inhibitory-excitatory coupling value $g_{IE}/{\widehat{g}}_{IE}=1$. As ν_K increases all these boundaries shift to lower $I_{E_{\mathit{inj}}}$ values. This contraction of the firing and bistable regions upon increase in extracelluar potassium is readily illustrated by the hysteretic dynamics in response to a symmetric ramp function, computed from integrating the ODEs as before, shown in Figure 7B. Note also that the transition times observed in the hysteresis correspond well with the analytically computed boundaries, which in turn assume steady state inputs and dynamics. Note also that at this modestly increased potassium level - equivalent to increasing from [K]_o = 3.5 mM to [K]_o = 6.2 mM - that the region of bistability is very small.

In Figure 7C is shown this state space for the decreased inhibitory-excitatory coupling value $g_{IE}/{\widehat{g}}_{IE}=0.5$, with accompanying hysteresis curves for ramp input shown in Figure 7D. Although at nominal extracellular potassium the input current needed to drive the network into DB is only modestly lower (about 18% lower) than for the nominal coupling, at the elevated potassium level is lowered by another factor of about 30% compared to the nominal coupling at elevated potassium.

The state space as a function of $g_{IE}/{\widehat{g}}_{IE}$ vs. $I_{E_{\mathit{inj}}}$ is shown in Figure 8 for values of Δν_K = 0, 15 mV. Here again the boundaries between behaviors are smooth, monotonically increasing functions. Higher ν_K leads to smaller areas of the dynamical regions F, BS, and RAE. Beyond the point where the RAE boundary (in green) crosses the DB boundary, the E mass goes into DB right after bistablity. Hence, we shift the DB (in red) boundary to the RAE boundary in the figures.

FIGURE 8

FIGURE 8. Shift in analytical Bounds as a function of coupling constant: The figures show the shift in dynamical region boundaries as the ratio $\frac{g_{IE}}{{\widehat{g}}_{IE}}$ is varied from 0.1 to 10, for Δν_k = 0 mV in (A), and Δν_k = 15 mV in (B). The grey horizontal lines depict the two coupling constants for which the region boundaries are shown in Figures 7A,C. The dotted blue and red lines in both the figures represent firing onset and DB onset boundaries of the E mass, when it is not coupled to the I mass.

We note that because of the smooth monotonically increasing shape of these boundaries, that the choice of our nominal value for ${\widehat{g}}_{IE}$ and comparison with a halved value would have lead to the same pattern of decreases in region space as discussed in the previous paragraphs.

4 Discussion and Conclusion

In the present work our aim was to establish and demonstrate a method to develop computational models of neuronal ensembles - neural masses - that stem from biophysical and coupling characteristics of membrane dynamics resolved cellular models of neurons, whose states and output could quantitatively be mapped across modeling scales, whose bifurcations realistically included depolarization block, and whose dynamics realistically bidirectionally couple to the neurons’ environment.

The developed neural mass models are parametrized through simple mathematical functions, show physiologically interpretable behaviors and dynamical transitions from one state to another, as a function of key parameters of neural environment.

In this work, we have demonstrated this development process for both an inhibitory neural mass, based on the Wang-Buzsáki (Wang and Buzsáki, 1996) (WB), and an accommodating excitatory neuron (SEAN), based on a simplification of the Pinsky-Rinzel (Pinsky and Rinzel, 1994) model. The masses are characterized in terms of time-averaged outputs of the corresponding neuron models, and represent the mean dynamics of populations of these neurons.

These mass descriptions can serve as plug-in replacements for existing NM elements used in other models, in that they produce firing rates as a function of synaptic input from other masses.

These masses also receive as input slow variables in the form of extracellular potassium, for both, and accommodation state for SEAN, and as output provide the coupling to those slow variables in the form of average membrane potential and transmembrane in our models, instead of using a sigmoidal function as the base relation between the input and firing rate, we use a fit for actual firing rates measured from ODE implementations of the original compartment models with each’s firing range, and explicitly the thresholds for firing and depolarization block as a function of the slow variables.

Inclusion of depolarization block has not been included in the classic NMMs (Freeman, 1972a; Freeman, 1972b; Wilson and Cowan, 1972; Lopes da Silva et al., 1974; Freeman, 1975; Jansen and Rit, 1995; Cona et al., 2014) before. Depolarization block was introduced into continuum Wilson-Cowen dynamics (Meijer et al., 2015) in a roughly ad-hoc fashion, but not dependent on [K]_o.

The explicit parametrization of depolarization block enables these new implementation of neural masses to express the underlying mechanisms of depolarization blocks and runaway excitation (McCormick and Contreras, 2001) needed for networks built from these elements to express realistic transitions to seizure and spreading depolarization. The sensitivity of both firing rates responses and DB transitions to extracellular potassium concentration then embodies a critical components that are known to play roles in epileptic dynamics and demonstrates the method to incorporate related volume transmitted signaling (Agnati et al., 1992).

Because these mass elements are built to parameterize these from the original cellular models, these variables map directly to measurable variables in brain. We investigated the dynamics of a simple two element network formed from coupling an excitatory and an inhibitory mass in which input was delivered only to the excitatory mass. This network itself displays key features that underlie epilepsy. In particular, as a function of input to the network it transitions in a hysteretic manner from firing to runaway excitation to depolarization block. Importantly, as extracellular potassium is increased, the input threshold to unstable dynamics is decreased. More importantly, as the inhibitory coupling to the excitatory network is decreased, this sensitivity to increases in extracellular potassium is greatly enhanced. This latter observation is consistent with the observation that epileptic networks are not seizing all the time, but are more susceptible to seizures and spreading depolarizations.

4.1 Limitations of the Model

The dynamics of the masses we developed are computationally efficient and can be used as direct plug-ins in existing NMM networks, yet their properties (firing rate vs input and environment) are closer to realistic.

This manuscript was built from two well known and established Hodgkin-Huxley style compartment models—that by Wang-Buzsáki (Wang and Buzsáki, 1996) and a reduced version of the Pinsky-Rinzel model (Pinsky and Rinzel, 1994) of a pyramidal neuron. This was done not because these models are the most accurate, but because they are recognizable and understood.

With this choice comes the limitations of these models. For example, we achieve firing rates from our SEAN model that far exceed the maximum firing rate of read pyramical neurons, which is a feature also noted in (Pinsky and Rinzel, 1994) for their soma-only reduction. In addition, by reducing this to a single compartment model we have eliminated burst-firing from this model as it is more complex.

We note that both the WB and SEAN models follow canonical Type 1 transition behavior, with roughly square-root firing rate behavior following. In our modeling, we have only parametrized that behavior and the subsequent transition to depolarization block. We leave for future work parametrization of the host of different transitions that have been articulated for example in (Izhikevich, 2000). Such efforts will further need to deal with the relationship between within burst dynamics and rate of axonally transmitted action potentials and resulting post-synaptic current generation. Likewise, a method to deal with hysteretic firing onset transitions such as observed with the simplest type 1 neurons (Izhikevich, 2000) may be needed.

For the current generated NMMs, we have parametrized the average response of a population of homogeneous neuron types formed with identical detailed parameters. Even when subdivided into common cell types, neurons in real biology will have a distribution of shapes, sizes and even ion channel or pump densities. These different within cell differences will lead to changes in the detailed input/output dynamics that we have parametrized in our models. Given a distribution of such parameters, it is straightforward to build the average mass response based on an average of the parameterized responses given such cell-parameter variations. As long as such cell-parameter variations do not change the firing rate from the square-root shape observed, these will translate simply to shifts in the positions of ISS_F and ISS_DB and the maximum firing rate as function of ν_K. As illustrated in Supplementary Figure S4, such averages primarily smooth out the distinct transitions at firing onset and DB offset. We do not expect that such changes will substantially alter the hysteretic network dynamics of firing, runaway excitation, and depolarization block. We anticipate the hardest part of incorporating the heterogenetity into such models is to establish and justify what distribution of cell parameters are realistic and should be used.

4.2 Future Directions

An interesting future work with these neural masses will include linking them to spatial networks of extracellular space that include tracking and diffusion of potassium as well as glial buffering and cell swelling. These elements should reveal components that spread seizure and SD events as well as express normal physiological rhythms. Through the procedure used to include other slow variables including intracellular sodium concentration, synaptic vesicle reserve, and oxygen dependent ATP production for driving the ionic pumps, it will also be possible to investigate the role of mutated channel dynamics at the network level.

Data Availability Statement

All data was generated within the publicly available code contained within the article.

Author Contributions

BJG conceived and designed the study. RT implemented the methodology. BJG and RT carried out the formal analysis. Both the authors wrote and reviewed the manuscript.

Funding

BJG efforts were supported through NIH grant R01 EB014641. RT was financially supported for this project by overseas research fellowship of IITGN. RT was partially funded by the Center of Advanced Systems Understanding (CASUS), which is financed by Germany’s Federal Ministry of Education and Research (BMBF) and by the Saxon Ministry for Science, Culture, and Tourism (SMWK) with tax funds on the basis of the budget approved by the Saxon State Parliament.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

RT is grateful to Pennsylvania State University (PSU) and Indian Institute of Technology Gandhinagar (IITGN) for facilitating this research work.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnetp.2022.911090/full#supplementary-material

References

Aarabi, A., and He, B. (2014). Seizure Prediction in Hippocampal and Neocortical Epilepsy Using a Model-Based Approach. Clin. Neurophysiol. 125, 930–940. doi:10.1016/j.clinph.2013.10.051

PubMed Abstract | CrossRef Full Text | Google Scholar

Agnati, L. F., Bjelke, B., and Fuxe, K. (1992). Volume Transmission in the Brain. Am. Sci. 80, 362–373.

Google Scholar

Aiba, I., and Noebels, J. L. (2015). Spreading Depolarization in the Brainstem Mediates Sudden Cardiorespiratory Arrest in Mouse SUDEP Models. Sci. Transl. Med. 7. doi:10.1126/scitranslmed.aaa4050

PubMed Abstract | CrossRef Full Text | Google Scholar

Bahari, F., Ssentongo, P., Liu, J., Kimbugwe, J., Curay, C., Schiff, S. J., et al. (2018). Seizure-associated Spreading Depression Is a Major Feature of Ictal Events in Two Animal Models of Chronic Epilepsy. bioRxiv. [Preprint]. doi:10.1101/455519

CrossRef Full Text | Google Scholar

Basu, I., Crocker, B., Farnes, K., Robertson, M. M., Paulk, A. C., Vallejo, D. I., et al. (2018). A Neural Mass Model to Predict Electrical Stimulation Evoked Responses in Human and Non-human Primate Brain. J. Neural Eng. 15, 066012. doi:10.1088/1741-2552/aae136

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhattacharya, B. S. (2013). Implementing the Cellular Mechanisms of Synaptic Transmission in a Neural Mass Model of the Thalamo-Cortical Circuitry. Front. Comput. Neurosci. 7, 1–17. doi:10.3389/fncom.2013.00081

PubMed Abstract | CrossRef Full Text | Google Scholar

Brennan, K. C., and Pietrobon, D. (2018). A Systems Neuroscience Approach to Migraine. Neuron 97, 1004–1021. doi:10.1016/j.neuron.2018.01.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Brown, T. H., and Johnston, D. (1983). Voltage-clamp Analysis of Mossy Fiber Synaptic Input to Hippocampal Neurons. J. Neurophysiology 50, 487–507. doi:10.1152/jn.1983.50.2.487

PubMed Abstract | CrossRef Full Text | Google Scholar

Butler, C. R., Graham, K. S., Hodges, J. R., Kapur, N., Wardlaw, J. M., and Zeman, A. Z. J. (2007). The Syndrome of Transient Epileptic Amnesia. Ann. Neurol. 61, 587–598. doi:10.1002/ana.21111

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, B., Choi, H., Hirsch, L. J., Katz, A., Legge, A., Buchsbaum, R., et al. (2017). Psychiatric and Behavioral Side Effects of Antiepileptic Drugs in Adults with Epilepsy. Epilepsy & Behav. 76, 24–31. doi:10.1016/j.yebeh.2017.08.039

CrossRef Full Text | Google Scholar

Cona, F., Lacanna, M., and Ursino, M. (2014). A Thalamo-Cortical Neural Mass Model for the Simulation of Brain Rhythms during Sleep. J. Comput. Neurosci. 37, 125–148. doi:10.1007/s10827-013-0493-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Corcoran, R., and Thompson, P. (1992). Memory Failure in Epilepsy: Retrospective Reports and Prospective Recordings. Seizure 1, 37–42. doi:10.1016/1059-1311(92)90053-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Costa, M. S., Born, J., Claussen, J. C., and Martinetz, T. (2016). Modeling the Effect of Sleep Regulation on a Neural Mass Model. J. Comput. Neurosci. 41, 15–28. doi:10.1007/s10827-016-0602-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Deschle, N., Ignacio Gossn, J., Tewarie, P., Schelter, B., and Daffertshofer, A. (2021). On the Validity of Neural Mass Models. Front. Comput. Neurosci. 118. doi:10.3389/fncom.2020.581040

CrossRef Full Text | Google Scholar

Diniz Behn, C. G., and Booth, V. (2010). Simulating Microinjection Experiments in a Novel Model of the Rat Sleep-Wake Regulatory Network. J. Neurophysiology 103, 1937–1953. doi:10.1152/jn.00795.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Fleshner, M., Booth, V., Forger, D. B., and Diniz Behn, C. G. (2011). Circadian Regulation of Sleep-Wake Behaviour in Nocturnal Rats Requires Multiple Signals from Suprachiasmatic Nucleus. Phil. Trans. R. Soc. A 369, 3855–3883. doi:10.1098/rsta.2011.0085

PubMed Abstract | CrossRef Full Text | Google Scholar

Freeman, W. J. (1972). Waves, Pulses, and the Theory of Neural Masses. Prog. Theor. Biol. 2, 1–10. doi:10.1016/b978-0-12-543102-6.50010-8

CrossRef Full Text | Google Scholar

Freeman, W. J. (1972). Linear Analysis of the Dynamics of Neural Masses. Annu. Rev. Biophys. Bioeng. 1, 225–256. doi:10.1146/annurev.bb.01.060172.001301

PubMed Abstract | CrossRef Full Text | Google Scholar

Freeman, W. J. (1975). Mass Action in the Nervous System, Vol. 2004. New York: Citeseer.

Google Scholar

Izhikevich, E. M. (2000). Neural Excitability, Spiking and Bursting. Int. J. Bifurc. Chaos 10, 1171–1266. doi:10.1142/S0218127400000840

CrossRef Full Text | Google Scholar

Jansen, B. H., and Rit, V. G. (1995). Electroencephalogram and Visual Evoked Potential Generation in a Mathematical Model of Coupled Cortical Columns. Biol. Cybern. 73, 357–366. doi:10.1007/BF00199471

PubMed Abstract | CrossRef Full Text | Google Scholar

Jansen, B. H., Zouridakis, G., and Brandt, M. E. (1993). A Neurophysiologically-Based Mathematical Model of Flash Visual Evoked Potentials. Biol. Cybern. 68, 275–283. doi:10.1007/BF00224863

PubMed Abstract | CrossRef Full Text | Google Scholar

Kager, H., Wadman, W. J., and Somjen, G. G. (2000). Simulated Seizures and Spreading Depression in a Neuron Model Incorporating Interstitial Space and Ion Concentrations. J. Neurophysiology 84, 495–512. doi:10.1152/jn.2000.84.1.495

PubMed Abstract | CrossRef Full Text | Google Scholar

Kameneva, T., Ying, T., Guo, B., and Freestone, D. R. (2017). Neural Mass Models as a Tool to Investigate Neural Dynamics during Seizures. J. Comput. Neurosci. 42, 203–215. doi:10.1007/s10827-017-0636-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Köksal Ersöz, E., Modolo, J., Bartolomei, F., and Wendling, F. (2020). Neural Mass Modeling of Slow-Fast Dynamics of Seizure Initiation and Abortion. PLoS Comput. Biol. 16, e1008430. doi:10.1371/journal.pcbi.1008430

CrossRef Full Text | Google Scholar

Lanzone, J., Ricci, L., Assenza, G., Ulivi, M., Di Lazzaro, V., and Tombini, M. (2018). Transient Epileptic and Global Amnesia: Real-Life Differential Diagnosis. Epilepsy & Behav. 88, 205–211. doi:10.1016/j.yebeh.2018.07.015

CrossRef Full Text | Google Scholar

Leao, A. A. P. (1944). Spreading Depression of Activity in the Cerebral Cortex. J. Neurophysiology 7, 359–390. doi:10.1152/jn.1944.7.6.359

CrossRef Full Text | Google Scholar

Leo, A. A. P., and Morison, R. S. (1945). Propagation of Spreading Cortical Depression. J. Neurophysiology 8, 33–45. doi:10.1152/jn.1945.8.1.33

CrossRef Full Text | Google Scholar

Leo, A. A. P. (1944). Pial Circulation and Spreading Depression of Activity in the Cerebral Cortex. J. Neurophysiology 7, 391–396. doi:10.1152/jn.1944.7.6.391

CrossRef Full Text | Google Scholar

Liu, C., Zhou, C., Wang, J., Fietkiewicz, C., and Loparo, K. A. (2021). Delayed Feedback-Based Suppression of Pathological Oscillations in a Neural Mass Model. IEEE Trans. Cybern. 51, 5046–5056. doi:10.1109/tcyb.2019.2923317

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, C., Zhu, Y., Liu, F., Wang, J., Li, H., Deng, B., et al. (2017). Neural Mass Models Describing Possible Origin of the Excessive Beta Oscillations Correlated with Parkinsonian State. Neural Netw. 88, 65–73. doi:10.1016/j.neunet.2017.01.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, F., Wang, J., Liu, C., Li, H., Deng, B., Fietkiewicz, C., et al. (2016). A Neural Mass Model of Basal Ganglia Nuclei Simulates Pathological Beta Rhythm in Parkinson's Disease. Chaos 26, 123113. doi:10.1063/1.4972200

PubMed Abstract | CrossRef Full Text | Google Scholar

Loonen, I. C. M., Jansen, N. A., Cain, S. M., Schenke, M., Voskuyl, R. A., Yung, A. C., et al. (2019). Brainstem Spreading Depolarization and Cortical Dynamics during Fatal Seizures in Cacna1a S218L Mice. Brain 142, 412–425. doi:10.1093/brain/awy325

PubMed Abstract | CrossRef Full Text | Google Scholar

Lopes da Silva, F. H., Hoeks, A., Smits, H., and Zetterberg, L. H. (1974). Model of Brain Rhythmic Activity. Kybernetik 15, 27–37. doi:10.1007/bf00270757

PubMed Abstract | CrossRef Full Text | Google Scholar

López-Cuevas, A., Castillo-Toledo, B., Medina-Ceja, L., and Ventura-Mejía, C. (2015). State and Parameter Estimation of a Neural Mass Model from Electrophysiological Signals during the Status Epilepticus. NeuroImage 113, 374–386. doi:10.1016/j.neuroimage.2015.02.059

PubMed Abstract | CrossRef Full Text | Google Scholar

Löscher, W., and Schmidt, D. (2011). Modern Antiepileptic Drug Development Has Failed to Deliver: Ways Out of the Current Dilemma. Epilepsia 52, 657–678. doi:10.1111/j.1528-1167.2011.03024.x

PubMed Abstract | CrossRef Full Text | Google Scholar

McCormick, D. A., and Contreras, D. (2001). On the Cellular and Network Bases of Epileptic Seizures. Annu. Rev. Physiol. 63, 815–846. doi:10.1146/annurev.physiol.63.1.815

PubMed Abstract | CrossRef Full Text | Google Scholar

Meador, K. J., and Loring, D. W. (2016). Developmental Effects of Antiepileptic Drugs and the Need for Improved Regulations. Neurology 86, 297–306. doi:10.1212/WNL.0000000000002119

PubMed Abstract | CrossRef Full Text | Google Scholar

Meijer, H. G., Eissa, T. L., Kiewiet, B., Neuman, J. F., Schevon, C. A., Emerson, R. G., et al. (2015). Modeling Focal Epileptic Activity in the Wilson-cowan Model with Depolarization Block. J. Math. Neurosci. 5, 7–17. doi:10.1186/s13408-015-0019-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Noebels, J. L. (2019). Brainstem Spreading Depolarization: Rapid Descent into the Shadow of SUDEP. Brain a J. Neurology 142, 231–233. doi:10.1093/brain/awy356

PubMed Abstract | CrossRef Full Text | Google Scholar

Pietrobon, D., and Moskowitz, M. A. (2014). Chaos and Commotion in the Wake of Cortical Spreading Depression and Spreading Depolarizations. Nat. Rev. Neurosci. 15, 379–393. doi:10.1038/nrn3770

PubMed Abstract | CrossRef Full Text | Google Scholar

Pinsky, P. F., and Rinzel, J. (1994). Intrinsic and Network Rhythmogenesis in a Reduced Traub Model for CA3 Neurons. J. Comput. Neurosci. 1, 39–60. doi:10.1007/bf00962717

PubMed Abstract | CrossRef Full Text | Google Scholar

Poh, M. Z., Loddenkemper, T., Reinsberger, C., Swenson, N. C., Goyal, S., Madsen, J. R., et al. (2012). Autonomic Changes with Seizures Correlate with Postictal EEG Suppression. Neurology 78, 1868–1876. doi:10.1212/WNL.0b013e318258f7f1

PubMed Abstract | CrossRef Full Text | Google Scholar

Rajakulendran, S., and Nashef, L. (2015). Postictal Generalized EEG Suppression and SUDEP: A Review. J. Clin. neurophysiology 32, 14–20. doi:10.1097/WNP.0000000000000147

PubMed Abstract | CrossRef Full Text | Google Scholar

Rall, W., Burke, R. E., Smith, T. G., Nelson, P. G., and Frank, K. (1967). Location of Synapses and Mechanisms for the Monosynaptic in Motoneurons Possible. J. Neurophysiol. 30, 1169–1193. doi:10.1152/jn.1967.30.5.1169

PubMed Abstract | CrossRef Full Text | Google Scholar

Rodrigues, S., Chizhov, A. V., Marten, F., and Terry, J. R. (2010). Mappings between a Macroscopic Neural-Mass Model and a Reduced Conductance-Based Model. Biol. Cybern. 102, 361–371. doi:10.1007/s00422-010-0372-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Ryvlin, P., Nashef, L., Lhatoo, S. D., Bateman, L. M., Bird, J., Bleasel, A., et al. (2013). Incidence and Mechanisms of Cardiorespiratory Arrests in Epilepsy Monitoring Units (MORTEMUS): A Retrospective Study. Lancet Neurology 12, 966–977. doi:10.1016/S1474-4422(13)70214-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Schellenberger Costa, M., Weigenand, A., Ngo, H.-V. V., Marshall, L., Born, J., Martinetz, T., et al. (2016). A Thalamocortical Neural Mass Model of the EEG during NREM Sleep and its Response to Auditory Stimulation. PLoS Comput. Biol. 12, e1005022–20. doi:10.1371/journal.pcbi.1005022

PubMed Abstract | CrossRef Full Text | Google Scholar

Segneri, M., Bi, H., Olmi, S., and Torcini, A. (2020). Theta-Nested Gamma Oscillations in Next Generation Neural Mass Models. Front. Comput. Neurosci. 14. doi:10.3389/fncom.2020.00047

PubMed Abstract | CrossRef Full Text | Google Scholar

Somjen, G. G. (2004). Ions in the Brain: Normal Function, Seizures, and Stroke. Oxford, New York: Oxford University Press.

Google Scholar

Ssentongo, P., Robuccio, A. E., Thuku, G., Sim, D. G., Nabi, A., Bahari, F., et al. (2017). A Murine Model to Study Epilepsy and SUDEP Induced by Malaria Infection. Sci. Rep. 7, 1–14. doi:10.1038/srep43652

PubMed Abstract | CrossRef Full Text | Google Scholar

Stefanescu, R. A., and Jirsa, V. K. (2008). A Low Dimensional Description of Globally Coupled Heterogeneous Neural Networks of Excitatory and Inhibitory Neurons. PLoS Comput. Biol. 4, e1000219. doi:10.1371/journal.pcbi.1000219

PubMed Abstract | CrossRef Full Text | Google Scholar

Svalheim, S., Sveberg, L., Mochol, M., and Taubøll, E. (2015). Interactions between Antiepileptic Drugs and Hormones. Seizure 28, 12–17. doi:10.1016/j.seizure.2015.02.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Thurman, D. J., Logroscino, G., Beghi, E., Hauser, W. A., Hesdorffer, D. C., Newton, C. R., et al. (2017). The Burden of Premature Mortality of Epilepsy in High-Income Countries: A Systematic Review from the Mortality Task Force of the International League against Epilepsy. Epilepsia 58, 17–26. doi:10.1111/epi.13604

PubMed Abstract | CrossRef Full Text | Google Scholar

Tolner, E. A., Chen, S.-P., and Eikermann-Haerter, K. (2019). Current Understanding of Cortical Structure and Function in Migraine. Cephalalgia 39, 1683–1699. doi:10.1177/0333102419840643

PubMed Abstract | CrossRef Full Text | Google Scholar

Traub, R. D., and Miles, R. (1991). Neuronal Networks of the hippocampus, Vol. 777. Cambridge, New York: Cambridge University Press.

PubMed Abstract | Google Scholar

Wang, X.-J., and Buzsáki, G. (1996). Gamma Oscillation by Synaptic Inhibition in a Hippocampal Interneuronal Network Model. J. Neurosci. 16, 6402–6413. doi:10.1523/jneurosci.16-20-06402.1996

PubMed Abstract | CrossRef Full Text | Google Scholar

Wei, Y., Ullah, G., and Schiff, S. J. (2014). Unification of Neuronal Spikes, Seizures, and Spreading Depression. J. Neurosci. 34, 11733–11743. doi:10.1523/jneurosci.0516-14.2014

PubMed Abstract | CrossRef Full Text | Google Scholar

Wendling, F., Bartolomei, F., Bellanger, J. J., and Chauvel, P. (2002). Epileptic Fast Activity Can Be Explained by a Model of Impaired Gabaergic Dendritic Inhibition. Eur. J. Neurosci. 15, 1499–1508. doi:10.1046/j.1460-9568.2002.01985.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Wendling, F. (2008). Computational Models of Epileptic Activity: a Bridge between Observation and Pathophysiological Interpretation. Expert Rev. Neurother. 8, 889–896. doi:10.1586/14737175.8.6.889

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilson, H. R., and Cowan, J. D. (1972). Excitatory and Inhibitory Interactions in Localized Populations of Model Neurons. Biophysical J. 12, 1–24. doi:10.1016/s0006-3495(72)86068-5

CrossRef Full Text | Google Scholar