The inhibition mechanism of epilepsy disease in a computational model

Abstract

BACKGROUND:

The mechanism of prevention and treatment of epilepsy is a hot issue in theoretical research.

OBJECTIVE:

In this paper, we studied the control mechanism of the generalized spike-and-wave discharges (GSWD) by different types of external electrical stimulation acting on the subthalamic nucleus (STN) in a computational model.

METHODS:

Firstly, we analyzed the pathological mechanism of seizures, which were induced by different parameters in the thalamocortical (TC) circuit. Then, a voltage V was exerted in the STN. At last, we used the sine wave and square wave current stimulation in the STN.

RESULTS:

We found that seizures can be inhibited by tuning stimulus intensity into suitable range, and the direction of adjustment depended on the size of the parameter. We observed that the seizure can also be inhibited by tuning different parameters in current.

CONCLUSIONS:

Different inhibition mechanisms can be explained in this model, which may provide theoretical evidences for selecting the optimal treatment scheme in the clinical.

1. Introduction

Absence epilepsy (AE) is the most common type of primary generalized epilepsy in children [1]. The clinical features of AE are convulsive unconsciousness attacks, and the generalized spike-and-wave discharges (GSWD) (2–4 Hz) appear in the background of normal electroencephalogram (EEG) [2]. Recently, many studies based on the functional MRI localization (fMRI) of EEG synchronous have found that the thalamus structure and cerebral cortex closely relate to the onset of the GSWD [3, 4], and the theoretical mechanism of which has been explored in many dynamic models [5, 6]. Amiri et al. found that disruption of the signaling function of astrocytes may induce seizure activities in a modified TC model [5]. Haghighi et al. studied the transition mechanism between ictal and interictal states based on dynamic analysis in a TC network [6]. Shojaee et al. observed that the change of coupling strength in TC network can be taken as detector to predict epilepsy seizures [7]. These results may provide reference for the prevention of seizures.

Antiepileptic drug therapy is a common method to relieve seizures [8], drugs control the regulation mechanism inside the body. However, the treatment and mitigation of many intractable seizures must resort to surgical brain stimulation technique, such as the deep brain stimulation (DBS). The DBS has been proved to be very effective in relieving epilepsy symptoms and the subthalamic nucleus (STN) is a common stimulus area [9, 10, 11]. Handforth et al. found that subthalamic DBS as an adjunct treatment for refractory epilepsy can partly control seizures [11]. Usui et al. shown that unilateral STN DBS has significant effect in inhibiting limbic seizures [12]. Ke et al. pointed out that electrical stimulation on dorsal STN can control melanocortinergic signaling in astrocytes [13]. Kase et al. indicated that the STN is involved in the rhythm maintenance and STN DBS can inhibit the GSWD in AE by tuning frequency [14]. Wille et al. shown that the STN may be an effective stimulus target for progressive myoclonic epilepsy, especially for less impaired patients [15]. Prabhu et al. pointed out that STN DBS could be used as an experimental therapy for motor epilepsy [16]. Recently, Wang et al. reported a favorable outcome of STN DBS for focal myoclonic seizures [17]. And, the STN DBS has been proved to be successful in treating many other nervous system diseases, such as Parkinson’s disease [18, 19] and dystonia [20]. The effect of the DBS for controlling drug-resistant seizures in childhood can refer to a systematic review in the reference [21]. However, the mechanism of the DBS is still unclear and theoretical research evidences for epilepsy are rare now.

In this paper, we used a computational model [22, 23] to explore the inhibition mechanism of the GSWD by different electrical stimulation acting in the STN. The results may provide theoretical evidences for experimental studies in inhibiting of AE, and can inspire control mechanism research in other types of epilepsy.

2. Network model and computing method

2.1 Mathematical mean-field model

The network connection structure is shown in Fig. 1 [22, 23]. Different nerve nuclei are abbreviated as follows: $\varsigma =$ STN; $p_1=$ substantia nigra pars reticulata (SNr); $p_2=$ external globus palladium (GPe); $d_1=$ striatal D1 neurons (SD1); $d_2=$ striatal D2 neurons (SD2); e $=$ excitatory pyramidal neurons (EPN); i $=$ inhibitory interneurons (IIN); r $=$ thalamic reticular nucleus (TRN); s $=$ specific relay nuclei (SRN). The globus pallidus internal (GPi) and SNr are taken as a single population in this model [22, 23]. In this connection structure diagram, arrowheads represent excitatory projections mediated by glutamate; The black lines with round heads are inhibitory connections adjusted through $\gamma$ -aminobutyric acid-A ( ${\text{𝐺𝐴𝐵𝐴}}_A$ ) and $\gamma$ -aminobutyric acid-B ( ${\text{𝐺𝐴𝐵𝐴}}_B$ ) receptors. In this model, the basal ganglia (BG) sends inhibitory outputs to the TC network by the pathway “SNr $\to$ SRN” and “SNr $\to$ TRN”, and if the stimulation is exerted in BG, the relieve of seizures may be realized by the regulation between these two pathways [22]. Therefore, we infer that seizures may be inhibited by employing the external stimulation in the STN in this model. The TC network mainly sends excitatory projections to the BG, i.e., the cortex exerts excitatory inputs to the striatum and STN, and the SRN has excitatory feedbacks to the striatum.

Figure 1.

The connection structure diagram of the basal ganglia-thalamocortical model (BGTC) [22, 23].

The coupling relation of the model is described by the following first-order mean-field differential equations, the specific biological significance of which has been discussed in previous research [22, 23, 24, 25, 26]:

${\displaystyle \frac{d{\varphi }_e(t)}{dt}}={\dot{\varphi }}_e(t)$ ${\displaystyle \frac{d{\dot{\varphi }}_e(t)}{dt}}={\gamma }_e^2[-{\varphi }_e(t)+F(V_e(t))]-2{\gamma }_e{\dot{\varphi }}_e(t)$ ${\displaystyle \frac{dX(t)}{dt}}=\dot{X}(t)$ $X(t)={[V_e(t),V_{d_1}(t),V_{d_2}(t),V_{p_1}(t),V_{p_2}(t),V_{\zeta }(t),V_r(t),V_s(t)]}^T$ ${\displaystyle \frac{d{\dot{V}}_e(t)}{dt}}=\alpha \beta ({\nu }_{ee}{\varphi }_e+{\nu }_{ei}F(V_e(t))+{\nu }_{es}F(V_s(t))-V_e(t))-(\alpha +\beta ){\dot{V}}_e(t)$ ${\displaystyle \frac{d{\dot{V}}_{d_1}(t)}{dt}}=\alpha \beta ({\nu }_{d_{1}e}{\varphi }_e+{\nu }_{d_1{d}_1}F(V_{d_1}(t))+{\nu }_{d_{1}s}F(V_s(t))-V_{d_1}(t))-(\alpha +\beta ){\dot{V}}_{d_1}(t)$ ${\displaystyle \frac{d{\dot{V}}_{d_2}(t)}{dt}}=\alpha \beta ({\nu }_{d_{2}e}{\varphi }_e+{\nu }_{d_2{d}_2}F(V_{d_2}(t))+{\nu }_{d_{2}s}F(V_s(t))-V_{d_2}(t))-(\alpha +\beta ){\dot{V}}_{d_2}(t)$ ${\displaystyle \frac{d{\dot{V}}_{p_1}(t)}{dt}}=\alpha \beta ({\nu }_{p_1{d}_1}F(V_{d_1}(t))+{\nu }_{p_1{p}_2}F(V_{p_2}(t))+{\nu }_{p_1\zeta }F(V_{\zeta }(t))-V_{p_1}(t))-(\alpha +\beta ){\dot{V}}_{p_1}(t)$ ${\displaystyle \frac{d{\dot{V}}_{p_2}(t)}{dt}}=\alpha \beta ({\nu }_{p_2{d}_2}F(V_{d_2}(t))+{\nu }_{p_2{p}_2}F(V_{p_2}(t))+{\nu }_{p_2\zeta }F(V_{\zeta }(t))-V_{p_2}(t))-(\alpha +\beta ){\dot{V}}_{p_2}(t)$ ${\displaystyle \frac{d{\dot{V}}_{\zeta }(t)}{dt}}=\alpha \beta ({\nu }_{\zeta e}{\varphi }_e+{\nu }_{\zeta p_2}F(V_{p_2}(t))-V_{\zeta }(t)+V)-(\alpha +\beta ){\dot{V}}_{\zeta }(t)$ ${\displaystyle \frac{d{\dot{V}}_r(t)}{dt}}=\alpha \beta ({\nu }_{re}{\varphi }_e+{\nu }_{r{p}_1}F(V_{p_1}(t))+{\nu }_{rs}F(V_s(t))-V_r(t))-(\alpha +\beta ){\dot{V}}_r(t)$ ${\displaystyle \frac{d{\dot{V}}_s(t)}{dt}}=\alpha \beta ({\nu }_{se}{\varphi }_e+{\nu }_{s{p}_1}F(V_{p_1}(t))+{\nu }_{sr}^{A}F(V_r)+{\nu }_{sr}^{B}F(V_r(T-\tau ))-V_s(t)+{\varphi }_n)$ $\mathrm{\hspace{1em}}-(\alpha +\beta ){\dot{V}}_s(t)$

2.2 Numerical simulation technology and data

All simulations were obtained by the Runge-Kutta algorithm and bifurcation analysis [22, 26]. The default data were derived from previous studies [22, 23, 24, 25, 26, 27, 28, 29], which were listed as follows: The unit in mV s: ${\nu }_{es}={\nu }_{ei}=$ 1.8, ${\nu }_{p_2\varsigma }=$ 0.45, ${\nu }_{d_{2}e}=$ 0.7, ${\nu }_{d_2{d}_2}={\nu }_{p_2{d}_2}=$ 0.3, ${\nu }_{d_1{d}_1}=$ 0.2, ${\varphi }_n=$ 2, ${\nu }_{ee}={\nu }_{d_{1}e}=$ 1, ${\nu }_{p_1{p}_2}=$ 0.03, ${\nu }_{\varsigma e}={\nu }_{d_{1}s}={\nu }_{p_1\varsigma }={\nu }_{p_1{d}_1}=$ 0.1, ${\nu }_{re}=$ 0.05, ${\nu }_{s{p}_1}={\nu }_{r{p}_1}=$ 0.035, ${\nu }_{\varsigma p_2}=$ 0.04, ${\nu }_{rs}=$ 0.5, ${\nu }_{p_2{p}_2}=$ 0.075, ${\nu }_{d_{2}s}=$ 0.05, ${\nu }_{se}=$ 2.2; The unit in mV: ${\theta }_{p_2}=$ 9, ${\theta }_{d_1},{\theta }_{d_2}=$ 19, ${\theta }_{p_1}={\theta }_{\varsigma }=$ 10, ${\theta }_r={\theta }_s={\theta }_e,{\theta }_i=$ 15, $\sigma =$ 6; The unit in Hz: $Q_{p_2}^{\max }=$ 300, $Q_{\varsigma }^{\max }=$ 500, $Q_{d_1}^{\max },Q_{d_2}^{\max }=$ 65, $Q_r^{\max }=Q_{p_1}^{\max }=Q_s^{\max }=Q_e^{\max },Q_i^{\max }=$ 250, ${\gamma }_e=$ 100; $\alpha =$ 50 s ${}^{-1}$ , $\beta =$ 200 s ${}^{-1}$ , $\tau =$ 50 ms.

Figure 2.

(a): A state bifurcation process induced by $-{\nu }_{ei}$ . (b)–(e): Four oscillation states obtained by setting different values of $-{\nu }_{ei}$ . We take $-{\nu }_{sr}=$ 1.5, ${\nu }_{re}=$ 0.06, ${\nu }_{d_{1}e}=$ 1.3, ${\nu }_{d_{1}s}=$ 0.2, ${\nu }_{d_{2}e}=$ 0.78, ${\nu }_{d_{2}s}=$ 0.055, $-{\nu }_{p_1{d}_1}=$ 0.15, $-{\nu }_{p_1{p}_2}=$ 0.05, $-{\nu }_{p_2{d}_2}=$ 0.37, ${\nu }_{p_2\varsigma }=$ 0.5, $-{\nu }_{\varsigma p_2}=$ 0.046, ${\nu }_{\varsigma e}=$ 0.3, ${\nu }_{p_1\varsigma }=$ 0.23, ${\varphi }_n=$ 2.3 and $\tau =$ 0.058 in simulations. (f)–(i): The control effects obtained in the plane ( ${\nu }_{se},V$ ), ( ${\nu }_{es},V$ ), ( ${\nu }_{rs},V$ ) and ( $T,V$ ), respectively. Seizures induced by ${\nu }_{se}$ , ${\nu }_{es}$ , ${\nu }_{rs}$ and T can be inhibited by V, as indicated by arrows. Here, we take $-{\nu }_{sr}=$ 1.

3. Main results

Coupling weights are main parameters in the model, which represent connection strengths between different populations and the change of them is easy to induce state transition in the TC network. Figure 2a is a typical state bifurcation process, which was induced by $-{\nu }_{ei}$ . ${\nu }_{ei}$ is the inhibitory coupling weight from the IIN to the EPN, the increase of which will inhibit the firing ability of the EPN. When the activation level (ACL) of the EPN is high, it will appear maximum firing state D with firing rate of 250 Hz, as shown in Fig. 2b. As $-{\nu }_{ei}$ increasing, the ACL of the EPN reduces, and the SWD state C appears, as shown in Fig. 2c. They are typical absence epileptiform activities with four poles in a period, similar waveforms are often observed in the EEG of absence epilepsy patients. The state B represents simple periodic oscillation with firing rate is lower than the state C (as shown in Fig. 2d). When $-{\nu }_{ei}$ increases to large enough, oscillatory activities disappear and the EPN was pushed into the low firing state A (as shown in Fig. 2e). Therefore, the onset mechanism of seizures can be understood in this model together with numerical simulation. Of course, in addition to the $-{\nu }_{ei}$ , there are also many other pathways in the thalamocortical network that may induce seizures, such as “EPN $\to$ SRN” [30], “SRN $\to$ EPN” [23], “SRN $\to$ TRN” [23] and delay [22, 30]. We infer that the dynamic mechanisms of other pathological factors in the TC network are similar in this model, and the control mechanism of which will be discussed in the follow.

Figure 2f–i show four control analysis diagrams in ( ${\nu }_{se},V$ ), ( ${\nu }_{es},V$ ), ( ${\nu }_{rs},V$ ) and ( $T,V$ ), respectively. They show that seizures induced by ${\nu }_{se}$ , ${\nu }_{es}$ , ${\nu }_{rs}$ and T in the thalamocortical system can also be inhibited by V, as indicated by arrows, which implies that the control method is effective in theory.

Figure 3.

(a): The sine wave current. Here, we set the parameter P in (1) to 30 s. (b): The control effect when the current Isin acts on STN, with P $=$ 30. (c): The square wave current. Here, we set the parameter P in (2) to be 3 s and D $=$ 2 s. (d): The control effect diagram obtained by adjusting the D. Here, we set A $=$ 50, T $=$ 40. (e): The state transition diagram obtained by tuning the T. Here, we set simulation calculation points in the interval [0.0375, 3] for T is 80, A $=$ 40, D $=$ 0.0006 and $-{\nu }_{sr}=$ 1.5. We take the resistance $R=1\mathrm{\Omega }$ .

In addition to uniform voltage stimulation, other types of variable – frequency currents, such as sine wave (1) and square wave (2) [31, 32] stimulation, are often used in clinical trials. Where A represents the amplitude, P is a parameter associated with the period of the current, and D is the effective duration of current in a period. Figure 3a describes the sine wave current. Figure 3b is the inhibition effect of the seizures by employing the sine wave current on the STN. Indicated by arrows, state C may be relieved via tuning the current amplitude A into suitable values. However, the control effect depends on $-{\nu }_{sr}$ . As $-{\nu }_{sr}$ is small, the excitatory projection from the SRN to the EPN is strong; thus, the ACL of the EPN is high and the EPN may only stay in the seizure state or the saturation state. Therefore, the seizure transfers into the saturation state by enhancing the A. With increasing of $-{\nu }_{sr}$ , the ACL of the EPN decreases gradually, and the state C, B and A may appear in the EPN. In this case, the seizure state may only be pushed into the simple firing state or the low firing state, indicated by the downward arrow. Figure 3c describes the square wave current. Figure 3d is the inhibition effect by employing the square wave current on the STN. As indicated by arrows, state C disappears by tuning the current parameter D into suitable values. Similar to Fig. 3b, the control effect also depends on $-{\nu }_{sr}$ . From Fig. 3c, we can infer that the larger D, the stronger the stimulation effect in a period. Thus, the inhibited mechanism in Fig. 3d may be similar to that in Fig. 3b. Figure 3e represents a state transition obtained by tuning the T, which implies that the seizure state disappears when T falls into appropriate ranges through periodically transferring with the low firing state. The mechanism of Fig. 3e may be difficult to explain by this model, but we infer that it depends on the periodic characteristic of the current itself.

$I_{\sin }=(A/2)*(\sin (2*\pi *t/P)+1)$(1) $I_{\text{𝑠𝑞𝑢𝑎𝑟𝑒}}=A*\text{ℎ𝑒𝑎𝑣𝑖𝑠𝑖𝑑𝑒}(\sin (2*\pi *t/P))*(1-\text{ℎ𝑒𝑎𝑣𝑖𝑠𝑖𝑑𝑒}(\sin (2*\pi *(t+D)/P)))$(2)

4. Conclusion

In this paper, we used a basal ganglia – thalamocortical computational model to study the control mechanism of absence epileptiform activities by different electrical stimulation acting in the STN. Firstly, we analyzed the pathological mechanism of seizures, which were caused by the coupling weight ${\nu }_{ei}$ in the TC loop. We found that the origin mechanism can be understood in this model. Then, we exerted a voltage in the STN and turned the stimulus intensity. We found that seizures disappeared as the voltage strength falling into appropriate intervals, and the direction of adjustment depended on the size of pathological parameters (such as, ${\nu }_{se},{\nu }_{es},{\nu }_{rs}$ and $\tau$ ). Lastly, we used the sine wave and square wave current in the STN and observed that the seizure can also be inhibited by changing the amplitude, period and effective duration in current. The regulatory mechanisms of these control schemes are different, which can be well understood in this computational model. The results may provide theoretical support for selecting the optimal control scheme in treating of epilepsy.