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

Abstract

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

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 and ν_K, respectively, the change in it is expressed as Δν_K.

Here . [.]_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 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.

FIGURE 2

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

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:

where,

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 , , and the nominal Nernst potentials for these ions are . The Nernst potential for sodium ion is assumed to be constant, i.e., . 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:

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:

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:

where,

where ρ = 1.25 is the pump strength. mM, 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:

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: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:

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)):

In practice, the mNM parametrization is computed just using the fast dynamics I_K−AHP,fast, where we have chosen . 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: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 ) 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:

The hatted versions of the Nernst potentials ( and ) 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 , depolarization threshold , maximum , maximum , and maximum , as functions of for the I mass are:

Similarly, the firing threshold , depolarization threshold , maximum , maximum , and maximum functions (of ) that determine the threshold values for the E mass are:

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 , mean potassium flux functions (KF_I) for the I mass are:where, C1 = −60, , , and are the maximum of the respective quantities as determined for a given from Eq. 14. Q_I is the normalized total current that the mass receives and is expressed as:

Notice that, for the isolated mass this total current is just the injected current . 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:

The mean firing rate function (FR_E), mean membrane potential function , mean potassium flux functions (KF_E) for the E mass are:where, C1 = −65 and C2 = 0.1 and normalized total current (Q_E) is expressed as:

Here again, the total current 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).

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:

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:In this form, any postsynaptic current in neuron b from action potentials in i will be characterized as,

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),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:

Note that under this parametrization, under steady state conditions F_a(t) = F_a and 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

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)) . 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 is constant, and replaced with an overall coupling constant g_a,b.

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, , and ) at three different ν_K values. We also find out these three quantities (⟨FR⟩_mass, , and ), 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

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: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: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 and . 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

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 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 . This is better observed in Figure 6B where the firing rates are plotted as a function of , and the sign of 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

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 - Δν_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

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

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

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

These four boundaries separate the Δν_K vs. 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 . As ν_K increases all these boundaries shift to lower 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 , 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 vs. 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

We note that because of the smooth monotonically increasing shape of these boundaries, that the choice of our nominal value for 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.

References

• 1

  AarabiA.HeB. (2014). Seizure Prediction in Hippocampal and Neocortical Epilepsy Using a Model-Based Approach. Clin. Neurophysiol.125, 930–940. 10.1016/j.clinph.2013.10.051

  • CrossRef
  • Google Scholar

• 2

  AgnatiL. F.BjelkeB.FuxeK. (1992). Volume Transmission in the Brain. Am. Sci.80, 362–373.

  • Google Scholar

• 3

  AibaI.NoebelsJ. L. (2015). Spreading Depolarization in the Brainstem Mediates Sudden Cardiorespiratory Arrest in Mouse SUDEP Models. Sci. Transl. Med.7. 10.1126/scitranslmed.aaa4050

  • CrossRef
  • Google Scholar

• 4

  BahariF.SsentongoP.LiuJ.KimbugweJ.CurayC.SchiffS. J.et al (2018). Seizure-associated Spreading Depression Is a Major Feature of Ictal Events in Two Animal Models of Chronic Epilepsy. bioRxiv. [Preprint]. 10.1101/455519

  • CrossRef
  • Google Scholar

• 5

  BasuI.CrockerB.FarnesK.RobertsonM. M.PaulkA. C.VallejoD. 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. 10.1088/1741-2552/aae136

  • CrossRef
  • Google Scholar

• 6

  BhattacharyaB. 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. 10.3389/fncom.2013.00081

  • CrossRef
  • Google Scholar

• 7

  BrennanK. C.PietrobonD. (2018). A Systems Neuroscience Approach to Migraine. Neuron97, 1004–1021. 10.1016/j.neuron.2018.01.029

  • CrossRef
  • Google Scholar

• 8

  BrownT. H.JohnstonD. (1983). Voltage-clamp Analysis of Mossy Fiber Synaptic Input to Hippocampal Neurons. J. Neurophysiology50, 487–507. 10.1152/jn.1983.50.2.487

  • CrossRef
  • Google Scholar

• 9

  ButlerC. R.GrahamK. S.HodgesJ. R.KapurN.WardlawJ. M.ZemanA. Z. J. (2007). The Syndrome of Transient Epileptic Amnesia. Ann. Neurol.61, 587–598. 10.1002/ana.21111

  • CrossRef
  • Google Scholar

• 10

  ChenB.ChoiH.HirschL. J.KatzA.LeggeA.BuchsbaumR.et al (2017). Psychiatric and Behavioral Side Effects of Antiepileptic Drugs in Adults with Epilepsy. Epilepsy & Behav.76, 24–31. 10.1016/j.yebeh.2017.08.039

  • CrossRef
  • Google Scholar

• 11

  ConaF.LacannaM.UrsinoM. (2014). A Thalamo-Cortical Neural Mass Model for the Simulation of Brain Rhythms during Sleep. J. Comput. Neurosci.37, 125–148. 10.1007/s10827-013-0493-1

  • CrossRef
  • Google Scholar

• 12

  CorcoranR.ThompsonP. (1992). Memory Failure in Epilepsy: Retrospective Reports and Prospective Recordings. Seizure1, 37–42. 10.1016/1059-1311(92)90053-4

  • CrossRef
  • Google Scholar

• 13

  CostaM. S.BornJ.ClaussenJ. C.MartinetzT. (2016). Modeling the Effect of Sleep Regulation on a Neural Mass Model. J. Comput. Neurosci.41, 15–28. 10.1007/s10827-016-0602-z

  • CrossRef
  • Google Scholar

• 14

  DeschleN.Ignacio GossnJ.TewarieP.SchelterB.DaffertshoferA. (2021). On the Validity of Neural Mass Models. Front. Comput. Neurosci.118. 10.3389/fncom.2020.581040

  • CrossRef
  • Google Scholar

• 15

  Diniz BehnC. G.BoothV. (2010). Simulating Microinjection Experiments in a Novel Model of the Rat Sleep-Wake Regulatory Network. J. Neurophysiology103, 1937–1953. 10.1152/jn.00795.2009

  • CrossRef
  • Google Scholar

• 16

  FleshnerM.BoothV.ForgerD. B.Diniz BehnC. G. (2011). Circadian Regulation of Sleep-Wake Behaviour in Nocturnal Rats Requires Multiple Signals from Suprachiasmatic Nucleus. Phil. Trans. R. Soc. A369, 3855–3883. 10.1098/rsta.2011.0085

  • CrossRef
  • Google Scholar

• 17

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

  • CrossRef
  • Google Scholar

• 18

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

  • CrossRef
  • Google Scholar

• 19

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

  • Google Scholar

• 20

  IzhikevichE. M. (2000). Neural Excitability, Spiking and Bursting. Int. J. Bifurc. Chaos10, 1171–1266. 10.1142/S0218127400000840

  • CrossRef
  • Google Scholar

• 21

  JansenB. H.RitV. G. (1995). Electroencephalogram and Visual Evoked Potential Generation in a Mathematical Model of Coupled Cortical Columns. Biol. Cybern.73, 357–366. 10.1007/BF00199471

  • CrossRef
  • Google Scholar

• 22

  JansenB. H.ZouridakisG.BrandtM. E. (1993). A Neurophysiologically-Based Mathematical Model of Flash Visual Evoked Potentials. Biol. Cybern.68, 275–283. 10.1007/BF00224863

  • CrossRef
  • Google Scholar

• 23

  KagerH.WadmanW. J.SomjenG. G. (2000). Simulated Seizures and Spreading Depression in a Neuron Model Incorporating Interstitial Space and Ion Concentrations. J. Neurophysiology84, 495–512. 10.1152/jn.2000.84.1.495

  • CrossRef
  • Google Scholar

• 24

  KamenevaT.YingT.GuoB.FreestoneD. R. (2017). Neural Mass Models as a Tool to Investigate Neural Dynamics during Seizures. J. Comput. Neurosci.42, 203–215. 10.1007/s10827-017-0636-x

  • CrossRef
  • Google Scholar

• 25

  Köksal ErsözE.ModoloJ.BartolomeiF.WendlingF. (2020). Neural Mass Modeling of Slow-Fast Dynamics of Seizure Initiation and Abortion. PLoS Comput. Biol.16, e1008430. 10.1371/journal.pcbi.1008430

  • CrossRef
  • Google Scholar

• 26

  LanzoneJ.RicciL.AssenzaG.UliviM.Di LazzaroV.TombiniM. (2018). Transient Epileptic and Global Amnesia: Real-Life Differential Diagnosis. Epilepsy & Behav.88, 205–211. 10.1016/j.yebeh.2018.07.015

  • CrossRef
  • Google Scholar

• 27

  LeaoA. A. P. (1944). Spreading Depression of Activity in the Cerebral Cortex. J. Neurophysiology7, 359–390. 10.1152/jn.1944.7.6.359

  • CrossRef
  • Google Scholar

• 28

  LeoA. A. P.MorisonR. S. (1945). Propagation of Spreading Cortical Depression. J. Neurophysiology8, 33–45. 10.1152/jn.1945.8.1.33

  • CrossRef
  • Google Scholar

• 29

  LeoA. A. P. (1944). Pial Circulation and Spreading Depression of Activity in the Cerebral Cortex. J. Neurophysiology7, 391–396. 10.1152/jn.1944.7.6.391

  • CrossRef
  • Google Scholar

• 30

  LiuC.ZhouC.WangJ.FietkiewiczC.LoparoK. A. (2021). Delayed Feedback-Based Suppression of Pathological Oscillations in a Neural Mass Model. IEEE Trans. Cybern.51, 5046–5056. 10.1109/tcyb.2019.2923317

  • CrossRef
  • Google Scholar

• 31

  LiuC.ZhuY.LiuF.WangJ.LiH.DengB.et al (2017). Neural Mass Models Describing Possible Origin of the Excessive Beta Oscillations Correlated with Parkinsonian State. Neural Netw.88, 65–73. 10.1016/j.neunet.2017.01.011

  • CrossRef
  • Google Scholar

• 32

  LiuF.WangJ.LiuC.LiH.DengB.FietkiewiczC.et al (2016). A Neural Mass Model of Basal Ganglia Nuclei Simulates Pathological Beta Rhythm in Parkinson's Disease. Chaos26, 123113. 10.1063/1.4972200

  • CrossRef
  • Google Scholar

• 33

  LoonenI. C. M.JansenN. A.CainS. M.SchenkeM.VoskuylR. A.YungA. C.et al (2019). Brainstem Spreading Depolarization and Cortical Dynamics during Fatal Seizures in Cacna1a S218L Mice. Brain142, 412–425. 10.1093/brain/awy325

  • CrossRef
  • Google Scholar

• 34

  Lopes da SilvaF. H.HoeksA.SmitsH.ZetterbergL. H. (1974). Model of Brain Rhythmic Activity. Kybernetik15, 27–37. 10.1007/bf00270757

  • CrossRef
  • Google Scholar

• 35

  López-CuevasA.Castillo-ToledoB.Medina-CejaL.Ventura-MejíaC. (2015). State and Parameter Estimation of a Neural Mass Model from Electrophysiological Signals during the Status Epilepticus. NeuroImage113, 374–386. 10.1016/j.neuroimage.2015.02.059

  • CrossRef
  • Google Scholar

• 36

  LöscherW.SchmidtD. (2011). Modern Antiepileptic Drug Development Has Failed to Deliver: Ways Out of the Current Dilemma. Epilepsia52, 657–678. 10.1111/j.1528-1167.2011.03024.x

  • CrossRef
  • Google Scholar

• 37

  McCormickD. A.ContrerasD. (2001). On the Cellular and Network Bases of Epileptic Seizures. Annu. Rev. Physiol.63, 815–846. 10.1146/annurev.physiol.63.1.815

  • CrossRef
  • Google Scholar

• 38

  MeadorK. J.LoringD. W. (2016). Developmental Effects of Antiepileptic Drugs and the Need for Improved Regulations. Neurology86, 297–306. 10.1212/WNL.0000000000002119

  • CrossRef
  • Google Scholar

• 39

  MeijerH. G.EissaT. L.KiewietB.NeumanJ. F.SchevonC. A.EmersonR. G.et al (2015). Modeling Focal Epileptic Activity in the Wilson-cowan Model with Depolarization Block. J. Math. Neurosci.5, 7–17. 10.1186/s13408-015-0019-4

  • CrossRef
  • Google Scholar

• 40

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

  • CrossRef
  • Google Scholar

• 41

  PietrobonD.MoskowitzM. A. (2014). Chaos and Commotion in the Wake of Cortical Spreading Depression and Spreading Depolarizations. Nat. Rev. Neurosci.15, 379–393. 10.1038/nrn3770

  • CrossRef
  • Google Scholar

• 42

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

  • CrossRef
  • Google Scholar

• 43

  PohM. Z.LoddenkemperT.ReinsbergerC.SwensonN. C.GoyalS.MadsenJ. R.et al (2012). Autonomic Changes with Seizures Correlate with Postictal EEG Suppression. Neurology78, 1868–1876. 10.1212/WNL.0b013e318258f7f1

  • CrossRef
  • Google Scholar

• 44

  RajakulendranS.NashefL. (2015). Postictal Generalized EEG Suppression and SUDEP: A Review. J. Clin. neurophysiology32, 14–20. 10.1097/WNP.0000000000000147

  • CrossRef
  • Google Scholar

• 45

  RallW.BurkeR. E.SmithT. G.NelsonP. G.FrankK. (1967). Location of Synapses and Mechanisms for the Monosynaptic in Motoneurons Possible. J. Neurophysiol.30, 1169–1193. 10.1152/jn.1967.30.5.1169

  • CrossRef
  • Google Scholar

• 46

  RodriguesS.ChizhovA. V.MartenF.TerryJ. R. (2010). Mappings between a Macroscopic Neural-Mass Model and a Reduced Conductance-Based Model. Biol. Cybern.102, 361–371. 10.1007/s00422-010-0372-z

  • CrossRef
  • Google Scholar

• 47

  RyvlinP.NashefL.LhatooS. D.BatemanL. M.BirdJ.BleaselA.et al (2013). Incidence and Mechanisms of Cardiorespiratory Arrests in Epilepsy Monitoring Units (MORTEMUS): A Retrospective Study. Lancet Neurology12, 966–977. 10.1016/S1474-4422(13)70214-X

  • CrossRef
  • Google Scholar

• 48

  Schellenberger CostaM.WeigenandA.NgoH.-V. V.MarshallL.BornJ.MartinetzT.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. 10.1371/journal.pcbi.1005022

  • CrossRef
  • Google Scholar

• 49

  SegneriM.BiH.OlmiS.TorciniA. (2020). Theta-Nested Gamma Oscillations in Next Generation Neural Mass Models. Front. Comput. Neurosci.14. 10.3389/fncom.2020.00047

  • CrossRef
  • Google Scholar

• 50

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

  • Google Scholar

• 51

  SsentongoP.RobuccioA. E.ThukuG.SimD. G.NabiA.BahariF.et al (2017). A Murine Model to Study Epilepsy and SUDEP Induced by Malaria Infection. Sci. Rep.7, 1–14. 10.1038/srep43652

  • CrossRef
  • Google Scholar

• 52

  StefanescuR. A.JirsaV. K. (2008). A Low Dimensional Description of Globally Coupled Heterogeneous Neural Networks of Excitatory and Inhibitory Neurons. PLoS Comput. Biol.4, e1000219. 10.1371/journal.pcbi.1000219

  • CrossRef
  • Google Scholar

• 53

  SvalheimS.SvebergL.MocholM.TaubøllE. (2015). Interactions between Antiepileptic Drugs and Hormones. Seizure28, 12–17. 10.1016/j.seizure.2015.02.022

  • CrossRef
  • Google Scholar

• 54

  ThurmanD. J.LogroscinoG.BeghiE.HauserW. A.HesdorfferD. C.NewtonC. 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. Epilepsia58, 17–26. 10.1111/epi.13604

  • CrossRef
  • Google Scholar

• 55

  TolnerE. A.ChenS.-P.Eikermann-HaerterK. (2019). Current Understanding of Cortical Structure and Function in Migraine. Cephalalgia39, 1683–1699. 10.1177/0333102419840643

  • CrossRef
  • Google Scholar

• 56

  TraubR. D.MilesR. (1991). Neuronal Networks of the hippocampus, Vol. 777. Cambridge, New York: Cambridge University Press.

  • Google Scholar

• 57

  WangX.-J.BuzsákiG. (1996). Gamma Oscillation by Synaptic Inhibition in a Hippocampal Interneuronal Network Model. J. Neurosci.16, 6402–6413. 10.1523/jneurosci.16-20-06402.1996

  • CrossRef
  • Google Scholar

• 58

  WeiY.UllahG.SchiffS. J. (2014). Unification of Neuronal Spikes, Seizures, and Spreading Depression. J. Neurosci.34, 11733–11743. 10.1523/jneurosci.0516-14.2014

  • CrossRef
  • Google Scholar

• 59

  WendlingF.BartolomeiF.BellangerJ. J.ChauvelP. (2002). Epileptic Fast Activity Can Be Explained by a Model of Impaired Gabaergic Dendritic Inhibition. Eur. J. Neurosci.15, 1499–1508. 10.1046/j.1460-9568.2002.01985.x

  • CrossRef
  • Google Scholar

• 60

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

  • CrossRef
  • Google Scholar

• 61

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

  • CrossRef
  • Google Scholar