2.1 Membrane potential dynamics
The pyramidal neuron is described by the modified Hodgkin-Huxley model[35–38]. The membrane potential \(V\) of neuron are as follows[37]:
$$C\dot{V}={I}_{ext}-{I}_{Na}-{I}_{K}-{I}_{Cl}$$
$${I}_{Na}={G}_{Na}{m}^{\text{3}}h\left(V- {E}_{Na}\right)+{G}_{Nal}(V-{E}_{Na})$$
$${I}_{K}={G}_{K}{n}^{\text{4}}\left(V- {E}_{K}\right)+{G}_{Kl}(V-{E}_{K})$$
$${I}_{Cl}={G}_{ClL}(V-{E}_{Cl})$$
$$\dot{q}={\alpha }_{q}\left(\text{1}-q\right)-{\beta }_{q}q, q=m,n,h$$
1
where \({I}_{ext}\) represents external input current, \({I}_{Na}\) and \({I}_{K}\) are sodium and potassium ion currents respectively and containing their respective leakage currents, and \({I}_{Cl}\) is chloride leakage current. \({G}_{Na}\) and \({G}_{K}\) represent the channel conductance of \({Na}^{+}\) and \({K}^{+}\). \({G}_{Nal}\), \({G}_{Kl}\) and \({G}_{ClL}\) represent the conductance of the three leak currents. \(m\) and \(n\) describe the sodium and potassium activation degree, \(h\) describes the sodium inactivation, and all vary between 0 and 1. \({\alpha }_{q}\) and \({\beta }_{q}\) (\(q=m,n,h\)) describe the opening and closing probability of ion channels, these parameters are described as follows[39, 40]:
$${\alpha }_{m}=\text{0.32}(\text{54}+V)/\left\{\text{1}-exp\left[-(V+\text{54})/\text{4}\right]\right\}$$
$${\beta }_{m}=\text{0.28}(\text{27}+V)/\left\{exp\left[(V+\text{27})/\text{5}\right]-\text{1}\right\}$$
$${\alpha }_{h}=\text{0.128}exp\left[-(V+\text{50})/\text{18}\right]$$
$${\beta }_{h}=\text{4}/\left\{\text{1}+exp\left[-(V+\text{27})/\text{5}\right]\right\}$$
$${\alpha }_{n}=\text{0.032}(\text{52}+V)/\left\{\text{1}-exp\left[-(V+\text{52})/\text{5}\right]\right\}$$
$${\beta }_{n}=\text{0.5}exp\left[-(V+\text{57})/\text{40}\right]$$
2
In equation 1, \({E}_{Na}\), \({E}_{K}\), and \({E}_{Cl}\) represent the equilibrium potential of \({Na}^{+}\), \({K}^{+}\) and \({Cl}^{-}\) that dependent on Nernst equations, as follows:
$${E}_{Na}=\text{26.64}\text{l}\text{n}\left({\left[{Na}^{+}\right]}_{o}/{\left[{Na}^{+}\right]}_{i}\right)$$
$${E}_{K}=\text{26.64}\text{l}\text{n}\left({\left[{K}^{+}\right]}_{o}/{\left[{K}^{+}\right]}_{i}\right)$$
$${E}_{Cl}=\text{26.64}\text{l}\text{n}\left({\left[{Cl}^{-}\right]}_{i}/{\left[{Cl}^{-}\right]}_{o}\right)$$
3
where \({\left[s\right]}_{o}\) and \({\left[s\right]}_{i}\) (\({\text{s}={Na}^{+},K}^{+},{Cl}^{-}\)) are ion concentration inside and outside the neuron respectively, and the ion concentrations are changing with neuronal activity except chloride concentration.
2.2 Ion concentration dynamics
In the model, ions concentration are changed follow neuronal activity. \({\left[{K}^{+}\right]}_{o}\) is influenced by \({K}^{+}\) currents (\({I}_{K}\)), \({Na}^{+}\)-\({K}^{+}\) pumps (\({I}_{pump}\) from neuron, \({I}_{gliapump}\) from glia), diffusion of \({K}^{+}\) from interstitial volume (\({I}_{diff}\)), and uptake by glial cell (\({I}_{glia}\))[35, 37], the specific form of \({\left[{K}^{+}\right]}_{o}\) is as follows:
$$\dot{{\left[{K}^{+}\right]}_{o}}=\text{0.0445}\beta {I}_{K}-\text{2}\beta {I}_{pump}-{I}_{diff}-{I}_{glia}-\text{2}{I}_{gliapump}$$
4
where \(\beta\) is the ratio of the intracellular volume to the outer volume. We supposed the quantity of \({Na}^{+}\) into the cell is equal to \({K}^{+}\) out of the cell, so \({\left[{K}^{+}\right]}_{i}\) can be described as:
$${\left[{K}^{+}\right]}_{i}=\text{140}\text{m}\text{M}+(\text{18}\text{m}\text{M}-{\left[{Na}^{+}\right]}_{i})$$
5
\({\left[{Na}^{+}\right]}_{i}\) is determined by \({Na}^{+}\) currents (\({I}_{Na}\)) and \({Na}^{+}\)-\({K}^{+}\) pumps, and \({\left[{Na}^{+}\right]}_{o}\) are updated based on \({\left[{Na}^{+}\right]}_{i}\), the specific form of \({[Na}^{+}]\) are as follows:
$$\dot{{\left[{Na}^{+}\right]}_{i}}=-\text{0.0445}{I}_{Na}-\text{3}{I}_{pump}$$
$${\left[{Na}^{+}\right]}_{o}=\text{144}\text{m}\text{M}-{\beta }({\left[{Na}^{+}\right]}_{i}-\text{18}\text{m}\text{M})$$
6
The \({Na}^{+}\)-\({K}^{+}\) pumps (\({I}_{pump}\) from neuron, \({I}_{gliapump}\) from glia) are modeled by:
$${I}_{pump}=\frac{\rho }{\text{1}+\text{e}\text{x}\text{p}[\left(\text{25}-{\left[{Na}^{+}\right]}_{i}\right)/\text{3}]}\times \frac{\text{1}}{\text{1}+\text{e}\text{x}\text{p}(\text{5.5}-{\left[{K}^{+}\right]}_{o})}$$
$${I}_{gliapump}=\frac{\text{1}}{\text{3}}\frac{\rho }{1+\text{e}\text{x}\text{p}[\left(\text{25}-{\left[{Na}^{+}\right]}_{gi}\right)/\text{3}]}\times \frac{1}{\text{1}+\text{e}\text{x}\text{p}(\text{5.5}-{\left[{K}^{+}\right]}_{o})}$$
7
where
$$\rho ={\rho }_{max}/(\text{1}+\text{e}\text{x}\text{p}[(\text{20}-{\left[{O}_{2}\right]}_{o})/\text{3}\left]\right)$$
8
where \({\rho }\) is \({Na}^{+}\)-\({K}^{+}\) pump rate which is updated based on extracellular oxygen concentration (\({\left[{O}_{2}\right]}_{o}\)).
The ability of glial cell handling \({\left[{K}^{+}\right]}_{o}\) is modeled by:
$${I}_{glia}={G}_{glia}/(\text{1}+\text{e}\text{x}\text{p}[\left(\text{18}-{\left[{K}^{+}\right]}_{o}\right)/\text{2.5}\left]\right)$$
9
where \({G}_{glia}\) is the glial uptake rate of \({\left[{K}^{+}\right]}_{o}\). The diffusion of \({K}^{+}\) from interstitial volume (\({I}_{diff}\)) is modeled by:
$${I}_{diff}={\epsilon }_{k}({\left[{K}^{+}\right]}_{o}-{\left[{K}^{+}\right]}_{ves})$$
10
where \({\epsilon }_{k}\) is \({K}^{+}\) diffusion constant, and \({\left[{K}^{+}\right]}_{ves}\) represents \({[K}^{+}]\) in blood vessel.
2.3 Oxygen concentration dynamics
In the work, we introduced and modified the blood vessel model[41], so that it can simulate blood oxygen metabolism in the state of ischemic stroke. The specific form of the blood vessel model are as follows:
$$BFlow=\eta {BFlow}_{ves}$$
$$\dot{OEF}=({OEF}_{ave}-{OEF}_{var}\text{t}\text{a}\text{n}\text{h}(\frac{BFlow-\text{0.7}}{\text{0.4}})-OEF)/{\tau }_{OEF}$$
$${\left[{O}_{2}\right]}_{ves}=OEF*BFlow\text{*}{[O}_{ini}]$$
11
where \(BFlow\) represents blood flow in blood vessel, \({\eta }\) describes the blockage degree of the blood vessel (\({\eta }\)=1 when the blood vessel in healthy state and \({\eta }\)=0 when the blood vessel is completely concluded), and \({\text{B}\text{F}\text{l}\text{o}\text{w}}_{ves}\) represents maximum blood flow in healthy blood vessel (\({\text{B}\text{F}\text{l}\text{o}\text{w}}_{ves}\)=1). \(OEF\) is oxygen extraction factor, \({OEF}_{ave}\) is normal oxygen extraction factor, \({OEF}_{var}\) represents the maximum possible variation of \(OEF\), and \({\tau }_{OEF}\) is an adaptation time that \(OEF\) adapt to \(\text{B}\text{F}\text{l}\text{o}\text{w}\) changes. \({\left[{O}_{2}\right]}_{ves}\) represents the concentration of oxygen supplied to the neuron by blood vessel, Figure 1 describes the oxygen concentration \({\left[{O}_{2}\right]}_{ves}\) varies with the blockage degree of the blood vessel \(\eta\), and \(\left[{O}_{ini}\right]\) represents the normal oxygen concentration in healthy blood vessel. The extracellular oxygen concentration (\({\left[{O}_{2}\right]}_{o}\)) is modeled by:
$$\dot{{\left[{O}_{\text{2}}\right]}_{O}}=-\alpha \left({I}_{pump}+{I}_{gliapump}\right)+{\epsilon }_{o}({\left[{O}_{\text{2}}\right]}_{ves}-{\left[{O}_{\text{2}}\right]}_{o})$$
12
where \(\alpha\) is a conversion factor from changing current to oxygen concentration[37]. \({\epsilon }_{o}\) is diffusion rate of oxygen.
2.4 Coupled population model
To study the propagation of epilepsy, we built a coupled population model based on the simple network model and consider the astrocyte feedback, including 150 pyramidal neurons, 150 astrocytes and a blood vessel. In the model, pyramidal neurons are connected by synapses, and astrocytes are connected by gap junctions, pyramidal neurons and astrocyte cells are connected with the bidirectional feedback, and equal concentration of oxygen are transported to every neurons and astrocytes by blood vessel. The schematic of model is shown in Figure 2. In this model, the membrane potential \(V\) of neuron is changed as follows:
$$C\dot{V}={I}_{Na}-{I}_{K}-{I}_{Cl}-{I}_{syn}-{I}_{as}$$
13
where
$${I}_{syn}={g}_{se}(V-{V}_{e})s$$
14
where \({I}_{syn}\) represents synaptic current, \({g}_{se}\) describes the activation level of glutamate receptors on the postsynaptic neuron, and \({V}_{e}\) is reversal potential of synapse. The synaptic input \(s\) from previous pyramidal neuron is modeled by:
$${\left[T\right]}_{i-\text{1}}=\frac{\text{1}}{\text{1}+\text{exp}\left(-\frac{V+\text{9}}{\text{8}}\right)}$$
$${\dot{s}}_{i}=\text{5}{p}_{neu}{\left[T\right]}_{i-\text{1}}\left(\text{1}-{s}_{i}\right)-{s}_{i}$$
15
where \({\left[T\right]}_{i-1}\) represents the released neurotransmitter by the \(i-1\)th neuron, since pyramidal neurons are excitatory neurons, the neurotransmitter released is mainly glutamate, so \({\left[T\right]}_{i-1}\) refers to glutamate in the model. \({p}_{neu}\) represents the synaptic input intensity to postsynaptic neuron.
In coupled population model, \({\left[{K}^{+}\right]}_{o}\) also is influenced by the \({K}^{+}\) diffusion from neighboring neurons, the specific form of the lateral diffusion term is as follows:
$${I}_{lat}=\frac{{D}_{k}}{{\varDelta x}^{\text{2}}}({{\left[{K}^{+}\right]}_{o}}^{i-\text{1}}+{{\left[{K}^{+}\right]}_{o}}^{i}+{{\left[{K}^{+}\right]}_{o}}^{i+\text{1}}-\text{3}{{\left[{K}^{+}\right]}_{o}}^{i})$$
16
where \({D}_{k}\) is the potassium diffusion coefficient, and \(\varDelta x\) is the distance between two neurons. So the final manifestation of \({\left[{K}^{+}\right]}_{o}\) is as follows:
$$\dot{{\left[{K}^{+}\right]}_{o}}=\text{0.0445}\beta {I}_{K}-\text{2}\beta {I}_{pump}-{I}_{diff}-{I}_{glia}-\text{2}{I}_{gliapump}+{I}_{lat}$$
17
The binding neurotransmitters to receptors on adjacent astrocytes will cause the astrocytes to produce \({IP}_{3}\), and eventually cause the concentration of \({Ca}^{2+}\) increasing in astrocytes. In order to describe this process, we introduced the improved Li-Rinzel model[42–44] as follows:
$$[\dot{I{P}_{\text{3}}]}=\frac{\left({[IP}_{\text{3}}^{*}]-{[IP}_{\text{3}}]\right)}{{\tau }_{{ip}_{\text{3}}}}+{r}_{{ip}_{\text{3}}}{p}_{as}\left[T\right]+{k}_{g}({\left[{IP}_{\text{3}}\right]}_{i+\text{1}}+{\left[{IP}_{\text{3}}\right]}_{i-\text{1}}-\text{2}{\left[{IP}_{\text{3}}\right]}_{i})$$
$$\dot{{[Ca}^{\text{2}+}]}={J}_{chan}+{J}_{leak}-{J}_{pump}$$
$$\dot{q}={\alpha }_{q}\left(\text{1}-q\right)-{\beta }_{q}q$$
18
where
$${J}_{chan}={c}_{\text{1}}{V}_{\text{1}}{p}_{\infty }^{\text{3}}{n}_{\infty }^{\text{3}}{q}^{\text{3}}\left({\left[{Ca}^{\text{2}+}\right]}_{ER}-\left[{Ca}^{\text{2}+}\right]\right)$$
$${J}_{leak}={c}_{\text{1}}{V}_{\text{2}}\left({\left[{Ca}^{\text{2}+}\right]}_{ER}-\left[{Ca}^{\text{2}+}\right]\right)$$
$${J}_{pump}=\frac{{V}_{\text{3}}{\left[{Ca}^{\text{2}+}\right]}^{\text{2}}}{\left[{Ca}^{\text{2}+}\right]+{k}_{\text{3}}^{\text{2}}}$$
19
where
$${p}_{\infty }=\frac{\left[{IP}_{\text{3}}\right]}{\left[{IP}_{\text{3}}\right]+{d}_{1}}$$
$${n}_{\infty }=\frac{\left[{Ca}^{\text{2}+}\right]}{\left[{Ca}^{\text{2}+}\right]+{d}_{\text{5}}}$$
$${\alpha }_{q}={a}_{\text{2}}{d}_{\text{2}}\frac{\left[{IP}_{\text{3}}\right]+{d}_{\text{1}}}{\left[{IP}_{\text{3}}\right]+{d}_{\text{3}}}$$
$${\beta }_{q}={a}_{\text{2}}\left[{Ca}^{\text{2}+}\right]$$
$${\left[{Ca}^{\text{2}+}\right]}_{ER}=\frac{{c}_{\text{0}}-\left[{Ca}^{\text{2}+}\right]}{{c}_{\text{1}}}$$
20
where \({[IP}_{3}]\) is the \({IP}_{3}\) concentration in astrocyte cytosolic, \({[IP}_{3}^{*}]\) is the balanced concentration\(\text{o}\text{f} {IP}_{3}\), \(\text{a}\text{n}\text{d} {p}_{as}\) represents the intensity of synaptic input to neighboring astrocyte, \({k}_{g}\) refers to the coupling coefficient of astrocyte gap junctions. \(\left[{Ca}^{2+}\right]\) is the \({Ca}^{2+}\) concentration in astrocyte cytosolic, \({J}_{chan}\) and \({J}_{leak}\) refer to \({Ca}^{2+}\) flow and leakage flow from endoplasmic reticulum (ER) to cytosol, and \({J}_{pump}\) refers to the pump flow from ER to cytosol.\({\left[{Ca}^{2+}\right]}_{ER}\) is the \({Ca}^{2+}\) concentration in astrocyte ER.
When \({[Ca}^{2+}]\) exceeds a threshold, astrocytes release gliotransmitters into synapses to regulate neuronal activity. According this process, a dynamical variable \(f\) is introduced to picture the astrocyte feedback[44] which has the following form:
$$\dot{f}=\frac{-f}{{\tau }_{{Ca}^{\text{2}+}}}+(\text{1}-f)\kappa {\Phi }(\left[{Ca}^{\text{2}+}\right]-{\left[{Ca}^{\text{2}+}\right]}_{th})$$
21
and combine the strength of the astrocyte feedback \({\lambda }\) which describes the activation of the metabotropic glutamate receptors on the presynaptic neuron, a final equation about astrocyte feedback current is described as follow:
$${I}_{as}={\lambda }f$$
22
The bidirectional feedback between astrocytes and neighbor neurons is achieved by “tripartite synapse” which consist of pre- and postsynaptic neurons and astrocyte[45–47]. It’s reported that gliotransmitters inhibit the release of neurotransmitters when gliotransmitters act on presynaptic neurons and inhibit synaptic activity [48, 49]. So in the model, the astrocyte feedback is inhibitory.