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., 2014, 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 Qamax in the Equation (1) ensures the mean firing rate Q_a 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, 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., 2014, 2015):

where ∇² 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, 2004; 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., 2014, 2015; Hu 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, 2004; 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., 2014, 2015; Hu 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 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.

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, 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).

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 cortico-thalamic 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, 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.

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 τ₀/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.

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 (-vsrB,-vsrA) panel. 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 -vsrA, 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 (-vsrB,-vsrA) panel (see Figure 5B). In comparison with the GABA_A inhibition, our results suggest that GABA_B-mediated 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.

These results indicate that GABA_A- and GABA_B-mediated inhibitions in the TRN-SRN pathway might play combination roles in controlling absence seizures. An intermediate relative strength of GABA_B and GABA_A inhibitions might trigger the onset of absence seizures in the CT system, and the dominant frequency of SWDs highly depends on the level of GABA_B inhibition.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.