Control of Absence Seizures by the Thalamic Feed-Forward Inhibition

As a subtype of idiopathic generalized epilepsies, absence epilepsy is believed to be caused by pathological interactions within the corticothalamic (CT) system. Using a biophysical mean-field model of the CT system, we demonstrate here that the feed-forward inhibition (FFI) in thalamus, i.e., the pathway from the cerebral cortex (Ctx) to the thalamic reticular nucleus (TRN) and then to the specific relay nuclei (SRN) of thalamus that are also directly driven by the Ctx, may participate in controlling absence seizures. In particular, we show that increasing the excitatory Ctx-TRN coupling strength can significantly suppress typical electrical activities during absence seizures. Further, investigation demonstrates that the GABAA- and GABAB-mediated inhibitions in the TRN-SRN pathway perform combination roles in the regulation of absence seizures. Overall, these results may provide an insightful mechanistic understanding of how the thalamic FFI serves as an intrinsic regulator contributing to the control of absence seizures.


INTRODUCTION
Absence epilepsy is a common subtype of idiopathic generalized epilepsies and mainly occurs in the childhood years (Durón et al., 2005;Tolaymat et al., 2015). This chronic neurological brain disorder is characterized by recurrent absence seizures, often leading to sudden periods of impaired consciousness (Crunelli and Leresche, 2002). During its ictal period, the bilaterally synchronous 2-4 Hz spike and wave discharges (SWDs) can be observed on the electroencephalogram (EEG) of patients, which are regarded as the electrophysiological hallmark of absence seizures (Crunelli and Leresche, 2002). Despite that onset mechanisms of absence seizures still remain debated, accumulating evidence indicates that the emergence of 2-4 Hz SWDs during absence seizures critically depends on cerebral cortex (Ctx) and thalamus (Blumenfeld, 2005;Bazhenov et al., 2008;Lytton, 2008;Avoli, 2012;Holt and Netoff, 2013;Depaulis et al., 2016). It has been postulated that properly disrupting abnormal interactions within corticothalamic (CT) system can cause SWD suppression, thus preventing absence seizures.
Previous studies have shown that direct neuronal modulations within the CT system may destabilize the pathological oscillations during absence seizures (Suffczynski et al., 2004;Lüttjohann and van Luijtelaar, 2013;Paz et al., 2013;Taylor et al., 2014Taylor et al., , 2015Fan et al., 2015;Liu et al., 2016). The most frequently reported direct stimulation site is the thalamus. Past experimental evidence has suggested that both closed-loop optogenetic control of thalamus and high frequency thalamic stimulation can dramatically inhibit the occurrence of epileptic seizures (Lüttjohann and van Luijtelaar, 2013;Paz et al., 2013). In addition to these direct interventions, indirect neuronal modulations from several external brain regions have also been identified to be effective in controlling absence seizures Paz et al., 2007;Chen et al., 2014Chen et al., , 2015Kros et al., 2015a,b). As an intermediate bridge between the Ctx and thalamus, the basal ganglia (BG) are widely considered as an ideal modulation site to prevent absence seizures (Depaulis et al., 1994;Deransart et al., 1998;Paz et al., 2007;Luo et al., 2012). This hypothesis is supported by both experimental and computational evidence, showing that the BG might multiply control absence seizures through the GABAergic nigro-thalamic and pallidocortical pathways Paz et al., 2005Paz et al., , 2007Luo et al., 2012;Chen et al., 2014Chen et al., , 2015. Besides the BG, another ideal candidate to modulate abnormal CT oscillations is the cerebellar nuclei (CN; Kros et al., 2015a,b). Using a mouse model of absence epilepsy, a recent study has found that increasing the frequency and regularity of CN neuronal firing can significantly reduce the occurrence of SWDs and terminate absence seizures (Kros et al., 2015b).
Moreover, recent intensive statistical investigations on neuronal connectivity have revealed that the brain contains some recurring nontrivial patterns of interconnected neurons, termed as "microcircuit motifs" (Sporns and Kötter, 2004;Womelsdorf et al., 2014;Paz and Huguenard, 2015). Interestingly, these microcircuit motifs are found to be linked to each other in a way that does not spoil the independent function of each microcircuit motif, thus postulating to work as basic building blocks of nervous systems and perform critical functional roles in regulating neurodynamics at the circuit level (Sporns and Kötter, 2004;Li, 2008;Guo and Li, 2009;Potjans and Diesmann, 2014;Womelsdorf et al., 2014;Paz and Huguenard, 2015;Yuste, 2015). As one of the most significant microcircuit motifs, the feed-forward inhibition (FFI) has been observed in multiple brain regions, including the neocortex, thalamus, hippocampal, and basal ganglia (Sloviter, 1991;Sun et al., 2005;Assisi et al., 2007;Tepper et al., 2008;Kee et al., 2015;Paz and Huguenard, 2015;Khubieh et al., 2016;Lindahl and Hellgren Kotaleski, 2016). Emerging evidence has suggested that significant loss of FFI in the brain can cause an imbalance between excitation and inhibition of neurons toward excessive excitation (Paz and Huguenard, 2015). Such imbalanced excitation-inhibition might lead to hypersynchronous neural firing activities that are commonly involved in various types of epileptic seizures.
Neurons in the thalamic reticular nucleus (TRN) receive excitatory signals from the Ctx, and send both GABA A -and GABA B -mediated inhibitory signals to the specific relay nuclei (SRN) of thalamus, which are also driven by the excitatory signals from the Ctx. Anatomically, these neural projections together form the so-called thalamic FFI microcircuit motif (Pinault, 2004;Paz and Huguenard, 2015). Although past experimental studies have suggested that restoring cortical excitation in TRN neurons may prevent absence seizures (Paz et al., 2011;Paz and Huguenard, 2015), so far the detailed biophysical mechanisms of SWD suppression induced by this thalamic FFI is still poorly understood. In the present study, we utilize a meanfield CT network model to address this question. We find that both the excitatory Ctx-TRN coupling strength and the relative strength between the GABA B -and GABA A -mediated inhibitions in the TRN-SRN pathway play critical roles in preventing pathological 2-4 Hz SWDs generated by the CT system. These findings highlight the functional importance of thalamic FFI in controlling absence seizures, and might provide insights into the treatment of absence epilepsy.

Computational Model of the CT Network
Previous studies have suggested that the generation of 2-4 Hz SWDs during absence seizures may be caused by abnormal interactions between the Ctx and thalamus (Crunelli and Leresche, 2002;Robinson et al., 2002;Suffczynski et al., 2004;Blumenfeld, 2005;Breakspear et al., 2006;Lytton, 2008;Marten et al., 2009a). To investigate possible roles of the thalamic FFI in controlling absence seizures, we establish a macroscopic-level computational model that describes the population dynamics of the CT network. As shown in Figure 1, our CT network comprises four neural populations, which are indicated as follows: e = excitatory pyramidal neurons, i = inhibitory interneurons, r = TRN and s = SRN. Inspired by anatomical data, three types of neural projections are considered in our CT model (Figure 1). The blue solid and dashed lines with round heads are used to represent the inhibitory projections mediated by GABA A and GABA B synapses, respectively. The red lines with arrow heads are employed to denote the glutamate-mediated excitatory projections. Similar to previous modeling studies, a non-specific external input φ n is also injected to the SRN to mimic the sensory input (Robinson et al., 2002;Breakspear et al., 2006;Marten et al., 2009a;Chen et al., 2014Chen et al., , 2015. Obviously, our CT network contains the thalamic FFI that lies at the Ctx-TRN-SRN pathway (see the gray shaded pathways in Figure 1).
In the present study, each neural population in the CT network is simulated using the mean-field model proposed by Robinson and his colleagues (Robinson et al., 2002;Breakspear et al., 2006). For a given neural population a, the population dynamics can be characterized by three state variables: the mean membrane potential V a , the mean firing rate Q a , and the presynaptic activity φ a . The transformation from the mean membrane potential V a to the mean firing rate Q a obeys a sigmoid function (Robinson et al., 2002;Breakspear et al., 2006;Marten et al., 2009a;Freyer et al., 2011;Yang et al., 2016): where subscripts a ∈ A = {e, i, r, s} indicate the neural populations considered in the CT network. r represents the spatial position in the brain, θ a is the threshold of the mean firing, and σ denotes the standard deviation of the firing rate. Due to the boundedness of the sigmoid function, the maximum firing rate Q max a in the Equation (1) ensures the mean firing rate Q a Frontiers in Computational Neuroscience | www.frontiersin.org There are four neural populations that the excitatory pyramidal neural population (e), the inhibitory interneuron population (i), the TRN (r), and SRN (s) in the corticothalamic network. A non-specific external input to the SRN is denoted as φ n . The connections within the network are mainly mediated by excitatory glutamate described by red lines with arrow heads and inhibitory GABA denoted by blue lines with round heads. The solid and dashed lines from the TRN to SRN are used to distinguish the projections mediated by GABA A (v A sr ) and GABA B (v B sr ) synapses, respectively. The gray shaded pathways show the thalamic feed-forward inhibition (FFI).
within the physiological range as changing the mean membrane potential V a .
Theoretically, the filtered incoming postsynaptic potentials from presynaptic neural populations in the dendrites lead to fluctuations of the mean membrane potential V a at the position r. This dynamical process can be modeled as (Robinson et al., 2002;Breakspear et al., 2006;Marten et al., 2009a;Freyer et al., 2011): Here, the D αβ is a differential operator describing the synaptic and dendritic filtering of incoming signals. The inverses of α and β denote the decay and rise time constants of membrane potential caused by incoming signals, respectively. v ab describes the coupling strength from the neural population b to a, and φ b (r, t) represents the presynaptic activity of the neural population b. Similar to previous modeling studies (Chen et al., 2014(Chen et al., , 2015, the transmission delays among most neural populations are ignored in our model. Note that, however, we introduce a delay parameter τ to the GABA B -related incoming pulse rate [i.e., φ r (r, t − τ )] to mimic the slow kinetics of GABA B receptor-mediated inhibition.
In the Ctx, the propagation effects of excitatory firing rates Q e along the cortical surface satisfy the damped wave equation (Robinson et al., 2002;Breakspear et al., 2006;Marten et al., 2009a;van Albada and Robinson, 2009;van Albada et al., 2009;Freyer et al., 2011;Chen et al., 2014Chen et al., , 2015: where ∇ 2 represents the Laplacian operator (the second spatial derivative) and φ e denotes the electrical field arising from the cortical excitatory neurons firing propagating along axons. The temporal damping rate of pulses is determined by the γ e = v e /r e , where v e is the conduction velocity and r e denotes the characteristic range of axons of excitatory neurons in the Ctx. For other neural populations, their axons are too short to support wave propagation on the relevant scales. This indicates φ c = F(V c ) (c = i, r, s). On the other hand, because absence seizures are believed to occur simultaneously throughout the brain, it is reasonable to assume that the spatial activities during absence seizures are uniform in our CT model. We therefore ignore the spatial derivative in Equation (4), and the propagation effect of cortical excitatory axonal field can be rewritten as (Breakspear et al., 2006;Marten et al., 2009a;Freyer et al., 2011): We can further simplify our model by assuming that the mean membrane potential V i and mean firing rate Q i of the inhibitory interneuron population satisfy V i = V e and Q i = Q e . Note that this assumption is reasonable, because the intracortical connections are suggested to be proportional to the number of synapses involved (Robinson et al., 2002(Robinson et al., , 2004Breakspear et al., 2006;Marten et al., 2009a;van Albada and Robinson, 2009;van Albada et al., 2009;Freyer et al., 2011;Chen et al., 2014Chen et al., , 2015. Consequently, we can rewrite our mean-field CT model in the first-order form, as given in the Appendix. It should be noted that the model parameters in the established CT network are adapted from previous studies (Robinson et al., 2002(Robinson et al., , 2004Breakspear et al., 2006;Marten et al., 2009a;van Albada and Robinson, 2009;van Albada et al., 2009;Freyer et al., 2011;Chen et al., 2014Chen et al., , 2015, which are systematically summarized in Table 1. These default parameters are estimated from real physiological data. Unless stated otherwise, we use these default parameter values in the following studies.
In simulations, our computational model is implemented with custom codes written in Matlab (MathWorks). We use the standard fourth-order Runge-Kutta method for numerically integrating normal differential equations, and employ an algorithm adapted from Matlab dde23 for calculating the delay differential equations (Shampine and Thompson, 2001). The temporal resolution of integration is fixed at 0.05 ms. Additional simulations confirm that this time step is sufficiently small to Model parameters employed in the present study are adapted from previous studies (Chen et al., 2014(Chen et al., , 2015 and references cited therein. Unless otherwise noted, we use these default parameter values for simulations. Note that PY and IIN denote excitatory pyramidal neurons and inhibitory interneurons in the cerebral cortex, respectively.
ensure an accurate simulation of our established CT model. All simulations are carried out for sufficiently long time to collect stable data for further analysis. The computer codes implementing our CT model will be available upon request via email.

Data Analysis
In the present study, we analyze the simulated data generated by our model from various aspects. Note that these data analysis techniques have been described in our previous studies (Chen et al., 2014(Chen et al., , 2015. To classify the dynamical states of cortical oscillations, we determine the stable local maximum and minimum values of cortical excitatory axonal fields occurring in one cortical oscillatory period for each experimental setting. For a specific parameter, the state analysis of cortical oscillations is plotted as a function of this critical parameter (for example, see Figure 3). By employing this state analysis technique, we can identify the regions of different dynamical states in the two-dimensional parameter space (for example, see Figure 2).
To evaluate the oscillation frequency of cortical oscillations, we compute the power spectral density of the time series φ e using the Fast Fourier transform. For an experimental setting, the maximum peak frequency emerging in the power spectral density is defined as the dominant frequency of cortical oscillations. By combining both the state and frequency analysis techniques, the boundaries of the typical 2-4 Hz SWDs regions can be further outlined in the two-dimensional parameter space (see the asterisk regions in Figure 2). To understand the biophysical mechanisms underlying the control of absence seizures by the thalamic FFI, we calculate the long-term mean firing rates (MFRs) of three populations that the Ctx, the TRN, and SRN. To quantify impacts of Ctx and SRN on TRN, we further compute the mean potential that from Ctx and SRN to TRN, denoted as V Ctx-TRN and V SRN-TRN in Figure 3D. In this work, V Ctx-TRN is defined as the product of the v re and the MFRs of Ctx, whereas V SRN-TRN represents the product of the v rs and the MFRs of SRN. The summation of V Ctx-TRN and V SRN-TRN is thus defined as the overall effect of Ctx and SRN on TRN (V TRN ).

Our CT Model Can Generate Typical 2-4 Hz SWDs by Introducing Pathological Mechanisms
Accumulating experimental data reveal that the generation of absence seizures is tightly associated with abnormal interactions within the CT network (Snead, 1995;Danober et al., 1998). Specifically, both the strong excitatory projections from the Ctx to the SRN and the slow kinetics of GABA B synapses along the TRN-SRN pathway are believed to be critical pathological factors in absence seizures (Destexhe, 1998;Robinson et al., 2002;Breakspear et al., 2006;Marten et al., 2009a,b;Han et al., 2012). Theoretically, the strong coupling strength of the corticothalamic pathway tends to induce the high amplitude oscillations, and the slow kinetics of GABA B synapses act as a delayed low-pass filter and can be linked to the production of slow waves. Thus, simultaneous modulation of these two parameters may destroy the normal oscillatory pattern generated in the CT system, triggering the onset of absence seizures.
To examine whether these pathological factors are also available in our CT model, we employ two related parameters: the excitatory coupling strength v se and the GABA B delay τ . Figures 2A,B show the state and frequency analysis in the combined parameter space (v se , τ), respectively. Consistent with previous findings (Chen et al., 2014(Chen et al., , 2015, the model exhibits four types of dynamical states with different frequency characteristics. The first dynamical state is the saturation firing that occurs for sufficiently long τ and strong v se (region I in Figure 2A). Under this condition, inhibitions from the TRN to SRN cannot successfully suppress the excessive excitation caused by strong cortico-thalamic interactions. Such excessive excitation drives the firing of cortical neurons to increase from the low firing to the maximum firing in a short time (Figure 2C). For a long GABA B delay τ , the SWD oscillation state appears provided that the excitatory coupling strength v se is appropriately strong. In this case, the TRN-SRN inhibitions mediated by GABA A and GABA B receptors effectively shape SRN firing at different instants. This dual suppression creates double firing peaks for SRN neurons, further leading to the generation of cortical SWDs (region II in Figure 2A). Importantly, most of SWD oscillation region is contained in the typical 2-4 Hz frequency range (asterisk regions in Figures 2A,B), which can be comparable with EEG signals of real patients during absence seizures. On the other hand, the decrease in delay τ mixes the GABA A -and GABA B -induced inhibitions, weakening their double peak shaping effect. For a fixed v se , we therefore observe that the network dynamics can transit from the SWD oscillation state to the simple oscillation state by lowing the GABA B delay (region III in Figure 2A). Theoretically, the boundary between dual time scales and mixed time scale highly depends on the decay and rise time constants of membrane potential in our model (Chen et al., 2014). For each GABA A -induced inhibition caused by TRN neurons, the SRN neurons require a certain level of recovery time to restore their neuronal firing to the rising state. If this recovery time is shorter than the GABA B delay, another firing peak can be created to SRN neurons due to the latter GABA B -induced inhibition. Further, frequency analysis reveals that the dominant frequency of simple oscillation state relies on the GABA B delay. The shorter the delay parameter τ , the faster the simple oscillations generated by our model (Figures 2E,F). However, if the excitatory coupling strength v se is too weak, neural oscillations within the CT network cannot be maintained. In this case, the network dynamics is kicked into the low firing state (region IV in Figures 2A,G).
These findings demonstrate that our model can replicate different dynamical states of human brain. Remarkably, the model successfully generates the typical 2-4 Hz SWDs by introducing previous pathological mechanisms. In the following studies, we set v se = 2.4 mV s and τ = 50 ms as default, and explore the effects of the thalamic FFI on controlling absence seizures.

Roles of the Ctx-TRN Pathway in Controlling Absence Seizures
In the brain, the FFI microcircuits are believed to serve as fundamental regulators in balancing neural excitation and inhibition. It is widely accepted that dysfunction of FFI microcircuits might lead to excessive and/or hypersynchronous neural activities and evoke epileptic seizures (Paz and Huguenard, 2015). Importantly, a recent experimental study also showed a specific reduction in Ctx-TRN strength in the GluA4-deficient (Gria4 −/− ) mice model of absence epilepsy (Paz et al., 2011). We thus hypothesize that the classical thalamic FFI microcircuit, i.e., the Ctx-TRN-SRN pathway (Figure 1, gray lines), may participate in the modulation of absence seizures.
To theoretically validate whether this hypothesis is correct, we concentrate on the TRN and modulate its activation level by changing the coupling strength of the excitatory Ctx-TRN pathway. Figure 3A shows the one-dimensional state analysis as a function of the Ctx-TRN strength v re . Indeed, we find that the Ctx-TRN strength performs an active role in terminating absence seizures. A slight enhancement of the Ctx-TRN strength from its default value can suppress the SWD oscillations and push the model dynamics into the simple oscillation state (Figure 3A). For sufficiently strong strength v re , the model dynamics is even kicked into the low firing state (Figure 3A), implying a strong suppression effect due to the strengthened Ctx-TRN pathway. To characterize the frequency feature of model dynamics, the dominant frequency of neural oscillations is estimated by using the spectral analysis technique. For equivalent GABA A -and GABA B -mediated inhibitions in the TRN-SRN pathway, we observe that the dominant frequency of neural oscillations is stable and insensitive to the excitatory Ctx-TRN strength ( Figure 3B).
To understand in detail how the Ctx-TRN pathway inducing SWD suppression arises, we further plot the MFR vs. the excitatory Ctx-TRN pathway strength v re for key neural populations in Figure 3C. As expected, the enhancement of excitatory coupling strength v re tends to excite TRN neurons ( Figure 3D, green, indicted by diamonds). However, increasing the TRN activation suppresses the firing of SRN, which subsequently inhibits cortical neurons. Conversely, however, inactivation of SRN neurons will introduce a collision to TRN neurons and reduce their firing in a significant way (Figure 3D, black and red, indicated by circles and squares, respectively). By strengthening the Ctx-TRN pathway, this collision destabilizes the pathological balance and causes an overall inhibition effect in the CT system ( Figure 3C). Such inhibition effect on TRN weakens the double peak shaping effect due to the slow kinetics of GABA B receptors. For a relatively strong excitatory Ctx-TRN strength, this GABA B weakening effect is considerable and therefore the generation of SWDs is terminated in our model.
Past studies have demonstrated that several other pathological factors may also lead to typical 2-4 Hz SWD oscillations in the CT system (Robinson et al., 2002;Breakspear et al., 2006). Thus, a natural question is whether such control feature caused by the thalamic FFI is a generalized regulatory mechanism for absence seizures. We try to answer this question by introducing another SWD generation mechanism into our model. In literature, it has been widely reported that long transmission delay between the Ctx and thalamus can drive the CT system to produce the 2-4 Hz SWDs. To apply this pathological factor in our model, we block the GABA B pathway from TRN to SRN, and consider a bidirectional transmission delay τ 0 /2 between the Ctx and thalamus. Note that several coupling strengths are also adapted to ensure the CT system to generate the SWD oscillation. Our results presented in Figure 4 confirm that the increase in the excitatory coupling strength v re can successfully suppress the typical 2-4 Hz SWDs due to abnormal transmission delay between the Ctx and thalamus, further emphasizing the generality and functional importance of the thalamic FFI in regulating absence seizures.
Taken together, we show that the thalamic FFI indeed participates in the modulation of absence seizures. Our theoretical observations indicate that normal CT oscillations require a certain level of Ctx-TRN strength, and loss of the Ctx-TRN strength might lead to over-excitation in the CT system and thus cause absence seizures. Importantly, this control manner by the thalamic FFI is found to be a possibly generalized regulatory mechanism for absence seizures and may be extendable to other pathological factors.

Combination Roles of GABA A -and GABA B -Mediated Inhibitions in the TRN-SRN Pathway in Controlling Absence Seizures
The TRN-SRN pathway is another component of the thalamic FFI and might therefore also play essential roles in controlling absence seizures. Because neural projections from the TRN to the SRN are mediated by both GABA A and GABA B receptors, we next investigate effects of these two types of inhibitions on the modulation of absence seizures. Figure 5A depicts the two-dimensional state analysis in the Similarly, we identify that four types of dynamical states are distributed in these two-parameter space (Figure 5A), revealing that the dynamics of the CT system depends on both GABA A and GABA B inhibitions. At a low level of GABA A inhibition, the CT system exhibits the saturation firing state (Figure 5A, region I). For a relatively larger GABA B delay that we considered in this work (Destexhe, 1998;Marten et al., 2009a,b), this occurs even when the GABA B inhibition is at a sufficiently high level, indicating that strong disinhibition of GABA A receptors in the TRN-SRN pathway may cause explosion of high-rate activities (Chen et al., 2014). Note that, however, such effect caused by GABA A inhibition might be also dependent on the delay time in the GABA B inhibition. If both GABA A and GABA B inhibitions are too strong, the firing of SRN neurons is extremely inhibited. In this case, the network dynamics is pushed into the low firing state (Figure 5A, region IV). As a consequence, neural oscillations of the CT system can be observed for intermediate GABA A and GABA B inhibitions (Figure 5A, regions II and III). Moreover, we find that the emergence of SWDs is highly associated with a scale factor K, defined as the relative strength of GABA B and GABA A inhibitions in the TRN-SRN pathway. For a fixed GABA A strength −v A sr , an intermediate K might cause SWDs in the CT system (Figure 5C), thus might serve as a possible pathological factor of absence seizures.
To examine the dependence of neural oscillations on GABA A and GABA B inhibitions, we further perform dominant frequency analysis in the (−v B sr , −v A sr ) panel (see Figure 5B). In comparison with the GABA A inhibition, our results suggest that GABA Bmediated inhibition plays a dominant role in determining the frequency of neural oscillations generated by the CT system. At a given GABA A inhibition level, increasing the GABA B strength in the TRN-SRN pathway enhances and prolongs the firing depression of SRN neurons. Such strengthened suppression effect leads to a significant reduction in the frequency of neural oscillations. Accordingly, we observe the fast and slow simple neural oscillations at the low and high levels of GABA B inhibition, and the SWDs with frequency varying from 2 to 6.5 Hz at the intermediate level of GABA B inhibition (Figures 5A,B). Consistent with a previous modeling study (Destexhe, 1999), our frequency analysis suggests that the different oscillatory frequencies in the CT system may be mainly attributable to differences in the kinetics of GABA B -mediated inhibition in the TRN-SRN pathway.

DISCUSSIONS
Using a mean-field CT network model, we systematically investigated the roles of thalamic FFI in controlling absence seizures in this study. By computational modeling, we demonstrated that the enhancement of the excitatory Ctx-TRN strength can significantly inhibit typical 2-4 Hz SWDs during absence seizures. Theoretically, this SWD suppression is found to be due to the GABA B weakening effect caused by the collision in TRN neurons. More importantly, we showed that such control manner by the thalamic FFI might be a generalized regulatory mechanism for absence seizures and could be extendable to other pathological factors. Our findings highlight the functional importance of thalamic FFI, which might contribute to the treatment of absence epilepsy.
It has been experimentally reported that loss of the FFI both in thalamus and neocortex might lead to over-excitation of neurons and induce epileptic seizures (Paz and Huguenard, 2015). For example, a recent study has identified a pathophysiological mechanism for the generation of absence seizures in the Gria4 −/− mice model of absence epilepsy, for which the TRN receive reduced excitatory inputs from the direct Ctx-TRN neural projections (Paz et al., 2011;Paz and Huguenard, 2015). To a certain extent, this finding is in agreement with our computational observations, but our results established more insightful mechanistic understandings. Similarly, conditionally ablating Ca v 2.1 channel function in cortical parvalbumin (PV) + interneurons was found to cause a loss of the FFI in neocortex and result in the generalized absence seizures (Rossignol et al., 2013). In addition, a previous experimental study also showed that selective loss of the FFI in rat somatosensory barrel circuits induce epileptic seizures when thalamocortical afferents are activated (Sun et al., 2005). These experimental findings provide clear evidence that the FFI in thalamus and neocortex is of great importance to maintain the balance of excitation and inhibition, which might provide some inspirations to make anti-epileptic drugs (Paz and Huguenard, 2015).
From the anatomical perspective, the GABAergic inhibitions from the TRN to SRN are mediated by both GABA A and GABA B synapses. Several previous experimental studies have suggested that GABAergic transmission within thalamus might play critical roles in shaping abnormal CT oscillations during absence seizures (Liu et al., 1991(Liu et al., , 1992Snead, 1992), but so far the relevant biophysical mechanisms are still not completely established. Our computational findings presented here indicated that GABA A and GABA B -mediated inhibitions in the TRN-SRN pathway might play combination roles in regulating absence seizures. For an intermediate relative strength of GABA B and GABA A inhibitions, the double suppression caused by these two types of GABA receptors occurring at different time instants are effective and may provide an effective mechanism to create multiple firing patterns for different neural populations. Mechanistically, such double suppression effect is responsible for the onset of absence seizures in the CT system. Remarkably, our model further makes prediction that the dominant frequency of SWDs in the CT system is more sensitive to GABA B -mediated inhibition, and a high level of GABA B inhibition corresponds to a relatively slower oscillation frequency. This observation emphasizes that GABAergic transmission in the thalamic FFI serves as an important intrinsic regulator to modulate the oscillation frequency of SWDs.
The results presented in this work support the idea that significant loss of thalamic FFI might trigger absence seizures (Paz et al., 2011). However, it is worth to notice that the functional roles of the FFI in brain circuits might be diverse and brain region-dependent. In comparison with our current study that enhancing the thalamic FFI could terminate absence seizures, a recent computational study using a basal gangliathalamo-cortical network suggested that abnormally strong striatal FFI facilitates to maintain, but not interrupt, synchronous oscillations during absence seizures (Arakaki et al., 2016). Instead, the reduction in the striatal FFI activates medium spiny neurons in the striatum, which might provide potential roles to desynchronize abnormal oscillations and control absence seizures (Arakaki et al., 2016). Moreover, the FFI in neural circuits is not always acting as an opposite depressing the output excitation. Indeed, the FFI might enhance output excitation in some circumstances. For example, as a subset of GABAergic cortical interneurons, the chandelier cells that receive excitatory projections from pyramidal cells. However, it has been observed that these types of interneurons can provide excitatory inputs to their subsequent pyramidal cells (Molnár et al., 2008). Therefore, the functional roles of FFI might highly depend on the involvement of neural circuits and brain regions. Although here we focus on the functional role of thalamic FFI in suppressing absence seizures, it should be noted that the feedback from SRN to TRN might also impact our model dynamics. In additional simulations, we have shown that decreasing the coupling strength of SRN-TRN pathway could push the cortical firing pattern from the SWD oscillation state to the saturation firing state, which is a non-physiological brain state (see Supplementary Figure 1). Further, investigation is needed to uncover the importance of SRN-TRN pathway in controlling absence seizures.
Although our model is designed to capture basic modulations of absence seizures by thalamic FFI, it can be extended to investigate other related issues. For example, deep brain stimulation (DBS) has been demonstrated to be a powerful neurosurgical procedure to treat epileptic seizures (Lega et al., 2010;Fisher and Velasco, 2014;Pantoja-Jiménez et al., 2014). Recent experimental studies have shown that applying highfrequency DBS with suitable current intensity to both SRN and TRN can significantly reduce the occurrence of absence seizures (Lüttjohann and van Luijtelaar, 2013;Pantoja-Jiménez et al., 2014). Nevertheless, so far the most effective thalamic DBS target for treating absence seizures is still debated. In future studies, it is of importance and deserves to use our developed model to determine the most effective thalamic DBS target as well as its corresponding optimal stimulation parameters. On the other hand, several other mental illness, such as the generalized tonic-clonic epilepsy, Parkinson's disease, neurogenic pain and depression, are also reported to be highly associated with the CT dysrhythmia (Llinás et al., 1999). Theoretically, the CT model established in this work might provide a modeling framework for further exploring underlying functional roles of thalamic FFI in these mental illnesses. In addition, considering that suppression of absence seizures in the CT model is a dynamical control process, we will further try to investigate how to control this system from the view of dynamical control in our following work.
To summarize, we have performed mechanistic studies to examine underlying effects of thalamic FFI on the control and modulation of absence seizures. Our results showed that the FFI in thalamus not only contributes to suppressing absence seizures, but also serves as an intrinsic regulator to modulate the oscillation frequency of SWDs generated in the CT system. These observations established an insightful mechanistic understanding on how the thalamic FFI regulates typical 2-4 Hz SWDs during absence seizures, and might provide insightful physiological implications into the treatment of absence seizures.

AUTHOR CONTRIBUTIONS
MC, DG, YX, and DY conceived and designed the experiments. MC performed the experiments and data analysis. MC and DG wrote the paper.

FUNDING
This work was supported by National Natural Science Foundation of China (No. 81571770, No. 81330032, No. 61527815, No. 81371636