Bidirectional control of absence seizures by the basal ganglia: a computational evidence.

Absence epilepsy is believed to be associated with the abnormal interactions between the cerebral cortex and thalamus. Besides the direct coupling, anatomical evidence indicates that the cerebral cortex and thalamus also communicate indirectly through an important intermediate bridge-basal ganglia. It has been thus postulated that the basal ganglia might play key roles in the modulation of absence seizures, but the relevant biophysical mechanisms are still not completely established. Using a biophysically based model, we demonstrate here that the typical absence seizure activities can be controlled and modulated by the direct GABAergic projections from the substantia nigra pars reticulata (SNr) to either the thalamic reticular nucleus (TRN) or the specific relay nuclei (SRN) of thalamus, through different biophysical mechanisms. Under certain conditions, these two types of seizure control are observed to coexist in the same network. More importantly, due to the competition between the inhibitory SNr-TRN and SNr-SRN pathways, we find that both decreasing and increasing the activation of SNr neurons from the normal level may considerably suppress the generation of spike-and-slow wave discharges in the coexistence region. Overall, these results highlight the bidirectional functional roles of basal ganglia in controlling and modulating absence seizures, and might provide novel insights into the therapeutic treatments of this brain disorder.

Introduction

Absence epilepsy is a generalized non-convulsive seizure disorder of the brain, mainly occurring in the childhood years [1]. A typical attack of absence seizures is characterized by a brief loss of consciousness that starts and terminates abruptly, and meanwhile an electrophysiological hallmark, i.e. the bilaterally synchronous spike and wave discharges (SWDs) with a slow frequency at approximately 2–4 Hz, can be observed on the electroencephalogram (EEG) of patients [1], [2]. There is a broad consensus that the generation of SWDs during absence seizures is due to the abnormal interactions between cerebral cortex and thalamus, which together form the so-called corticothalamic system. The direct evidence in support of this view is based on simultaneous recordings of cortex and thalamus from both rodent animal models and clinical patients [3]–[5]. Recent computational modelling studies on this prominent brain disorder also approved the above viewpoint and provided more deep insights into the possible generation mechanism of SWDs in the corticothalamic system [6]–[13].

The basal ganglia comprise a group of interconnected subcortical nucleus and, as a whole, represent one fundamental processing unit of the brain. It has been reported that the basal ganglia are highly associated with a variety of brain functions and diseases, such as cognitive [14], emotional functions [15], motor control [16], Parkinson's disease [17], [18], and epilepsy [19], [20]. Anatomically, the basal ganglia receive multiple projections from both the cerebral cortex and thalamus, and in turn send both direct and indirect output projections to the thalamus. These connections enable the activities of the basal ganglia to influence the dynamics of the corticothalamic system. Therefore, it is naturally expected that the basal ganglia may provide an active role in mediating between seizure and non-seizure states for absence epileptic patients. Such hypothesis has been confirmed by both previous animal experiments [19], [21]–[23] and recent human neuroimage data [20], [24], [25]. Nevertheless, due to the complicated interactions between basal ganglia and thalamus, the underlying neural mechanisms on how the basal ganglia control the absence seizure activities are still remain unclear.

From the anatomical perspective, the substantia nigra pars reticulata (SNr) is one of the major output nucleus of the basal ganglia to thalamus. Previous experimental studies using various rodent animal models have demonstrated that suitable changes in the firing of SNr neurons can modulate the occurrence of absence seizures [21]–[23], [26]. Specifically, it has been found that pharmacological inactivation of the SNr by injecting -aminobutyric acids (GABA) agonists or glutamate antagonists suppresses absence seizures [21], [22]. Such antiepileptic effect was supposed to be attributed to the overall inhibitory effect of the indirect pathway from the SNr to thalamic reticular nucleus (TRN) relaying at superior colliculus [21], [22]. In addition to this indirect inhibitory pathway, it is known that the SNr also contains GABAergic neurons directly projecting to the TRN and specific relay nuclei (SRN) of thalamus [27], [28]. Theoretically, changing the activation level of SNr may also significantly impact the firing activities of SRN and TRN neurons [28], [29]. This contribution might further interrupt the occurrence of SWDs in the corticothalamic system, thus providing an alternative mechanism to regulate typical absence seizure activities. To our knowledge, however, so far the precise roles of these direct basal ganglia-thalamic pathways in controlling absence seizures are not completely established.

To address this question, we develop a realistic mean-field model for the basal ganglia-corticothalamic (BGCT) network in the present study. Using various dynamic analysis techniques, we show that the absence seizures are controlled and modulated either by the isolated SNr-TRN pathway or the isolated SNr-SRN pathway. Under suitable conditions, these two types of modulations are observed to coexist in the same network. Importantly, in this coexist region, both low and high activation levels of SNr neurons can suppress the occurrence of SWDs due to the competition between these two direct inhibitory basal ganglia-thalamic pathways. These findings clearly outline a bidirectional control of absence seizures by the basal ganglia, which is a novel phenomenon that has never been identified both in previous experimental and modelling studies. Our results, on the one hand, further improve the understanding of the significant role of basal ganglia in controlling absence seizure activities, and on the other hand, provide testable hypotheses for future experimental studies.

Materials and Methods

Model

We build a biophysically based model that describes the population dynamics of the BGCT network to investigate the possible roles of basal ganglia in the control of absence seizures. The network framework of this model is inspired by recent modelling studies on Parkinson's disease [30], [31], which is shown schematically in Fig. 1. The network totally includes nine neural populations, which are indicated as follows: excitatory pyramidal neurons (EPN); inhibitory interneurons (IIN); TRN; SRN; striatal D1 neurons; striatal D2 neurons; SNr; globus pallidus external (GPe) segment; subthalamic nucleus (STN). Similar to other modelling studies [30]–[32], we do not model the globus pallidus internal (GPi) segment independently but consider SNr and GPi as a single structure in the present study, because they are reported to have closely related inputs and outputs, as well as similarities in cytology and function. Three types of neural projections are contained in the BGCT network. For sake of clarity, we employ different line types and heads to distinguish them (see Fig. 1). The red lines with arrow heads denote the excitatory projections mediated by glutamate, whereas the blue solid and dashed lines with round heads represent the inhibitory projections mediated by and , respectively. It should be noted that in the present study the connections among different neural populations are mainly inspired by previous modelling studies [30], [31]. Additionally, we also add the connection sending from SNr to TRN in our model, because recent anatomical findings have provided evidence that the SNr also contains GABAergic neurons directly projecting to the TRN [27]–[29].

Figure 1

Framework of the basal ganglia-corticothalamic network.

Neural populations include excitatory pyramidal neurons (EPN); inhibitory interneurons (IIN); thalamic reticular nucleus (TRN); specific relay nuclei (SRN); striatal D1 neurons; striatal D2 neurons; substantia nigra pars reticulata (SNr); globus pallidus external (GPe) segment; subthalamic nucleus (STN). Note that we do not model the globus pallidus internal (GPi) segment independently but consider SNr and GPi as a single structure in this work. Red lines with arrow heads denote the excitatory projections mediated by glutamate receptors. Blue solid and dashed lines with round heads represent the inhibitory projections mediated by GABAA and GABAB receptors, respectively. Compared with the traditional model of corticothalamic system on absence seizures, the basal ganglia are also included in our biophysical model.

The dynamics of neural populations are characterized by the mean-field model [9], [11], [33]–[35], which was proposed to study the macroscopic dynamics of neural populations in a simple yet efficient way. The first component of the mean-field model describes the average response of populations of neurons to changes in cell body potential. For each neural population, the relationship between the mean firing rate and its corresponding mean membrane potential satisfies an increasing sigmoid function, given bywhere indicate different neural populations, denotes the maximum firing rate, r represents the spatial position, is the mean firing threshold, and is the threshold variability of firing rate. If exceeds the threshold , the neural population fires action potentials with an average firing rate . It should be noted that the sigmoid shape of is physiologically crucial for this model, ensuring that the average firing rate cannot exceed the maximum firing rate . The changes of the average membrane potential at the position r, under incoming postsynaptic potentials from other neurons, are modeled as [9], [33]–[36] where is a differential operator representing the dendritic filtering of incoming signals. and are the decay and rise times of cell-body response to incoming signals, respectively. is the coupling strength between neural populations of type and type . is the incoming pulse rate from the neural population of type to type . For simplicity, we do not consider the transmission delay among most neural populations in the present work. However, since the functions via second messenger processes, a delay parameter is introduced to its incoming pulse rate (i.e., ) to mimic its slow synaptic kinetics. This results in a delay differential equation in the final mathematical description of the BGCT model. Note that the similar modelling method has also been used in several previous studies [13], [37].

In our system, each neural population gives rise to a field of pulses, which travels to other neural population at a mean conduction velocity . In the continuum limit, this type of propagation can be well-approximated by a damped wave equation [9], [33]–[35], [38]:Here is the Laplacian operator (the second spatial derivative), is the characteristic range of axons of type , and governs the temporal damping rate of pulses. In our model, only the axons of cortical excitatory pyramidal neurons are assumed to be sufficiently long to yield significant propagation effect. For other neural populations, their axons are too short to support wave propagation on the relevant scales. This gives (). Moreover, as one of typical generalized seizures, the dynamical activities of absence seizures are believed to occur simultaneously throughout the brain. A reasonable simplification is therefore to assume that the spatial activities are uniform in our model, which has been shown as the least stable mode in models of this class [33], [34], [36]. To this end, we ignore the spatial derivative and set in Eq. (4). Accordingly, the propagation effect of cortical excitatory axonal field is finally given by [33], [34], [36]:where . For the population of cortical inhibitory interneurons, the BGCT model can be further reduced by using and , which is based on the assumption that intracortical connectivities are proportional to the numbers of synapses involved [9], [13], [33]–[36]. It has been demonstrated that by making these above reductions, the developed BGCT model becomes computationally more tractable without significant deteriorating the precision of numerical results.

We then rewrite above equations in the first-order form for all neural populations. Following above assumptions, we use Eqs. (1)–(3) and (5) for modelling the dynamics of excitatory pyramidal neurons, and Eqs. (1)–(3) for modelling the dynamics of other neural populations. This yields the final mathematical description of the BGCT model given as follows: where In Eq. (10), the superscript T denotes transposition. The detailed expression of for different neural populations is represented by , and , given bywith Here the variable in Eq. (15) denotes the delay and the parameter in Eq. (16) represents the constant nonspecific subthalamic input onto SRN.

The parameters used in our BGCT model are compatible with physiological experiments and their values are adapted from previous studies [9], [11], [13], [30], [31], [36]. Unless otherwise noted, we use the default parameter values listed in Table 1 for numerical simulations. Most of the default values of these parameters given in Table 1 are based on either their nominal values or parameter ranges reported in above literature. A small number of parameters associated with the basal ganglia (i.e., , and ) are adjusted slightly, but still within their normal physiological ranges, to ensure our developed model can generate the stable 2–4 Hz SWDs under certain conditions. Note that due to lack of quantitative data, the coupling strength of the SNr-TRN pathway needs to be estimated. Considering that the SNr sends GABAergic projections both to SRN and TRN and also both of these two nuclei are involved in thalamus, it is reasonable to infer that the coupling strengths of these two pathways are comparable. For simplicity, here we chose by default. In the following studies, we also change (decrease or increase) the value of several folds by employing a scale factor (see below) to examine how the inhibition from the SNr-TRN pathway regulates absence seizures. Additionally, during this study, several other critical parameters (i.e., , and ) are also varied within certain ranges to obtain different dynamical states and investigate their possible effects on the modulation of absence seizures.

Table 1

Default parameter values.

Symbol	Description	Value	References
,	Cortical maximum firing rate		[9], [11], [13]
,	Striatum maximum firing rate		[30], [31]
	SNr maximum firing rate		[30], [31]
	GPe maximum firing rate		[30], [31]
	STN maximum firing rate		[30], [31]
	SRN maximum firing rate		[9], [11], [13]
	TRN maximum firing rate		[9], [11], [13]
,	Mean firing threshold of cortical populations		[9], [11], [13], [36]
,	Mean firing threshold of striatum		[30], [31]
	Mean firing threshold of SNr		[30], [31]
	Mean firing threshold of GPe		[30], [31]
	Mean firing threshold of STN		[30], [31]
	Mean firing threshold of SRN		[9], [11], [13], [36]
	Mean firing threshold of TRN		[9], [11], [13], [36]
	Self-coupling strength of EPN		[11], [13]
	Coupling strength from IIN to EPN		[11], [13]
	Coupling strength from EPN to TRN		[9], [13]
	Coupling strength from SRN to TRN		[9], [13]
	Coupling strength from TRN to SRN		[9], [36]
	Coupling strength from EPN to striatal D1 neurons		[30], [31]
	Self-coupling strength of striatal D1 neurons		[30], [31]
	Coupling strength from SRN to striatal D1 neurons		[30], [31]
	Coupling strength from EPN to striatal D2 neurons		[30], [31]
	Self-coupling strength of striatal D2 neurons		[30], [31]
	Coupling strength from SRN to striatal D2 neurons		[30], [31]
	Coupling strength from striatal D1 neurons to SNr		[30], [31]
	Coupling strength from GPe to SNr		[30], [31]
	Coupling strength from STN to SNr		[30], [31]
	Coupling strength from striatal D2 neurons to GPe		[30], [31]
	Self-coupling strength of GPe		[30], [31]
	Coupling strength from STN to GPe		[30], [31]
	Coupling strength from GPe to STN		[30], [31]
	Coupling strength from SRN to EPN		[9], [13]
	Coupling strength from EPN to SRN		[9]
	Coupling strength from EPN to STN		[30], [31]
	Coupling strength from SNr to SRN		[30], [31]
	Coupling strength from SNr to TRN		Estimated
	Cortical damping rate		[9], [11], [13]
	Time delay due to slow synaptic kinetics of GABAB		[13]
	Synaptodendritic decay time constant		[9], [11], [13], [36]
	Synaptodendritic rise time constant		[9], [11], [13], [36]
	Threshold variability of firing rate		[9], [11], [36]
	Nonspecific subthalamic input onto SRN		[9], [11], [36]

Model parameters employed in the present study. Unless otherwise noted, we use these default parameter values for simulations.

Data analysis

In the present study, several data analysis methods are employed to quantitatively evaluate the dynamical states as well as the properties of SWDs generated by the model. To reveal critical transitions between different dynamical states, we perform the bifurcation analysis for several key parameters of the model. For one specific parameter, the bifurcation diagram is simply obtained by plotting the “stable” local minimum and maximum values of cortical excitatory axonal fields (i.e., ) over changes in this parameter [11], [39]. To this end, all simulations are executed for sufficiently long time (10 seconds of simulation time, after the system reaches its stable time series), and only the local minimum and maximum values obtained from the latter stable time series are used. Using the above bifurcation analysis, we can also easily distinguish different dynamical states for combined parameters. Such analysis technique allows us to further identify different dynamical state regions in the two-parameter space (for example, see Fig. 2D). On the other hand, the power spectral analysis is used to estimate the dominant frequency of neural oscillations. To do this, the power spectral density is obtained from the time series (over a period of 10 seconds) by using the fast Fourier transform. Then, the maximum peak frequency is defined as the dominant frequency of neural oscillations. It should be noted that, by combining the results of both the state and frequency analysis, we can outline the SWD oscillation region that falls into the 2–4 Hz frequency range in the two-parameter space (for example, see the asterisk region in Fig. 2E). Moreover, we calculate the mean firing rates (MFRs) for several key neural populations in some figures. To compute the MFRs, all corresponding simulations are performed up to 25 seconds and the data from 5 to 25 seconds are used for statistical analysis. To obtain convincing results, we carry out 20 independent simulations with different random seeds for each experimental setting, and report the averaged result as the final result. Finally, in some cases, we also compute the low and high triggering mean firing rates (TMFRs) for SNr neurons. In the following simulations, we find that the mean firing rate of SNr neurons is increased with the growth of the excitatory coupling strength , which serves as a control parameter to modulate the activation level of SNr in our work (see the Results section). Based on this property, the low and high TMFRs can be determined by the mean firing rates of SNr neurons occurring at the boundaries of the typical region of 2–4 Hz SWDs (for example, see the black dashed lines in Fig. 3B).

Figure 2

Absence seizure activities induced by the slow kinetics of

receptors in TRN. A, B: Bifurcation diagrams of as a function of the TRN-SRN inhibitory coupling strength (A) and the delay parameter (B), respectively. Four different dynamical states can be observed from the time series of , which are: the saturation state (I), the SWD oscillation state (II), the simple oscillation state (III) and the low firing state (IV). C: Typical time series of correspond to the above four dynamical states. Here we set and chose (I), (II), (III), (IV), respectively. The colors in bifurcation diagrams (A) and (B) correspond to the typical time series plotted in (C). D, E: The state analysis (D) and frequency analysis (E) in the () panel. Different colors in (D) represent different dynamical state regions, corresponding to those dynamical states given in (A), (B) and (C). The yellow asterisk regions surrounded by black dashed lines in (D) and (E) represent the SWD oscillation regions falling into the 2–4 Hz frequency range. The other symbols in (D) and (E) are linked to parameter values used for different typical time series in (C): I (white filled circle), II (red filled square), III (yellow filled triangle), and IV (purple filled star). For all simulations, we set .

Figure 3

Control of absence seizures by the isolated SNr-TRN pathway.

A, B: The state analysis (A) and frequency analysis (B) in the () panel. Here is the inhibitory coupling strength of the TRN-SRN pathway, whereas is the excitatory coupling strength of the STN-SNr pathway. Different colors in (A) represent different dynamical state regions: the saturation region (I), the SWD oscillation region (II), the simple oscillation region (III) and the low firing region (IV). The suppression of SWDs appears to the right of the white dashed line in (A), where the white down arrow indicates that the SWD oscillation can be inhibited by decreasing . The yellow asterisk region surrounded by black dashed lines in (B) denotes the typical 2–4 Hz SWD oscillation region. C, D: Bifurcation diagrams of as a function of for different . The strengths of the inhibitory projections from the TRN to SRN are set as (C) and (D), respectively. Different colors in (C) and (D) represent different dynamical state regions, corresponding to those in phase diagram (A). E: The mean firing rates (MFRs) of several key neural populations as a function of , with . Here four neural populations are considered: SNr (blue dot), excitatory pyramidal neurons (green asterisk), SRN (black circle) and TRN (red square). Two black dashed lines in (E) represent the occurring positions of the low and high triggering mean firing rates (TMFRs), respectively. F: The low (red filled circle) and high (green filled square) TMFRs as a function of . For all simulations, the SNr-SRN pathway is artificially blocked (i.e., ).

Numerical simulation

All network simulations are written and performed under the MATLAB environment. The aforementioned dynamical equations are integrated by using the standard fourth-order Runge-Kutta method, with a fixed temporal resolution of [40]. In additional simulations, it turns out that the chosen integration step is sufficiently small to ensure the numerical accuracy of our developed BGCT model. The computer codes used in the present study are available on ModelDB (https://senselab.med.yale.edu/ModelDB/showmodel.asp?model=152113). The fundamental implementation of the BGCT model is provided as supplementary information to this paper (we also provide a XPPAUT code for comparison [41]; see Text S1 and S2).

Results

Slow kinetics of receptors in TRN create absence seizure activities for the BGCT model

Previous studies have suggested that the slow kinetics of receptors in TRN are a candidate pathological factor contributing to the generation of absence seizures both in animal experiments and biophysical models of corticothalamic network [7], [13], [42], [43]. To explore whether this mechanism also applies to the developed BGCT model, we perform one-dimensional bifurcation analysis for the inhibitory coupling strength and the delay parameter , respectively. The corresponding bifurcation diagrams and typical time series of are depicted in Figs. 2A–2C, which reveal that different dynamical sates emerge in our system for different values of and . When the coupling strength is too weak, the inhibition from TRN cannot effectively suppress the firing of SRN. In this case, due to the strong excitation from pyramidal neurons, the firing of SRN rapidly reaches a high level after the beginning of the simulation. Such high activation level of SRN in turn drives the firing of cortical neurons to their saturation states within one or two oscillation periods (region I). As the coupling strength grows, the inhibition from TRN starts to affect the firing of SRN. For sufficiently long , this causes our model to successively undergo two different oscillation patterns. The first one is the SWD oscillation pattern, in which multiple pairs of maximum and minimum values are found within each periodic complex (region II). Note that this oscillation pattern has been extensively observed on the EEG recordings of real patients during absence seizures [1]. The other one is the simple oscillation pattern, in which only one pair of maximum and minimum values appears within each periodic complex (region III). However, if the coupling strength is too strong, the firing of SRN is almost completely inhibited by TRN. In this situation, the model is kicked into the low firing region and no oscillation behavior can be observed anymore (region IV). Additionally, we also find that the model dynamics are significantly influenced by the delay, and only sufficiently long can ensure the generation of SWDs in the developed model (see Fig. 2B).

To check whether our results can be generalized within a certain range of parameters, we further carry out the two-dimensional state analysis in the () panel. As shown in Fig. 2D, the whole () panel is divided into four state regions, corresponding to those regions identified above. Unsurprisingly, we find that the BGCT model can generate the SWD oscillation pattern only for appropriately intermediate and sufficiently long . This observation is in consistent with our above finding, demonstrating the generalizability of our above results. To estimate the frequency characteristics of different oscillation patterns, we compute the dominant frequency based on the spectral analysis in the () panel. For both the simple and SWD oscillation patterns, the dominant frequency is influenced by and , and increasing their values can both reduce the dominant frequency of neural oscillations (Fig. 2E). However, compared to , our results indicate that the delay may have a more significant effect on the dominant oscillation frequency (Fig. 2E). By combining the results in Figs. 2D and 2E, we roughly outline the SWD oscillation region that falls into the 2–4 Hz frequency range (asterisk region). It is found that most of, but not all, the SWD oscillation region is contained in this specific region. Here we emphasize the importance of this specific region, because the SWDs within this typical frequency range is commonly observed during the paroxysm of absence epilepsy in human patients [1], [2].

Why can the slow kinetics of receptors in TRN induce absence seizure activities? Anatomically, the SRN neurons receive the TRN signals from the inhibitory pathway mediated by both and receptors. Under suitable condition, the double suppression caused by these two types of GABA receptors occurring at different time instants may provide an effective mechanism to create multiple firing peaks for the SRN neurons (see below). Such firing pattern of SRN in turn impacts the dynamics of cortical neurons, thus leading to the generation of SWDs. It should be noted that, during the above processes, both and play critical roles. In each oscillation period, after the -induced inhibition starts to suppress the firing of SRN neurons, these neurons need a certain recovery time to restore their mean firing rate to the rising state. Theoretically, if this recovery time is shorter than the delay, another firing peak can be introduced to SRN neurons due to the latter -induced inhibition. The above analysis implies that our model requires a sufficient long delay to ensure the occurrence of SWDs. However, as described above, too long is also a potential factor which may push the dominant frequency of SWDs beyond the typical frequency range. For a stronger , the inhibition caused by is also strong. In this situation, it is obvious that the SRN neurons need a longer time to restore their firing rate. As a consequent, a relatively longer is required for the BGCT model to ensure the occurrence of SWDs for stronger (see Fig. 2D).

These findings provide consistent evidence that our developed BGCT model can replicate the typical absence seizure activities utilizing previously verified pathological mechanism. Because we do not change the normal parameter values for basal ganglia during above studies, our results may also indicate that, even though the basal ganglia operate in the normal state, the abnormal alteration within the corticothalamic system may also trigger the onset of absence epilepsy. Throughout the following studies, we set for all simulations. For this choice, the delay parameter is within the physiological range and modest, allowing the generation of SWD oscillation pattern while preserving its dominant frequency around 3 Hz in most considered parameter regions. It should be noted that, in additional simulations, we have shown that by slightly tuning the values of several parameters our developed BGCT model is also powerful to reproduce many other typical patterns of time series, such as the alpha and beta rhythms (see Fig. S1), which to a certain extent can be comparable with real physiological EEG signals [9], [36].

Control of absence seizures by the isolated SNr-TRN pathway

Using the developed BGCT model, we now investigate the possible roles of basal ganglia in controlling absence seizure activities. Here we mainly concentrate on how the activation level of SNr influence the dynamics generated by the model. This is because, on the one hand, the SNr is one of chief output nucleus of the basal ganglia to thalamus, and on the other hand, its firing activity has been found to be highly associated with the regulation of absence seizures [21], [22]. To this end, the excitatory coupling strength is employed to control the activation level of SNr and a three-step strategy is pursued in the present work. In this and next subsections, we assess the individual roles of two different pathways emitted from SNr to thalamus (i.e., the SNr-TRN and SNr-SRN pathways) in the control of absence seizures and discuss their corresponding biophysical mechanisms, respectively. In the final two subsections, we further analyze the combination effects of these two pathways on absence seizure control and extend our results to more general cases.

To explore the individual role of the SNr-TRN pathway, we estimate both the state regions and frequency characteristics in the () panel. Note that during these investigations the SNr-SRN pathway is artificially blocked (i.e., ). With this “naive” method, the modulation of absence seizure activities by the SNr-SRN pathway is removed and the effect caused by the SNr-TRN pathway is theoretically amplified to the extreme. Similar to previous results, we find that the whole () panel can be also divided into four different regions (Fig. 3A). These regions are the same as those defined above. For weak inhibitory coupling strength , increasing the excitatory coupling strength moves the model dynamics from the SWD oscillation state to the saturation state. Here we have to notice that the saturation state is a non-physiological brain state even though it does not belong to typical seizure activities. In strong region, the suppression of SWDs is observed by decreasing the excitatory coupling strength , suggesting that inactivation of SNr neurons may result in seizure termination through the SNr-TRN pathway (Fig. 3A, right side). For strong enough inhibitory coupling strength , such suppression effect is rather remarkable that sufficiently low activation of SNr can even kick the network dynamics into the low firing region (compare the results in Figs. 3C and 3D).

The SNr-TRN pathway induced SWD suppression is complicated and its biophysical mechanism is presumably due to competition-induced collision. On the one side, the decrease of excitatory coupling strength inactivates the SNr (Fig. 3E, top panel), which should potentially enhance the firing of TRN neurons. On the other side, however, increasing the activation level of TRN tends to suppress the firing of SRN, which significantly reduces the firing of cortical neurons and in turn inactivates the TRN neurons. Furthermore, the inactivation of cortical neurons also tends to reduce the firing level of TRN neurons. As the excitatory coupling strength is decreased, the collision caused by such complicated competition and information interactions finally leads to the inactivation for all the TRN, SRN, and cortical neurons (Fig. 3E, bottom panel), which potentially provides an effective mechanism to destabilize the original pathological balance within the corticothalamic system, thus causing the suppression of SWDs.

Indeed, we find that not only the dynamical state but also the oscillation frequency is greatly impacted by the activation level of SNr, through the SNr-TRN pathway. For both the simple and SWD oscillation patterns, increasing the excitatory strength can enhance their dominant frequencies. The combined results of Figs. 3A and 3B reveal that, for a fixed , whether the model can generate the SWDs within the typical 2–4 Hz is determined by at least one and often two critical values of (Fig. 3B, asterisk region). Because the activation level of SNr is increased with the growth of , this finding further indicates that, due to effect of the SNr-TRN pathway, the model might exist the corresponding low and high triggering mean firing rates (TMFRs) for SNr neurons (Fig. 3E, dashed lines). If the long-term mean firing rate of SNr neurons falls into the region between these two TMFRs, the model can highly generate typical 2–4 Hz SWDs as those observed on the EEG recordings of absence epileptic patients. In Fig. 3F, we plot both the low and high TMFRs as a function of the inhibitory coupling strength . With the increasing of , the high TMFR grows rapidly at first and then reaches a plateau region, whereas the low TMFR almost linearly increases during this process. Consequently, it can be seen that these two critical TMFRs approach each other as the inhibitory coupling strength is increased until they almost reach an identical value (Fig. 3F).

The above findings indicate that the SNr-TRN pathway may play a vital role in controlling the absence seizures and appropriately reducing the activation level of SNr neurons can suppress the typical 2–4 Hz SWDs. The similar antiepileptic effect induced by inactivating the SNr has been widely reported in previous electrophysiological experiments based on both genetic absence epilepsy rats and tottering mice [21]–[23], [26]. Note that, however, in literature such antiepileptic effect by reducing the activation of SNr is presumed to be accomplished through the indirect SNr-TRN pathway relaying at superior colliculus [21], [22]. Our computational results firstly suggest that such antiepileptic process can be also triggered by the direct SNr-TRN GABAergic projections. Combining these results, we postulate that for real absence epileptic patients both of these two pathways might work synergistically and together provide a stable mechanism to terminate the onset of absence epilepsy.

Control of absence seizures by the isolated SNr-SRN pathway

We next turn on the SNr-SRN pathway and investigate whether this pathway is also effective in the control of absence seizures. Similar to the previous method, we artificially block the SNr-TRN pathway (i.e., ) to enlarge the effect of the SNr-SRN pathway to the extreme. Fig. 4A shows the two-dimensional state analysis in the () panel, and again the whole panel is divided into four different state regions. Compared to the results in Fig. 3A, the suppression of SWDs appears in a relatively weaker region by increasing the excitatory coupling strength . This finding suggests that the increase in the activation of SNr can also terminate the SWDs, but through the SNr-SRN pathway. For relatively weak within the suppression region, the SNr-SRN pathway induced suppression is somewhat strong. In this case, the high activation level of SNr directly kicks the network dynamics into the low firing region, without undergoing the simple oscillation state (Fig. 4C and compare with Fig. 4C). Note that this type of state transition is a novel one which has not been observed in the SWD suppression caused by the SNr-TRN pathway. For relatively strong within the suppression region, the double peak characteristic of SWDs generated by our model is weak. In this situation, as the inhibitory coupling strength is increased, we observe that the network dynamics firstly transit from the SWD oscillation state to the simple oscillation state, and then to the low firing state (Fig. 4C).

Figure 4

Control of absence seizures by the isolated SNr-SRN pathway.

A, B: The state analysis (A) and frequency analysis (B) in the () panel. Here is the inhibitory coupling strength of the TRN-SRN pathway, whereas is the excitatory coupling strength of the STN-SNr pathway. Similar to previous results, four different dynamical state regions are observed: the saturation region (I), the SWD oscillation region (II), the simple oscillation region (III) and the low firing region (IV), which correspond to those defined in Fig. 2 (D). The region between two white dashed lines in (A) represents the suppression region of SWDs, where the white up arrow indicates that the SWD oscillation can be inhibited by increasing . The yellow asterisk region surrounded by black dashed lines in (B) denotes the typical 2–4 Hz SWD oscillation region. C: Two different types of SWD suppressions caused by linearly increasing . Top (C1): The value of as a function of time. Middle (C2): Corresponding trace for . Bottom (C3): Corresponding trace for . D: The MFRs of several key neural populations as a function of , with . Here four neural populations are considered: SNr (blue dot), excitatory pyramidal neurons (green asterisk), SRN (black circle) and TRN (red square). The black dashed line in (D) represents the occurring position of TMFR. E: The TMFR as a function of . Note that the SNr-TRN pathway is artificially blocked (i.e., ) for all simulations.

To understand how the SNr-SRN pathway induced SWD suppression arises, we present the mean firing rates of several key neural populations within the corticothalamic system, as shown in Fig. 4D. It can be seen that increasing the strength significantly improves the activation level of SNr (Fig. 4D, top panel), which in turn reduces the firing of SRN neurons (Fig. 4D, bottom panel). The inactivation of SRN neurons further suppresses the mean firing rates for both cortical and TRN neurons (Fig. 4D, bottom panel). These chain reactions lead to the overall inhibition of firing activities in the corticothalamic system, which weakens the double peak shaping effect due to the slow kinetics of receptors in TRN. For strong , such weakening effect is considerable, thus causing the suppression of SWDs. Our results provide the computational evidence that high activation of SNr can also effectively terminate absence seizure activities by the strong inhibition effect from the SNr-SRN pathway. Compared to the SWD suppression induced by the SNr-TRN pathway, it is obvious that the corresponding biophysical mechanism caused by the SNr-SRN pathway is simpler and more direct.

Moreover, our two-dimensional frequency analysis indicates that the dominant frequency of neural oscillations depends on the excitatory coupling strength (see Fig. 4B). For a constant , progressive increase of reduces the dominant frequency, but not in a very significant fashion. Thus, we find that almost all the SWD oscillation region identified in Fig. 4A falls into the typical 2–4 Hz frequency range (Fig. 4B, asterisk region). Unlike the corresponding results presented in previous subsection, the combination results of Figs. 4A and 4B demonstrate that the BGCT model modulated by the isolated SNr-SRN pathway only exhibits one TMFR for SNr neurons. For a suitably fixed strength , the generation of SWDs can be highly triggered when the mean firing rate of SNr neurons is lower than this critical firing rate (Fig. 4D, dashed line). With the increasing of , we observe that this TMFR rapidly reduces from a high value to a low value (Fig. 4E). Note that this decreasing tendency is in contrast with our previous finding based on the model modulated by the isolated SNr-TRN pathway (compared to the results in Fig. 3F).

Taken together, these observations suggest that increasing the activation level of SNr neurons also significantly influences the dynamics of the corticothalamic system and causes the suppression of absence seizure activities. To the best of our knowledge, this is a new finding that underscores the importance of the direct inhibitory SNr-SRN pathway in controlling and modulating absence seizure activities. It is reasonable to believe that several other external factors, which are able to enhance the activation level of SNr, may also result in the termination of absence seizures due to the similar mechanism.

Competition-induced bidirectional control of absence seizures by the basal ganglia

So far, we have confirmed that the absence seizure activities generated by the BGCT model can be inhibited by either the isolated SNr-TRN pathway or the isolated SNr-SRN pathway, through different biophysical mechanisms. In real brain, however, both of these pathways should be available and work together at the same time. Thus, an important and naturally arising question is whether these two types of seizure control can coexist in the same network, and if possible, whether this feature can be maintained in certain range of parameters. To address this issue, we introduce a scale factor and set the coupling strength . By utilizing this method, we can flexibly control the relative coupling strength between the SNr-TRN and SNr-SRN pathways and discuss the combination roles of these two inhibitory pathways in detail.

Fig. 5A shows the state analysis in the () panel with . Unlike the previous results, here we only discover three dynamical state regions, which correspond to: the SWD oscillation state (II), the simple oscillation state (III), and the low firing state (IV). The disappearance of the saturation state is at least due to the following two reasons: (1) the double suppression from the SNr-TRN and SNr-SRN pathways and (2) the relatively strong inhibitory effect from TRN to SRN. As shown in Fig. 5A, the phenomenon of the SWD suppression appears in the strong region. For relatively weak within this suppression region, both increasing and decreasing the activation of SNr neurons from the normal level effectively inhibit the generation of SWDs (see Fig. 5A). Such bidirectional suppression behavior can be attributed to the effective competition between the SNr-TRN and SNr-SRN pathways. As the scale factor is increased, the enhancement of the SNr-TRN pathway breaks the original competition balance between these two inhibitory pathways. During this process, the inhibition from the SNr-TRN pathway progressively dominates the model dynamics. Accordingly, for sufficiently strong , the suppression of SWDs is only found by lowing the activation of SNr (Fig. 5A), which is consistent with our previous critical observation given in Fig. 3A.

Figure 5

Control of absence seizures by the combination effects of the SNr-TRN and SNr-SRN pathways.

A, B: The state analysis (A) and frequency analysis (B) in the () panel. Here is the scale factor, and is the excitatory coupling strength of the STN-SNr pathway. Unlike previous results, only three dynamical state regions are observed in the phase diagram (A): the SWD oscillation region (II), the simple oscillation region (III) and the low firing region (IV), corresponding to the dynamical states defined in Fig. 2(D). For relatively weaker scale factor , both increase and decrease in the activation level of SNr can inhibit the SWDs (double arrow, bidirectional suppression); whereas for sufficiently strong , only reducing the activation level of SNr suppresses the SWDs (single arrow, unidirectional suppression). In (A), the white dashed line represents the boundary of suppression region, and the red dashed line stands for the demarcation between the bidirectional and unidirectional suppression regions. The yellow asterisk region surrounded by black dashed lines in (B) denotes the SWD oscillation region that falls into the 2–4 Hz frequency range. C: The MFRs of several key neural populations as a function of , with the scale factor . Here four neural populations are considered: SNr (blue dot), excitatory pyramidal neurons (green asterisk), SRN (black circle) and TRN (red square). Two black dashed lines in (C) represent the occurring positions of low and high triggering mean firing rates (TMFRs), respectively. D: The low (red filled circle) and high (green filled square) TMFRs as a function of . For all simulations, the coupling strength of the TRN-SRN pathway is set as .

However, as described above, the increase in the excitatory coupling strength significantly reduces the dominant frequency of SWDs, through the SNr-TRN pathway. Therefore, although enhancing the activation level of SNr neurons cannot inhibit the SWDs directly, it tends to push the dominant frequency of SWDs below 2 Hz (Fig. 5B). By combining the results of Figs. 5A and 5B, we successfully outline the SWD oscillation region that falls into the 2–4 Hz frequency range (asterisk region). In this typical SWD region, the model exhibits both the low and high TMFRs of the SNr neurons for a fixed (Figs. 5B and 5C). The generation of SWDs within the typical frequency range can be highly triggered if the mean firing rate of SNr neurons is between these two critical TMFRs (Fig. 5C). Nevertheless, due to the competition between the SNr-TRN and SNr-SRN pathways, the firing activities of the TRN, SRN, and cortical neurons become more complicated, compared to the cases induced by any isolated pathway. For a constant , such competition creates a bell-shaped MFR curve for each key neural population by tuning the excitatory coupling strength , (see Fig. 5C). To further investigate how the scale factor impacts the low and high TMFRs, we plot these two TMFRs as a function of in Fig. 5D. With the increasing of relative strength , both of these two critical TMFRs are rapidly changed but in different fashions. The high TMFR is increased from a relatively low value to saturation, whereas the low TMFR is reduced from a relatively high value to 0 (see Fig. 5D). Obviously, such opposite tendencies of these two TMFRs are attributed to the combination effects of the SNr-TRN and SNr-SRN pathways.

Furthermore, we also find that both the suppression of SWDs and the typical 2–4 Hz SWD region are shaped by the strength of inhibitory projections from the TRN to SRN. In Figs. 6A and 6B, we perform a series of two-dimensional state and frequency analysis in the () panel for different values of . When the inhibitory coupling strength is too weak, the SWDs generated by the corticothalamic system is mainly controlled by the SNr-SRN pathway. In this case, the suppression of SWDs is observed in the intermediate region and only increasing the activation level of SNr can effectively inhibit the SWDs (Fig. 6A). As the coupling strength is increased, the inhibition from SNr-TRN pathway starts to influence the model dynamics. This introduces the competition between the SNr-SRN and SNr-TRN pathways, leading to the emergence of the bidirectional suppression of SWDs for intermediate (Fig. 6A). It is obvious that the higher the coupling strength , the stronger the inhibition effect caused by the SNr-TRN pathway. With increasing the value of , such strengthened inhibition significantly moves the boundary of the low TMFR toward higher values of , and thus notably shrinks the region of SWDs within the typical frequency range of 2–4 Hz (Fig. 6B).

Figure 6

Effect of inhibitory coupling strength

on the control of absence seizures by the SNr-TRN and SNr-SRN pathways. Tow-dimensional state analysis (A) and corresponding frequency analysis (B) in the () panel for different values of . Here only three dynamical state regions are observed in (A): the SWD oscillation region (II), the simple oscillation region (III) and the low firing region (IV). In (A1)–(A4), the regions marked by red diamonds denote the whole suppression regions of SWDs, the white dashed lines represent the boundaries of suppression regions, and the red dashed lines stand for the demarcations between the bidirectional (double arrow) and unidirectional (single arrow) suppression regions. In (B1)–(B4), the yellow asterisk regions surrounded by black dashed lines denote the typical 2–4 Hz SWD oscillation regions. From left to right, the strengths of inhibitory projections from the TRN to SRN are: (A1, B1), (A2, B2), (A3, B3), and (A4, B4), respectively. For better showing and comparing the bidirectional suppression regions among different subfigures, here we mainly consider the scale factor for all values of within the same and small interval from 0 to 1. As an additional comparison, both the state analysis and frequency analysis for in a relatively larger interval is given in Fig. S2.

These above observations emphasize the importance of the combination role of both the SNr-TRN and SNr-SRN pathways on the control of absence seizure activities. Quite remarkably, we observe that the bidirectional suppression of SWDs emerges under suitable conditions, that is, both increasing and reducing the activation levels of SNr neurons effectively suppress the SWDs. Such bidirectional suppression is determined and modulated by both the relative strength of these two inhibitory pathways and the strength of inhibitory projections from the TRN to SRN. This novel finding indicates the possible bidirectional control of absence seizures by the basal ganglia, which is induced by the competition between the SNr-TRN and SNr-SRN pathways. In additional simulations, we have found that several other factors, which can effectively change the activation level of SNr, can also lead to the bidirectional suppression of SWDs due to the similar mechanism, further demonstrating the generality of our results.

Bidirectional control of absence seizures by the basal ganglia might be extended to other pathological factors

In above subsections, we focused on one specific pathological factor and have computationally shown that the basal ganglia may bidirectionally control and modulate the typical absence seizure activities (i.e., the SWDs) induced by the slow synaptic kinetics of in TRN. A natural question to ask is whether such bidirectional control feature caused by the basal ganglia is a generalized regulatory mechanism for absence seizures. We argue that this might be true. To check this postulation, at least one additional SWD generation mechanism should be introduced into our model, and we need to examine whether the bidirectional control of absence seizures by the basal ganglia is also available for the new pathological factor. In literature, there are several other theories that are associated with the generation mechanisms of SWDs. A boldly accepted one is related to the transmission delay between the cerebral cortex and thalamus and, specifically, it has been found that suitably choosing such transmission delay can drive the corticothalamic system produce the SWDs [11], [33]. To apply this pathological factor in our model, here we block the pathway from TRN to SRN, and consider a bidirectional transmission delay between the cerebral cortex and thalamus, as that used in previous studies [11], [33]. Additionally, several coupling strengths within the corticothalamic loop are also needed to be adapted, because such transmission delay induced SWDs require strong interactions between the cerebral cortex and thalamus [11], [33]. The new added and modified model parameters that we used in this subsection are as follows [11], [33]: , , , , and . For simplicity, we term the current model as the modified model in the following studies.

Figs. 7A and 7B show an example pair of state analysis and frequency analysis in the () panel, respectively. As expected, due to the competition between the SNr-TRN and SNr-SRN pathways, we observe the significant bidirectional control feature for intermediate scale factor (Fig. 7A, the region between dashed lines). In this bidirectional region, both enhancing and lowing the excitatory coupling strength push the model dynamics from the SWD oscillation state into the simple oscillation state (Fig. 7C), thus inhibiting the generation of SWDs. This finding supports our above hypothesis that under suitable conditions the basal ganglia may control and modulate the absence seizure activities bidirectionally. However, compared to the results in Fig.5, we find that the bidirectional region appears in a relatively larger region for the modified model, and increasing the coupling strength cannot kick the model dynamics into the low firing state as well. This is not so surprising because these two SWD generation mechanisms that we used are similar but not completely identical. As we introduced above, the SWDs induced by the corticothalamic loop transmission delay require relatively stronger interactions between the cerebral cortex and thalamus, which essentially weaken the inhibitory effect from the SNr neurons. This might lead that the firing of SRN neurons cannot be fully suppressed even when the activation level of SNr reaches its saturation state. Another important finding that we discover here is that the dynamics of the modified model become complicated for large and strong (see Fig. 7A, upper right). In additional simulations, we further perform a series of state analysis for the modified model using the same group of parameter values but different random initial conditions (see Fig. S3). The corresponding results indicate that the modified model shows bistability (the simple oscillation state or the saturation state) in the large and strong region, and final model dynamics significantly depend on the initial conditions.

Figure 7

An example illustrating that the bidirectional control of absence seizures by the basal ganglia is also available for the corticothalamic loop transmission delay induced SWDs.

A, B: Tow-dimensional state analysis (A) and corresponding frequency analysis (B) in the () panel. Here is the scale factor, and represents the excitatory coupling strength of the STN-SNr pathway. Four different dynamical states are observed in (A): the saturation state (I), the SWD oscillation state (II), the simple oscillation state (III) and the low firing state (IV), which also correspond to previous figures. In (A), the region between two red dashed lines denotes the main bidirectional suppression region of SWDs, where the double arrow represents both increasing and decreasing the excitatory coupling strength can inhibit the SWDs. In (B), the yellow asterisk region surrounded by black dashed lines denotes the typical 2–4 Hz SWD oscillation region. (C) Three typical time series of for different values of , with . Here we choose (top), (middle) and (bottom), respectively. In all simulations, we set , , , , and .

In conclusion, these findings further stress the combination role of the inhibitory SNr-TRN and SNr-SRN pathways on the control of absence seizure activities. By combining all of our results, we postulate that the bidirectional control by the basal ganglia is possible a generalized regulatory mechanism for absence seizures and may be extendable to other pathological factors, even though the detailed bidirectional control behaviors may not be completely identical for different pathological factors.

Discussion

Using a mean-field macroscopic model that incorporates the basal ganglia, cerebral cortex and thalamus, we presented here the first investigation on how the basal ganglia control the absence epilepsy through the projections directly emitted from the SNr to several key nuclei of thalamus. Through simulations, we demonstrated that the absence seizure activities induced by the slow synaptic kinetics of in TRN can be inhibited by either the isolated SNr-TRN pathway or the isolated SNr-SRN pathway via different biophysical mechanisms. More importantly, our results showed that under certain conditions these two types of seizure control can coexist in the same network, suggesting that both decreasing and increasing the activation levels of SNr may considerably suppress the generation of SWDs. Theoretically, such bidirectional control of absence seizures by basal ganglia is due to the effective competition between the SNr-TRN and SNr-SRN pathways, which might be a generalized mechanism for regulating absence seizure activities and can be extended to other pathological factors. In addition, our detailed frequency analysis also indicated that, depending on different system conditions, the developed model may exist low, high or both TMFRs for the SNr neural population for triggering the typical 2–4 Hz SWDs.

These results are partly in agreement with former experimental observations. Previously, experimental studies based on electrophysiological recordings have established the linkage that reducing the activation of SNr neurons from the normal level can effectively suppress the SWDs in different rodent animal models [21]–[23], [26]. Such antiepileptic effect was supposed to be attributed to the indirect pathway of the SNr to TRN relaying at superior colliculus. The results presented in this work also demonstrated that decreasing the SNr activity is an effective approach that terminates the SWDs. However, it is important to note that in our model the similar antiepileptic effect is triggered by the direct GABAergic projections from SNr to TRN. Presumably, this is because the SNr has overall inhibitory impacts on TRN via both indirect and direct pathways. In the brain of absence epileptic patients, both of these two pathways might work together and provide a stable and endogenous mechanism to terminate the paroxysm of absence epilepsy.

Our model further makes prediction that increasing the activation of SNr from the normal level may also suppress SWDs. In previous experimental studies, there still lacks sufficient evidence to support this viewpoint. We speculate that this might be because the relative strength of the SNr-TRN and SNr-SRN pathways for rodent animals is generally high or at least not too low. Under this condition, the BGCT system for most rodent animals does not operate in the bidirectional control region, thus the suppression of SWDs caused by activating SNr is difficult to be observed in normal experiments. According to our results, the activation of SNr induced SWD suppression might appear by suitably tuning the relative strength between these two pathways. Further experiments based on animal models will be necessary to validate this prediction and characterize the detailed nature of the SWD suppression induced by the SNr-SRN pathway. Even so, the above prediction from our computational study might provide an alternative approach for terminating absence seizures.

An interesting and important question is: can real brain utilize this bidirectional modulation mechanism to control the paroxysm of absence epilepsy? Our answer is that it is possible, if the real brain has some mechanisms to automatically adjust the balance between the SNr-TRN and SNr-SRN pathways. Theoretically, there are several possible biological mechanisms and one of them is discussed as follows. Experimental data have uncovered that synapses conduct signals in an unreliable fashion, which is due to the probabilistic neurotransmitter release of synaptic vesicles [44]–[46]. It has been shown that the transmission failure rate at a given synapse generally tends to exceed the fraction of successful transmission, and in some specific cases it can be even higher than 0.9 [44]–[46]. Interestingly, recent studies indicated that such synaptic unreliability may play critical functional roles in neural processing and computation [47]–[49]. For patients with absence seizures, the suitable competition between the SNr-TRN and SNr-SRN pathways can be theoretically achieved by properly tuning the synaptic transmission rate of these pathways. This might be an important underlying mechanism and has significant functional advantages, because it does not require the changes of related anatomical structures and connection densities in the brain. However, we should notice that such tuning is not easy. Specifically, it requires the cooperation among neurons and needs to change the synaptic transmission rates collectively for most of relevant synapses. From the functional perspective, this mechanism may be associated with the self-protection ability of brain. After a long time of evolution, it is reasonable to suppose that our brain might have the abilities to use this type of plasticity-like mechanism to achieve complicated self-protection function during absence seizures [49]. Nevertheless, additional well-designed experiments are still needed to test whether our proposed hypothesis is correct.

The current results highlight the functional roles of basal ganglia in the control of absence seizures and might offer physiological implications on different aspects. First, the termination of absence seizures by activating SNr neurons through the SNr-SRN pathway might inspire the treatment of refractory absence seizures. Although the absence epilepsy is one typical benign epilepsy, a considerable proportion of patients yet may fail to achieve freedom from absence seizures and become refractory to multiple antiepileptic drugs [50], [51]. For those patients, one possible reason is that the strengths of their direct and indirect SNr-TRN pathways are somewhat weak, compared to the other patients. In this case, reducing the activation level of SNr may lose efficacy in the suppression of SWDs and, contrarily, increasing the firing of SNr neurons may stop the absence seizures. Second, our results might be generalized to the GPi neurons. In addition to the SNr, the GPi is also an important output structure of basal ganglia. The SNr often works in unison with GPi, since they have closely related inputs and outputs, as well as similarities in cytology and function [30]–[32]. It is thus reasonable to infer that the activation level of GPi might serve a similar bidirectional role in modulating the absence seizures. Third, the results presented in this study may provide new insight into the deep brain stimulation therapy on absence seizures [52]. In previous studies, it has been demonstrated that absence seizures can be inhibited by suitably applying the deep brain stimulation to STN. Our results not only support the traditional viewpoint, but also indicate that the SNr may be a more direct therapeutic target, compared to STN. We thus infer that the inhibition of typical absence seizure activities can also be accomplished by appropriately stimulating the SNr neurons. Finally, the main goal of this study was to explore the control role of basal ganglia in absence seizures, but not limited to this specific epilepsy. Several other types of epilepsies, such as the juvenile myoclonic epilepsy [53] and the generalized tonic-clonic epilepsy [54], are also highly associated with the corticothalamic system and mediated by GABA receptors. If our above findings on absence seizures capture the real fact, such bidirectional control feature by basal ganglia may be also available for these types of epilepsies.

Although our developed model is powerful enough to suggest some functional roles of basal ganglia in controlling and modulating absence seizures, we have to admit that this biophysical model is simplified and idealized. The limitations of this model and possible extensions in future studies need to be discussed. First, more detailed models based on the spiking neurons are able to introduce much more complexities, such as ionic dynamics, firing variability and connection property, to the system. In previous studies, it has been observed that ion concentrations inside the cell and in the extracellular space as well as the firing dynamics of neurons are changed during epileptic seizures [55]–[57]. It has been also found recently that the connection property of neural systems in different scales (for example, neuronal networks at microscopic scales and brain region networks at macroscopic scales) also contribute to the generation of seizure-like activities [58], [59]. Therefore, we cannot exclude other possible interesting results by using more detailed modelling methodologies. Moreover, due to lack of necessary data, our model does not include the indirect pathway from the SNr to TRN, relaying at the intermediate and superficial layers of the caudal superior colliculus [60]. Although we have inferred that both the direct and indirect SNr-TRN pathways play the similar effects on the control of absence seizures, further studies based on detailed anatomical data need to be investigated in the future work. Finally, it should be emphasized that in the present study we mainly focus on the SWDs induced by the slow kinetics of receptors in TRN. Even though we also showed that the similar bidirectional control behavior can be observed for another pathological factor (i.e., the corticothalamic loop transmission delay), more detailed analysis on this aspect is still needed. As we know, several other pathological factors, such as the increased T-type Ca2+ current in thalamocortical neurons [2], may also lead to the generation of SWDs. Nevertheless, it is still unknown and deserves to be further examined whether our proposed bidirectional control by basal ganglia is also available for these other pathological factors.

Theoretically, a better choice for mimicking the slower dynamics of currents should use relatively smaller values of and at the relevant synapses. However, we have to notice that the mean-filed theory used in the present study is under an assumption of uniform and parameters for all incoming connections of one neural population (see Eq. 3). To our knowledge, introducing different and to different types of connections will make such simple form of differential operator presented in Eq. 3 become invalid, thus making the dendritic filtering of incoming signals mathematically intractable in the current mean-field framework. Considering this, we employed a delay parameter to represent the order of magnitude of rise and decay time of currents in our model, as the method used in a recent neural mean-field modelling study [13]. Essentially, more detailed models based on spiking neurons will allow us to consider different dynamical processes for different types of synapses [49], [61]. We are currently trying to extend our above-identified results to a relevant spiking neural network and the corresponding results will be presented elsewhere.

It should be also pointed out that traditional modelling studies on absence seizures mainly focus on the possible generation mechanism of the typical absence seizure activities, i.e., SWDs, within the corticothalamic system, but little research addresses the regulatory mechanisms of absence seizures by basal ganglia. Inspired by previous experimental observations [21]–[23], [26] and recent modelling studies on Parkinson's disease [30], [31], here we introduced the basal ganglia to the traditional mean-field model of corticothalamic system proposed by Robinson and his collaborators, and firstly explored how the basal ganglia control the SWDs through the computational approach. Theoretically, several other models that describe the macro-scale dynamics about the activity of neural ensembles, such as the Wilson-Cowan model [62]–[64] and the models based on dynamic causal modelling [65]–[67], can be also used to deal with the similar problem but may need to re-estimate parameters according to different model assumptions.

In summary, we have performed mechanistic studies and investigated the detailed roles of basal ganglia in the control of absence seizures. We have computationally demonstrated that the SWDs generated by the developed model can be terminated and modulated by both decreasing and increasing the activation levels of SNr. Our results provide the first evidence on the bidirectional control of absence seizures by basal ganglia. This finding might deepen our conventional understanding about the functional roles of basal ganglia in the controlling and modulating of absence seizures. For patients with absence epilepsy, these results in turn indicate that the loss of several basal ganglia functions, especially the functions related to the SNr neurons, might further aggravate the attacks of absence epilepsy. We hope that predictions from our systematic model investigation can not only inspire testable hypotheses for future electrophysiological experiments but also provide additional therapeutic strategies for absence seizures.

Several other typical time series of (left frames) and their corresponding spectra (right frames) generated by our developed BGCT model. To a certain extent, these model time series are comparable with real physiological EEG signals: eyes-open (A) and (B), alpha rhythm (C), beta rhythm (D), coexistence of alpha and beta rhythms (E), and polyspike and wave (F). Here we adjust the excitatory corticothalamic coupling strength and delay parameter for reproducing these time series. The detailed parameter values used in our simulations are: and (A), and (B), and (C), and (D), and (E), and and (F), respectively. In addition, we also introduce a certain level of gaussian white noise into , with the mean and standard deviation . Note that more types of model time series might be observed by further tuning these critical parameter values.

(TIF)

Click here for additional data file.

An example representation of the bidirectional control of absence seizures by the basal ganglia in a relatively larger scale factor interval. The state analysis (A) and frequency analysis (B) in the () panel, with the inhibitory coupling strength . Here is the scale factor, and is the excitatory coupling strength of the STN-SNr pathway used to control the activation level of SNr neurons. Similar to the results in Figs. 5A and 6A, only three dynamical state regions are observed in the phase diagram (A): the SWD oscillation region (II), the simple oscillation region (III) and the low firing region (IV). In (A), the region marked by red diamond denotes the whole suppression regions of SWDs, the white dashed line represents the boundary of suppression region, and the red dashed line stands for the demarcation between the bidirectional (double arrow) and unidirectional (single arrow) suppression regions. In (B), the red asterisk region surrounded by dashed lines denotes the typical 2–4 Hz SWD oscillation region. Compared to the results shown in Figs. 5 and 6, here we consider a relatively larger scale factor interval from 0 to 2.3.

(TIF)

Click here for additional data file.

A series of two-dimensional state analysis in the ( ) panel for the modified model. In (A)–(H), we use the same group of parameter values but different random initial conditions for simulations. Similar to previous results, four different dynamical states are observed: the saturation state (I), the SWD oscillation state (II), the simple oscillation state (III) and the low firing state (IV). In each subfigure, the region between two red dashed lines denotes the main bidirectional suppression region of SWDs, where the double arrow represents that both increasing and decreasing can inhibit the generation of SWDs. The results given in (A)–(H) indicate that the modified model shows bistability (the simple oscillation state or the saturation state) in the large and strong region, and the final dynamics of the modified model significantly depend on the initial conditions. Note that in all simulations, we set , , , , and .

(TIF)

Click here for additional data file.

Supporting information code.

(TXT)

Click here for additional data file.

Supporting information code.

(TXT)

Click here for additional data file.