Eliminating Absence Seizures through the Deep Brain Stimulation to Thalamus Reticular Nucleus

Deep brain stimulation (DBS) can play a crucial role in the modulation of absence seizures, yet relevant biophysical mechanisms are not completely established. In this paper, on the basis of a biophysical mean-field model, we investigate a typical absence epilepsy activity by introducing slow kinetics of GABAB receptors on thalamus reticular nucleus (TRN). We find that the region of spike and slow-wave discharges (SWDs) can be reduced greatly when we add the DBS to TRN. Furthermore, we systematically explore how the corresponding stimulation parameters including frequency, amplitude and positive input duration suppress the SWDs under certain conditions. It is shown that the SWDs can be controlled as key stimulation parameters are suitably chosen. The results in this paper can be helpful for researchers to understand the thalamus stimulation in treating epilepsy patients, and provide theoretical basis for future experimental and clinical studies.


INTRODUCTION
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 nonconvulsive 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.

DESCRIPTION OF MODELS
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(Robinson et al., , 1998, where a ∈ A = {e, i, r, s} refer to different neural populations, Q max a 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 Q max a . The change of the average membrane potential for neurons in type a can be modeled as van Albada et al. (2009), 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), where ∇ 2 is the Laplacian operator and Q e is the excitatory firing rates. γ e = v e r e 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 ∇ 2 = 0 in Equation (4). Accordingly, the propagation effect of φ e can be finally given by, 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) 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), 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 1 T 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, 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(Chen et al., , 2015, the delay differential equations in the model are solved by the standard fourthorder 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

Generation of the Typical Absence Epilepsy Activities by Introducing Pathological Mechanisms
Recent studies have reported that the slow kinetics of GABA B receptors on TRN is an important pathological factor in destroying the normal oscillatory patterns of the corticothalamic system (Hosford et al., 1992;Destexhe, 1998;Timofeev and Steriade, 2004). Inspired by these findings, we will explore the impact of the slow kinetics of GABA B receptors on epileptic seizures by means of the reasonable model as given above.
Initially, we plot one-dimensional bifurcation diagrams as the inhibitory coupling strength −v sr and the GABA B delay τ are changed, respectively (see Figures 3a,b). In particular, the bifurcation diagrams are obtained by plotting the stable local minimum and maximum values of cortical excitatory axonal fields φ e . As shown in Figure 3a, there are four different dynamical states including the saturation region (I), the SWD oscillation region (II), the simple oscillation region (III) and the low firing region (IV) as −v sr is changed. For more clear vision, some typical time series are depicted in Figures 3c1-c4 with some fixed −v sr , respectively, which support the results of Figure 3a. It should be noted that the saturation state is a nonphysiological brain state, whose threshold is set as 250 Hz in our model. In addition, when the coupling strength −v sr is too strong, the TRN almost inhibits the firing of SRN. In this case, the model transfers into the low firing states and no oscillation phenomenon can be found. With decreasing the coupling strength −v sr , we find an approximate 3 Hz simple oscillation at −v sr ≈ 1.16 mVs, which is similar to oscillatory EEG signals at the onset of absence seizures in some patients (Rodrigues et al., 2005). And their modes just contain approximately 3 Hz sine-like oscillations and don't include any spikes. Hence, the simple oscillation actually may lead to a pathological system. And the state transition between the simple and SWD oscillation has been shown in Figure 3a. When the coupling strength −v sr is too weak, the inhibition from TRN hardly suppresses the firing of SRN. The strong excitatory from cortical pyramidal neurons makes firing of the SRN arrive a high level immediately. The high excitatory level promotes the firing of cortical neurons, and reaches the saturation states (region I). From Figure 3b, we can observe two types of oscillation patterns as the delay τ is changed: the SWD oscillation (region II) with multiple pairs of maximums and minimums and the simple oscillation (region III) with only one pair of maximums and minimums at the fixed delay parameter. Furthermore, as shown in Figure 3d, four different state regions are displayed in two-parameter space (−v sr , τ ), and the SWD oscillation state is represented by region II. The dominant frequency is estimated using the power spectral analysis, which is obtained from the time series of the excitatory axonal fields φ e via fast Fourier transform. Meanwhile, the maximum peak frequency Frontiers in Computational Neuroscience | www.frontiersin.org is defined as the dominant frequency of neural oscillation. In particular, Figure 3e shows the frequency range of epileptic seizures within 2-4 Hz. It should be noted that the region surrounded by the red dashed lines in Figure 3d represents the 2-4 Hz SWD oscillation region, i.e., absence seizures region (pathological region). We also find that the system can generate the 2-4 Hz SWD activities when the range of −v sr belongs to (0.47, 1.04) mVs and the parameter τ is larger than 40 ms as shown in Figure 3d. Hence, when the system is located in the region II within the red line, the decrease of τ can lead a transition from the SWD to simple oscillation state. However, both decreasing and increasing the inhibitory coupling strength −v sr can also suppress the SWDs. In addition, increasing the parameters −v sr and τ simultaneously can decrease the dominant frequency of neural oscillations only for the simple and SWD oscillation states (Figure 3e). Therefore, we obtain that the GABA B delay might has a greater impact on the dominant frequency than the coupling strength −v sr . So, a relatively long delay is required to generate the SWD oscillation activities for stronger coupling. From what we have mentioned above, we simulate four different dynamical states of the cerebral cortex via the slow kinetics of GABA B receptors and coupling strength −v sr of TRN-SRN pathway. In particular, the typical absence epilepsy state can be successful replicated. Throughout the whole process, we set τ = 50 ms, which pertains to the physiological range, contributing to the generation of 2-4 Hz SWD oscillation patterns.

Control of Absence Seizures by the DBS to TRN
DBS of the thalamus is a new and developing therapy, which is effective for the refractory epilepsy patients who are not suitable for pharmacotherapy or resective surgery (Valentín et al., 2013). Previous researchers utilized many animal and clinical experiments to confirm the stimulating efficacy for different brain structures, such as the subthalamic nucleus (Chabardès et al., 2002), the centromedian thalamic nucleus (Valentín et al., 2013), Anterior nucleus of the thalamus (Mirski et al., 1997), thalamic reticular nucleus (Nanobashvili et al., 2003) and the cortical regions (Morrell, 2011). The thalamus as a crucial nervous center is closely related to epileptic seizures. Here, we study the effect of the DBS to TRN on the generation and elimination of absence seizures. As shown in Figure 1, the TRN neurons send the inhibitory signals mediated by both GABA A and GABA B receptors to SRN. After the firing of SRN is inhibited by the GABA A receptors, the SRN neurons need some time to recover to their initial state. Meanwhile, the GABA B -induced inhibition will generate the delay τ . Therefore if the delay τ is longer than the recovery time, the SRN neurons can induce another firing peak. Hence the relatively large delay τ is easy to produce absence epilepsy activities. For a stronger −v sr , the GABA A -induced inhibition is also strong. In this situation, the SRN neurons need a longer time to restore their firing rate. Therefore, a relatively long τ is required to ensure the occurrence of SWDs for stronger −v sr . In addition, the DBS mainly affects the interactions among different neuron populations in our model, i.e., the range of delay τ in pathological region remains nearly unchanged. Based on above discussions, if we want to control absence seizures, the most straightforward method is to stimulate TRN because the TRN stimulation contributes to a stronger inhibitory effect.
Similarly, one-dimensional bifurcation diagrams are plotted with regard to the inhibitory coupling strength −v sr (Figure 4a) and the delay parameter τ (Figure 4b), respectively. As shown in Figure 4a, we find four different firing states as the parameter −v sr is changed, and their corresponding time series are depicted in Figures 4c1-c4. It is clear that the transitions in Figures 4a,b are similar to those in Figures 3a,b without DBS. However, further investigations find that the range of parameters −v sr and τ for different dynamical states are changed. Especially, the absence epilepsy state within the range of parameters −v sr and τ is reduced significantly. For more details, we perform the state analysis in the the two-parameter plane (−v sr , τ ). As shown in Figure 4d, the region surrounded by red dashed lines represents the SWD oscillation region (II) in Figure 3d, and the region surrounded by pink dashed lines is the SWD oscillation region as the DBS is applied. Comparing the two regions of Figures 3d, 4d, we find that the region (II) has an obvious reduction in Figure 4d. Additionally, the dominant frequency is shown in Figure 4e to estimate the frequency characteristics of different oscillation patterns. For clearer vision, in Figure 4f, we enlarge the SWD oscillation region of Figure 4e, which indeed falls into 2-4 Hz frequency range. The stimulation parameters we used keep within a reasonable range, which have been applied in previous studies (Mirski et al., 1997;Usui et al., 2005). Therefore, we can conclude that the DBS can eliminate absence seizures.

Effects of DBS on SWD Oscillation Pattern with Different Parameters
So far, we have confirmed that absence seizures can be simulated and eliminated in our simplified mean-field model. The choice of the stimulation parameters including frequency, amplitude and positive input duration is also critical to the theoretical and clinical researches. Therefore, we continue to investigate the influence of different stimulation parameters on the SWD oscillation region in the following studies.
Currently, we will explore the control effect of DBS current on absence seizures when we change the frequency f . As shown in Figure 5a, two-dimensional state analysis has been obtained and the white dashed line represents the largest parameter interval of the inhibitory coupling strength −v sr in the SWD oscillation region. For the section below the white dotted line, when the frequency f is decreased, the parameter interval of −v sr in the SWD oscillation region diminishes slowly, i.e., the suppression of SWDs can be achieved slowly. Whereas for the section above the white dotted line, the parameter range of −v sr in the SWD oscillation region reaches the minimum value quickly as the f is increased , i.e., the SWDs can be inhibited by raising the parameter f . Here the double arrow represents the parameter interval of −v sr in the SWD oscillation region. Moreover, the corresponding dominant frequency of Figure 5a is shown in Figure 5b. Although the distribution of frequency seems unclear, the SWD oscillation state indeed falls into 2-4 Hz frequency The enlargement (f) here is a local focus on (e), with the purpose of observing a obvious phenomenon. Different colors in (d) represent the different dynamical states: the saturation state (I), the SWD oscillation state (II), the simple oscillation region (III), the low firing state (IV). The region surrounded by pink dashed lines in (d,e) represents the SWDs falling into the 2-4 Hz frequency range after adding the DBS to TRN, however the regions surrounded by red dashed lines stand for the 2-4 Hz SWDs without DBS. The SWD oscillation pattern has a obvious reduction, i.e., the absence seizures can be eliminated through the DBS to TRN. range. On the whole, DBS strongly reduce the persistent absence seizures for low-frequency (<10 Hz) and high-frequency (>60 Hz) while intermediate frequency stimulation (around 45 Hz) have no effect. These results are nearly similar with previous researches (Rubin and Terman, 2004;Mina et al., 2013).
We further discuss the stimulation parameter δ. As presented in Figure 5c, the bifurcation diagram looks very regular, especially, the SWD oscillation region is divided into three parts, each one is similar to a rectangle. For a relative weaker δ, the parameter range of −v sr in the SWD oscillation region achieves the maximum value and keeps unchanged. With increasing δ, the interval of −v sr in the SWD oscillation region obviously decreases. Especially, when the positive input duration δ exceeds about 0.75 ms, the interval of SWDs almost disappears. It is clear that the SWDs can be eliminated through increasing the parameter δ. Figure 5d illustrates the corresponding dominant frequency analysis of Figure 5c and the corresponding SWD oscillation state falls into 2-4 Hz frequency range. At last, we investigate the amplitude a. The two-dimensional bifurcation diagram and frequency characteristic are demonstrated in the plane (−v sr , a). As shown in Figure 5e, as a is increased, the parameter interval of the SWD oscillation region about inhibitory coupling strength −v sr becomes smaller before it disappears. As shown in Figure 5f, the typical absence seizures falls into 2-4 Hz frequency range. Hence, it is shown that the suitable stimulations can contribute to the elimination of absence seizures through the DBS to TRN.

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

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