In this paper, we propose a modified thalamocortical neural field model originally from Taylor et al. (2014), where we introduce the feedforward inhibition from TC to IN into the thalamocortical network. As seen in Figure 2, connection schematic of this modified neural field model is composed of the thalamocortical loop, which consists of the cortical EX-IN circuit and thalamic TC-RE circuit. The current improved model equations can be written as follows:

Each equation is a rate equation with the units of voltage (activity) per time, where EX, IN, TC, RE are firing rates, i.e., the fractional firing activity in each neuronal populations, the parameters h_e,i,t,r are additive input constants to EX, IN, TC and RE, respectively. τ_e,i,t,r are the timescale parameters having the unit of s⁻¹. C_{ee,ei,et,ie,it,te,tr,re,rt,rr} are the connectivity strengths within different neuronal populations whose linking rules essentially follow the experimentally known connection values. In particular, C_it and C_et represent the feedforward inhibition and excitation from TC to cortex. f[·] and g[·] are the activation functions of the cortical subsystem and thalamic subsystem, respectively:

with x = EX, IN, TC and RE, y = TC, RE, and θ determines the steepness.

The thalamocortical model can be represented as follows (Rodrigues et al., 2007, 2009):

where F = (F₁, F₂, F₃, F₄)∈C^k(ℝ⁴ × ℝ²¹, ℝ⁴) is the vector field for large enough k, X = (EX, IN, TC, RE), ν∈ℝ²¹ is the total parameter space, but in this present study, the parameters of interest are the feedforward excitation and inhibition of cortico-thalamic synaptic circuit, w = (C_et, C_it) ⊂ ν, due to its physiological relevance in studies of absence seizures. Because C_et, C_it are the only parameters we allow to vary, in fact, F ∈C^k(ℝ⁴ × ℝ², ℝ⁴) and the system can be rewritten as

Let X^*(t) = (EX^*, IN^*, TC^*, RE^*) be fixed point of the system (8), i.e., equilibrium state being determined by setting the RHS (Borisyuk and Kirillov, 1992) of vector field of system (8) equal to 0,

The stability of equilibrium points can be analyzed by ensuring that the linearized version of (8) satisfies the Hartman-Groβman theorem (Kuznetsov et al., 1998). Consequently, we consider the Jacobian matrix (linearization matrix) of (8):

Assume that for Cet=Cet* and Cit=Cit* at least one of the eigenvalues of the Jacobian matrix calculated at X^*(t) is zero. This is valid if and only if the matrix determinant Δ(X*(t),w*)=|J|X*(t),Cet*,Cit*=0. We then consider the following equations

which defines a curve in the six-dimensional space (EX, IN, TC, RE, C_et, C_it). The projection of this curve onto the plane w = (C_et, C_it) is called the fold point bifurcation curve.

In the 4-dimensional system of differential equations, the necessary and sufficient condition for a sum of two eigenvalues to be zero is that 3th Gurvitz determinant Δ_G(X(t), C_et, C_it) = 0 (Gantmakher, 1967). If there exists a pair of eigenvalues of JX*(t), λ₁ and λ₂, the sum of which is zero (i.e., either λ_{1, 2} = ±iω, or λ₁= −λ₂ are real), we consider the following equations,

which define a curve in the six-dimension space (EX, IN, TC, RE, C_et, C_it). We call the projection of the curve onto the plane w = (C_et, C_it) the Andronov-Hopf bifurcation curve.

The double limit cycle curve can be defined by considering a fixed point of the Poincare map. Here the condition that the point (EX, IN, TC, RE, C_et, C_it) belongs to the double limit cycle curve can be written as a system of five equations, where four equations define that the point is the fixed point of Poincare map, and the fifth one defines that the point is a multiple fixed point of this map. The projection of this curve onto the plane w = (C_et, C_it) is called the double limit cycle curve. When parameters pass through this curve two limit cycles appear or disappear.

In addition, setting the linear change of variables to X^=X-X* and substituting appropriately into equations we obtain the following X^˙=F(X^,w) where now F(0,w) = 0. Therefore, without loss of generalization, the state X = 0 is chosen to be the state of low level background activity. That is X = 0 is the steady state solutions of (1–4). h_e,i,t,r are external inputs into EX, IN, TC and RE neuronal populations, which can qualitatively induce the transitions between pathological state and non-seizure state of system. In the presence of external input (i.e., h_e,i,t,r ≠ 0), substituting X = 0 into the RHS of (1–4), we can quantitatively identify h_e,i,t,r as follows (Equation 13) through the interactive intensities among different neuronal populations, C_{ee,ei,et,ie,it,te,tr,re,rt,rr}, within the thalamocortical circuit,

Finally, epileptic absence seizures are characterized by SWD which is essentially the bursting oscillations. Bursting is one of the most important firing activities of neuronal systems. Rinzel (1985) firstly made some theoretical analysis on bursting and recognized that bursting exhibits a transition between a slow wave (or resting state) and a spiking state, owning to slow variation process modulating fast firing activities. In this system, TC can be considered as the slow variables in the model subsystem (3), which controls the dynamics of the fast subsystem (1), (2), and (4). Therefore, the fast-slow dynamics analysis plays an important role for the dynamical behavior of firing patterns (Rinzel, 1985; Izhikevich, 2000; Duan et al., 2008), especially bursting-like SWD discharges of absence seizure in epilepsy.

All model parameters used to make simulations in this manuscript are listed in Table 1. The standard fourth-order Runge-Kutta integration scheme under the MATLAB (MathWorks, USA) simulating environment was used with Equations (1–4). The model EEG is taken as the mean of two cortical populations, EX and IN. Both the bifurcation and frequency analysis (Chen et al., 2014, 2015, 2017) are utilized to characterize the critical state transitions and neural oscillations generated by our model. Firstly, the bifurcation analysis is performed by calculating the stable local minimum and maximum values of the mean for EX and IN. Furthermore, to evaluate the dominant frequency of neural oscillations, the power spectral density (PSD) is estimated using the fast Fourier transform (FFT) for the time series of the mean for EX and IN. Then the maximum peak frequency is defined as the dominant frequency of neural oscillations. We further investigated the interesting dynamical features of the modified model numerically using the software package XPPAut (Ermentrout, 2002). The bifurcation diagrams were computed by using AUTO as incorporated in XPPAut.

Tonic-absence seizures, originally termed by Shih and Hirsch (2003), consist of a tonic seizure with generalized paroxysmal fast activity (>13 Hz) which is followed by an absence seizure with generalized slow spike-wave activity (~2–4 Hz). Clinically and electrographically, this is a poorly characterized seizure type in patients with symptomatic generalized epilepsy. The typical tonic-absence EEG pattern can be found in the Figure 1A, modified from Figure 2A in the work of Shih and Hirsch (2003). In this section, we computationally investigate the effect of feedforward excitation from TC to EX, i.e., C_et, on the generation of tonic-absence seizures, based on the modified thalamocortical network model.

In order to investigate the modulation functions of feedforward inhibition on the feedforward excitation inducing tonic-absence seizures, we gradually increase the feedforward inhibition, i.e., C_it, to observe the state transitions from the various phases of tonic-absence seizures. Figure 7 captures the typical state transitions under the different functions of feedforward inhibitions with fixing C_it at 0.1 (Figure 7A), 0.2 (Figure 7B), 0.4 (Figure 7C), 0.7 (Figure 7D), 1.0 (Figure 7E), and 1.8 (Figure 7F), respectively. In what follows, for the sake of clearness, the investigations for the state transitions are mainly focused on the three parameters intervals of C_et, i.e., C_et ∈ [0, 0.2], C_et ∈ [0.2, 0.6] and C_et ∈ [0.6, 2], respectively.

In order to observe the overall dynamics transition of the modified Taylor model as shown in Figure 2, we depict the state transitions (Figure 8A) and corresponding dominant frequency (Figure 8B) as the feedforward excitation and inhibition from subcortical thalamus to cortex, C_et and C_it, are changed in 2-D plane [0, 2] × [0, 2]. It can be seen in Figure 8A that there are ten types of firing states, which are denoted as: LS/HS: low/high saturated firings, TO: tonic oscillations with high frequency (>10 Hz) and low amplitude, h/l-CO: clonic oscillations with high (>5 Hz & <10 Hz) and low (~2 Hz) frequencies, r-CO: reversed clonic oscillations, m-SWD: multi-spike and wave discharges (m = 2,3), SWD: spike and wave discharges (m = 1), r-SWD: reversed SWD discharges.

Specifically, for the both small values of C_et and C_it, the system exhibits the epileptic tonic-clonic oscillations. By contrast, the much large C_et and C_it can suppress the system to show low-frequency clonic oscillations. However, the unbalanced functions between feedforward excitation and inhibition can induce the epileptic spike-wave discharges. For example, for the small C_it (e.g., <1, i.e., the weak feedforward inhibition), the strong feedforward excitation (e.g., C_et > 1) can activate the system to display typical m-SWD (m = 2, 3) and SWD (m = 1) discharges of absence seizures in epilepsy. Conversely, for the small C_et (e.g., <0.6, i.e., the weak feedforward excitation), the stronger feedforward inhibition can induce the reversed SWD (r-SWD) discharges and clonic (r-CO) oscillations. In particular, the horizontal double black arrow and the left vertical double pink arrow represent the physiologically observed epileptic transitions from tonic-clonic oscillations to absence seizures (see Figure 1A), which corresponds to the results of Sections 3.1.1 and 3.2.2, respectively. The right vertical pink arrow represents the physiologically observed seizure transitions from absence seizures to clonic oscillations of epilepsy (see Figure 1B), which also corresponds to the result of Section 3.2.2. In addition, non-oscillation activities, i.e., saturated states, occur for the small C_et and large C_it, or the large C_et and small C_it, corresponding to regions LS and HS, respectively. What's more, as shown in Figure 8B, only both weak feedforward excitation and inhibition, corresponding the parameter region of small C_et and C_it, can result in the high-frequency oscillations, while the other parameter regions correspond to the low-frequency oscillations.

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.