Absence epilepsy, a generalized non-convulsive epilepsy, can be observed from the electroencephalogram (EEG). A sudden, brief loss of consciousness with bilaterally synchronous 2–4 Hz SWDs is the main feature of absence seizures, which mainly affects children and young adolescents (Crunelli and Leresche, 2002). Panayiotopoulos (2008) found that absence seizures propagated across the whole cortex and happened without any prediction. The possible involvement of the thalamus during epileptic seizures has been shown in the human recordings (Jasper and Kershman, 1941; Prevett et al., 1995). After that, a series of electrophysiological recordings in experimental models showed the complexity of their dynamical characteristics and corresponding biophysics mechanisms of so-called corticothalamic system (Avoli and Gloor, 1981; Steriade et al., 1993a,b; Seidenbecher et al., 1998). A large number of results on the generalized non-convulsive epilepsy have shown that the generation of SWDs is due to the abnormal interactions between the cerebral cortex and thalamus (Marescaux and Vergnes, 1995; Timofeev and Steriade, 2004). And some studies of computational modeling also have been utilized to investigate more deep insights into the possible generation mechanisms of SWDs in the corticothalamic system (Breakspear et al., 2006; Marten et al., 2009; Chen et al., 2014).

Mean-field modeling can be useful to study the dynamical behaviors of collective neurons. It has been increasingly used to investigate the evolution behaviors of epileptic seizures, especially for absence seizures. For examples, Robinson et al. (2002) focused on understanding the transition from the steady state to the typical SWD state by means of the constructed mean-field model. In addition to the classical SWD activities, Marten et al. (2009) demonstrated much more dynamical phenomena including poly-spike and wave, wave spike or even no discernible spike-wave onset in the recorded EEG data. Furthermore, Destexhe (1998) suggested that the SWDs could arise from thalamocortical loops in which the corticothalamic feedback indirectly evoked GABA_B-mediated inhibition in the thalamus. Based on the dynamical bifurcations of the mean-field model in the brain, Breakspear et al. (2006) explained critical features of the generalized epilepsies. Recently, Chen et al. (2014) discovered a computational evidence for the bilateral control of absence seizures by simulating the typical SWD morphology, which once again highlights the critical role of mean-field model.

There are many different methods in the treatment of epilepsy, and pharmacotherapy is currently the primary modality (Löscher, 2002). Another techniques such as surgical resection of epileptogenic foci and neuromodulatory techniques have also been used. Especially for the case of intractable epilepsy, surgical intervention has been developed rapidly because of technological advancement. In general, 86 percent of patients were satisfied with their postoperative results, which were supported by previous statistical work (Macrodimitris et al., 2011). Even when the patients suffered from the seizures again after the surgery, 90 percent of patients stated that the surgical treatment was the right decision (Reid et al., 2004). Although surgical intervention can be well accepted by epilepsy patients, complications of surgery are still the potential risks. By contrast, DBS sends electrical pulses to the implanted electrodes that are placed in specific areas of the brain to control seizures. This method was found to be effective to modulate the epileptic brain activity in some animal experiments, especially acting on the thalamus (Mina et al., 2013). Nanobashvili et al. (2003) showed that the stimulation of thalamic reticular nucleus (TRN) suppressed the development of seizures that was taken from Male Wistar rats. In addition, many similar studies have indicated that the high-frequency stimulation of TRN could obtain an anti-epileptogenic effect and interrupt abnormal electroencephalogram recordings in rats (Pantoja-Jiménez et al., 2014). Hence, the above animal experiments suggest a new approach for seizure control in treating epilepsy. DBS might become a very attractive and meaningful method to treat epilepsy patients, because it is relative safe and reversible (Huang and van Luijtelaar, 2016). Thus, compared with the other neurosurgical treatments, DBS needs the complicated implantation procedures and advanced techniques to localize the suitable target in the brain (Fuentes et al., 2009; Lega et al., 2010).

According to the above findings, we develop a mean-field model for the epileptic absence seizures through adding the DBS to TRN, and analyze its neurophysiological effects on epileptic cortical dynamics. By means of dynamical analysis, we explore how the slow kinetics of GABA_B receptors on TRN induce the typical absence epilepsy activities. Furthermore, we also study how the typical SWD oscillation activities are modulated and terminated via changing the range of relevant stimulation parameters (i.e., frequency, amplitude and positive input duration). These results in this paper can further highlight the critical roles of thalamus stimulation in eliminating absence seizures, and can provide useful therapeutic method for the brain disease.

In order to simulate the EEG activity during absence seizures, we build a model to describe the population dynamics in two distinct brain regions, namely the thalamus and the cerebral cortex. The network connection schematic of this model has been presented in Figure 1, which includes four neural populations: e = excitatory pyramidal neurons (EPN); i = inhibitory interneurons (IIN); r = thalamus reticular nucleus (TRN); s = specific relay nuclei (SRN) (Marten et al., 2009; Chen et al., 2014). There are three types of neural projections in the network. For the sake of simplicity, we distinguish them in different types of lines and heads. Excitatory projections are mediated by glutamate, which are denoted by the red lines with arrow heads, while inhibitory projections are mediated by GABA_A and GABA_B receptors, which are represented by the blue solid and dashed lines with round heads, respectively. Additionally, DBS denotes the deep brain stimulation, which is a neuromodulator technique.

Based on the important roles of cortex and thalamus during epileptic seizures (Marescaux and Vergnes, 1995; Timofeev and Steriade, 2004; Breakspear et al., 2006), we modify a mean-field model (Marten et al., 2009; van Albada and Robinson, 2009; Chen et al., 2014) by removing the basal ganglia and adding the DBS to TRN to explore the macroscopic dynamics of neural populations. For each neural population, the mean firing rate Q_a is described by an increasing sigmoid function of the mean membrane potential V_a, and given by (Robinson et al., 1997, 1998), (1)Qa(r,t)=F[Va(r,t)]=Qamax1+exp[−π3(Va(r,t)−θa)σ], where a ∈ A = {e, i, r, s} refer to different neural populations, Qamax represents the maximum firing rate, θ_a is the mean threshold potential, σ is the standard deviation of firing thresholds, and r denotes the spatial position. With an average firing rate Q_a, the neural population fires action potentials when V_a exceeds the threshold θ_a. As the sigmoid shape of Q_a is physiologically crucial for the model, the average firing rate cannot exceed the maximum firing rate Qamax. The change of the average membrane potential for neurons in type a can be modeled as van Albada et al. (2009), (2)DαβVa(r,t)=∑b∈Avab·ϕb(r,t), (3)Dαβ=1αβ[∂2∂t2+(α+β)∂∂t+αβ], where the differential operator D_αβ represents the synaptic and dendritic filtering of incoming signals, α and β denote the decay and rise rate of the cell-body potential, respectively (van Albada et al., 2009). ϕ_b(r, t) is the incoming pulse rate and v_ab is connection strength from the neural population of type b to type a. For the sake of simplicity, we ignore the transmission delay among some neural populations including the delay between cortex and thalamus in our model. Nevertheless, some animal experiments have shown that the slow kinetics of inhibitory GABA_B receptors on TRN can result in the generation of absence seizures (Destexhe, 1998; Chen et al., 2014). Hence, to investigate the effect of slow kinetics of inhibitory GABA_B receptors on absence seizures, we introduce two connections ϕ_b(r, t) and ϕ_b(r, t−τ) from TRN to SRN, which are used to denote the inhibitory GABA_A and GABA_B synapses. The postsynaptic currents of GABA_B functions develop on a slower time scales than the currents of GABA_A since the GABA_B functions via second messenger processes. Therefore, the GABA_B-mediated inhibitory projection will be delayed with respect to GABA_A when the thalamic neurons fire.

In the cerebral cortex, the propagation of cortical excitatory axonal field ϕ_e can be well-approximated by a damped wave equation (Jirsa and Haken, 1996), (4)1γe2[∂2∂t2+2γe∂∂t+γe2−ve2∇2]ϕe(r,t)=Qe(r,t), where ∇² is the Laplacian operator and Q_e is the excitatory firing rates. γe=vere is the temporal damping rate of pulses, consisting of the characteristic range r_e of axons and the mean conduction velocity v_e of excitatory neurons. Except for the axons of cortical excitatory pyramidal neurons, other axons are too short to provide wave propagation on the relevant scales. Here, ϕ_c = F(V_c)(c = i, r, s). Moreover, the absence seizures, typical generalized seizures, the dynamical activities are regarded as global brain activities. Thus a reasonable hypothesis can be proposed that the spatial activities are unified. For this purpose, the spatial derivative can be ignored and then we can set ∇² = 0 in Equation (4). Accordingly, the propagation effect of ϕ_e can be finally given by, (5)1γe2[∂2∂t2+2γe∂∂t+γe2]ϕe(r,t)=Qe(r,t).

For cortical neuronal populations, We assume that intracortical connectivities are proportional to the numbers of synapses involved (Robinson et al., 2002; van Albada et al., 2009), which implies V_i = V_e, and thus Q_a = F(V_a) implies Q_i = Q_e. Although the model has been simplified in many aspects, it doesn't affect the results of this work.

In the end, we will introduce the DBS, which can be described as a periodic step function with the following form (Fan et al., 2016), (6)DBS(t)=a×H(sin(2πtρ))(1−H(sin(2π(t+δ))ρ)).

Where H is the Heaviside bi-value step function, which satisfies that H(x) = 1 if x > 0 and H(x) = 0 if x ≤ 0. a and ρ is the amplitude and period, respectively. δ is the duration of positive input which determines duty circle of positive pulse. In order to simulate the non-regular property of DBS, the instantaneous frequency f of incoming pulses can be expressed as 1T according to the previous literature (So et al., 2012). A series of monophasic current pulses have been shown as in Figure 2.

When the DBS is applied on the TRN, we will rewrite above Equations (1–6) in the first-order form for all neural populations. The final mathematical description of the model is given as follows, (7)dϕe(t)dt=ϕ˙e(t), (8)dϕ˙e(t)dt=γe2[−ϕe(t)+F(Ve(t))]−2γeϕ˙e(t), (9)dVe(t)dt=V˙e(t), (10)dV˙e(t)dt=αβ[veeϕe(t)+veiF(Vi(t))+vesF(Vs(t))−Ve(t)]                    − (α+β)V˙e(t), (11)dVr(t)dt=V˙r(t), (12)dV˙r(t)dt=αβ[vreϕe(t)+vrsF(Vs(t))−Vr(t)]− (α+β)V˙r(t)                   + DBS(t), (13)dVs(t)dt=V˙s(t), (14)dV˙s(t)dt=αβ[vseϕe(t)+vsrAF(Vr(t))+vsrBF(Vr(t−τ))      − (Vs(t)−ϕn)]                     − (α+β)V˙s(t).

Here ϕ_n is the constant nonspecific subthalamic input onto SRN, τ is GABA_B delay from TRN to SRN.

Unless otherwise mentioned, the values of parameters we used have been shown in Table 1. The first nine different parameters are approximately estimated from physiological experiments, and they have been applied in some literatures (Rennie et al., 2000; Robinson et al., 2001; Zhao and Robinson, 2015). The range of delay τ is set as [0, 180] ms, which is estimated from a number of different subjects with absence epilepsy (Marten et al., 2009). However v_ab (a, b ∈ A = {e, i, r, s}) are less well known from physiology (Robinson et al., 2002), which leads to a quite wide range of v_ab. In addition, we will change some critical parameters within certain ranges (i.e., v_sr, τ, f, α, and δ) to acquire different dynamical features and explore the possible effects of the parameters on controlling absence seizures.

All of the numerical simulations are performed in the MATLAB (MathWorks, USA) environment. Similar to the previous works (Chen et al., 2014, 2015), the delay differential equations in the model are solved by the standard fourth-order Runge-Kutta method. The fixed temporal resolution of numerical integration is 0.05 ms. All the simulations are performed up to 15 s and the data from 5 to 15 s are used for statistical analysis. We don't consider the effect of noise. What's more, the step length of integrate is small enough to ensure the numerical accuracy and reliability in the simulations.

Taking advantage of the mean-field macroscopic model which associates with the thalamus and the cerebral cortex, we have studied how the deep brain stimulation eliminates the absence seizures. The obtain results have shown that the appearance of some transition states, definitely including absence seizures in our model, results from the change of the coupling strength of TRN-SRN pathway and the slow synaptic kinetics of GABA_B receptors on TRN. We have also found that a relatively long delay τ is required to assure the occurrence of SWDs for stronger −v_sr. Combining with some existing animal experiments, we choose the TRN as a suitable stimulation target in this paper. The detail dynamics and qualitative analysis have shown that adding the suitable DBS to TRN contributes to an obvious reduction of the SWD oscillation region. We have further investigated how the corresponding stimulation parameters including frequency, amplitude and positive input duration suppress the SWDs under certain conditions. The results have indicated that suitable stimulations can eliminate absence seizures as the DBS is applied to TRN. However, it is well known that a human's brain is quite different from a rodent's and clinical trials can't be replaced by animal models. Thus, the information from animal models only can guide clinician to explore how the treatment is applied to human. Hence, we hope that our computational results can provide a theoretical guidance for the treatment of real absence epilepsy patients.

ZW and QW designed and performed the research as well as wrote the paper.