Critical transitions in a model of a genetic regulatory system

1.   Introduction

Epilepsy is one of the most common disorders of the central nervous system, characterized by spontaneity, recurrence, and complex pathogenesis, seriously affecting brain function and even threatening life. Epilepsy develops when there is a broad activity in both hemispheres of the brain ^[1]. Eplieptic seizures can be divided into focal and generalized onset seizures. The types of generalized seizures in epilepsy include tonic, clonic, and absence seizures. In tonic-clonic seizures, a tonic activity is followed by a clonic activity. The two-way transition between absence seizures and tonic-clonic seizures has been observed through electrophysiological studies ^[2]. Dysfunction of the cortical or thalamic circuitries is believed to be the cause of the absence, clonic, and tonic seizures. The presence of SWDs or poly-spike and wave discharges (poly-SWDs) with a frequency range of 2–4 Hz is a common feature of typical absence seizures in electroencephalogram (EEG) recordings ^[3,4], while atypical absence seizures are characterized by a lower frequency of approximately 1–2.5 Hz ^[5]. Tonic seizures are a low amplitude, high frequency ($ > $14 Hz) fast-spiking activity, and clonic seizures are a high amplitude, low frequency (approximately 3 Hz) synchronized simple slow-wave activity. Although the pathogenesis of epilepsy remains controversial, studies have confirmed that seizures are closely linked to the cerebral cortex and thalamus ^[6,7]. Meeren et al. found that the cerebral cortex is the main trigger for sudden oscillations in corticothalamic circuits ^[8]. Avanzini and Sysoeva et al. found that seizures may be associated with the thalamus and there is a mutual coupling between the cortex and thalamus during epileptic seizures ^[9,10]. Thus, abnormal interactions between the cortex and the thalamus may be one of the main causes of absence seizures ^[11,12].

In general, an imbalance between excitatory and inhibitory functions is one of the triggers of epilepsy. Interneurons are commonly known as inhibitory neurons, but there are also excitatory neurons, so there are actually both inhibitory and excitatory interneurons in the central nervous system. Inhibitory interneurons use gamma-aminobutyric acid (GABA) as their neurotransmitter, while excitatory interneurons are spiny stellate cells found in the neocortex of the human brain and use glutamate as their neurotransmitter ^[13]. Experiments conducted on monkeys in the past have demonstrated that epileptic activity can be propagated in the brain through synchronous discharges of EIN. A study on the human brain has shown that the neurons in the neocortex, i.e., excitatory interneurons, play a crucial role in the production of SWDs ^[14,15]. In a study of the last two years, Tabatabaee et al. first investigated the key role of spiny stellate cells or excitatory interneurons in seizures ^[16], but didn't consider the excitatory effect of EIN on PY. Yan et al. ^[17,18] improved the thalamocortical model by introducing EIN to investigate the dominant role of EIN in the generation of SWDs. However, it is worth pointing out that excitatory interneurons have been less studied in thalamocortical models.

Recently, neural field models have gained popularity in the field of modeling. These models aim to simplify complex systems and represent the development of neuronal populations using minimal equations. In models of epileptic patients, thalamic-cortical interactions have been shown to play a key role in the production of SWDs ^[19,20]. Taylor et al. developed a neural field model of cortico-thalamic loop based on the Wilson-Cowan model and the Amari neural field model ^[21,22,23]. From a mathematical computational point of view, the heterogeneity of cortical connections in this model may play a role in cortical focusing. Following the work of Taylor et al., Tabatabaee et al. ^[16], and Yan et al. ^[17,18] considered for the first time the effect of excitatory coupling strength from EIN to PY during seizures. Fan et al. ^[24] found that feedforward excitation from TC to PY induced a transition from tonic oscillations to SWDs in the cortex. Cao et al. ^[25] altered the strength of the inhibitory coupling from IN to PY, and found that the system shifted intermittently from absence to tonic seizures. Based on the neurological dynamic diseases, many scholars performed bifurcation analyses of the transitions of different firing states to explore the nature of their generation processes ^[26].

Both physiological experiments and computational models show that the excitatory pyramidal neuronal population is involved in the induction of various firing patterns in the cortex ^{[27,28,29,30,31]}. On this basis, we improve the neural field model to investigate the excitatory effect of PY on seizures, which is directly controlled by the three nuclei, EIN, IN and TC, so it is particularly important to investigate the effect of these three pathways on seizures. This paper explores the dynamic transition phenomena caused by the excitatory projection of EIN on PY, the inhibitory projection of IN on PY, and the excitatory projection of TC on PY, respectively. The results show that all three pathways can induce seizure transitions. Section 2 provides a comprehensive explanation of the model and statistical indicators. The main findings are presented in Section 3, which includes a detailed analysis of the causes of the state transition using single-parameter and two-parameter bifurcation analysis. Finally, Section 4 discusses the implications of the results.

2.   Models and methods

The thalamocortical model consists of cortical subnetwork and thalamus subnetwork, as shown in Figure 1. The excitatory interneurons (EIN), the inhibitory interneurons (IN) and the excitatory pyramidal (PY) neuronal population compose the cortical subsystem, and the thalamic subsystem includes the specific relay nucleus (TC) and thalamic reticular nucleus (RE). Here, due to the dominant role of PY, we mainly consider the effects of EIN-PY, IN-PY and TC-PY pathways during seizures.

Figure 1.  Framework of the thalamocortical model. The cortical subnetwork contains the excitatory pyramidal neuronal population (PY), the excitatory interneuronal population (EIN) and the inhibitory interneuronal population (IN). The thalamic subnetwork contains the specific relay nucleus (TC) and thalamic reticular nucleus (RE). Red lines indicate excitatory projections, whereas green lines indicate inhibitory projections.

DownLoad: Full-Size Img PowerPoint

The following differential equations can be used to describe the improved model:

where $ \epsilon_1, ..., \epsilon_5 $ are additive constants, and $ \tau_1, ..., \tau_5 $ indicate different timescale parameters. $ C_{PY-PY}, ..., C_{RE-RE} $ represent the connectivity strengths between neuron populations, respectively. $ f[\cdot] $ and $ g[\cdot] $ are transfer functions, which are defined as follows:

where $ x = PY, IN, EIN, \text{and} \ TC $; $ y = TC\ \text{and}\ RE $; $ v $ is the steepness, $ v = 250, 000 $. Table 1 displays the parameter values for the model. All parameters are consistent with the literatures in Yan et al. ^[17], Fan et al. ^[24], Taylor et al. ^[32], and most of them are estimated from the original experimental data ^[33], while $ C_{EIN-PY}, C_{IN-PY} \ \text{and}\ C_{TC-PY} $ are considered as the primary bifurcation parameters.

Table 1.  Parameter values used in the model.

Parameter	Value	Parameter	Value	Parameter	Value	Parameter	Value
$ C_{PY-PY} $	1.8	$ C_{PY-TC} $	3	$ \tau_1 $	26	$ \epsilon_2 $	-3.4
$ C_{PY-EIN} $	0.1	$ C_{TC-PY} $	Varied	$ \tau_2 $	32.5	$ \epsilon_3 $	-0.1
$ C_{EIN-PY} $	Varied	$ C_{TC-RE} $	10.5	$ \tau_3 $	26	$ \epsilon_4 $	-2
$ C_{PY-IN} $	4	$ C_{RE-TC} $	0.6	$ \tau_4 $	2.6	$ \epsilon_5 $	-5
$ C_{IN-PY} $	Varied	$ C_{PY-RE} $	2	$ \tau_5 $	2.6	a	2.8
$ C_{IN-IN} $	0.05	$ C_{RE-RE} $	0.1	$ \epsilon_1 $	-0.5	b	0.5

---

 | Show Table

DownLoad: CSV

The numerical calculations are performed using MATLAB software, and the bifurcation diagrams are generated using AUTO, which is incorporated in XPPAUT. The differential equations of the system are solved using the standard fourth-order Runge-Kutta integration method with a fixed time step of 1 ms. To plot the extreme diagrams, the stable local maximum and minimum values of the time series of PY are calculated, and the spectrograms of the time series are estimated using fast Fourier transform (FFT). The dominant frequency is defined as the maximum peak frequency of the power spectral density (PSD).

3.   Results

The excitability of PY has been shown to play a major role in seizures. There are three neuronal nuclei in the model that have a direct effect on PY, which are EIN, IN and TC. PY is linked to EIN and TC via glutamatergic receptors and to IN via GABAergic receptors. Changes in the strength of synaptic connections are influenced by the release of neurotransmitters from excitatory and inhibitory neurons in the brain. Experimental studies have shown that in the rat cortex, changes in GABAergic and glutamatergic receptors function lead to absence, tonic and clonic seizures ^[34,35]. However, the transition dynamics between absence and tonic-clonic seizures caused by abnormal glutamatergic and GABAergic receptors in PY are unclear. We study the impact of the excitatory projection from EIN to PY, the inhibitory projection from IN to PY, and the excitatory projection from TC to PY on epileptic seizures in distinct sections. Here, we use $ C_{EIN-PY} $ to denote the strength of the excitatory connection from EIN to PY, $ C_{IN-PY} $ to denote the strength of the inhibitory connection from IN to PY, and $ C_{TC-PY} $ to denote the strength of the excitatory connection from TC to PY. Furthermore, the co-regulatory effects of $ C_{EIN-PY} $, $ C_{IN-PY} $ and $ C_{TC-PY} $ in seizures are also explored in detail. Single-parameter and two-parameter bifurcation analysis can uncover the dynamic mechanisms of state transitions.

The discharge patterns associated with epilepsy, as observed in EEG recordings, are classified based on the waveform and dominant frequency. These patterns include the seizure-free saturated state, 2-SWDs, SWDs, clonic, and tonic.

3.1.   Epilepsy activities modulated by the excitatory projection of EIN on PY

Tabatabaee et al. ^[16] first investigated the key role of EIN in seizures, but didn't consider the excitatory projection from EIN to PY, on the basis of which we consider the effect of changes in $ C_{EIN-PY} $ on epileptic seizures. As shown in Figure 3(a), with the increased excitatory action from EIN to PY, four different discharge patterns are produced: saturated state (I, Figure 2(a)), 2-SWDs state (II, Figure 2(b)), SWDs state (III, Figure 2(c)) and clonic state (IV, Figure 2(d)).

Figure 2.  Time series of PY: (a) Saturated state, ($ C_{EIN-PY} $, $ C_{IN-PY} $) = (0.0001, 1.5); (b) 2-SWDs, ($ C_{EIN-PY} $, $ C_{IN-PY} $) = (0.12, 1.5); (c) SWDs, ($ C_{EIN-PY} $, $ C_{IN-PY} $) = (0.3, 1.5); (d) Clonic, ($ C_{EIN-PY} $, $ C_{IN-PY} $) = (0.44, 1.5); (e) Tonic, ($ C_{EIN-PY} $, $ C_{IN-PY} $) = (0.8, 2.6). $ C_{TC-PY} $ is fixed at 1.

DownLoad: Full-Size Img PowerPoint

Figure 3.  (a) The extreme diagram induced by the projection from EIN to PY, $ C_{EIN-PY} $, red and blue dots represent the local maximum and minimum values, respectively. Regions are divided into four: (I) saturated state, (II) 2-SWDs, (III) SWDs and (IV) clonic. (b) Dominant frequency corresponding to (a), the frequency of the oscillatory state is in the range of 2–4 Hz, which is consistent with the frequency of absence and clonic seizures.

DownLoad: Full-Size Img PowerPoint

The total transition behaviors mediated by $ C_{EIN-PY} $ are shown in Figure 3(a), where red and blue dots indicate local maximum and minimum values, respectively. It is observed that the system exhibits saturated discharge state (I) when $ C_{EIN-PY} $ is sufficiently small. Increasing $ C_{EIN-PY} $, the system progresses to 2-SWDs (II), SWDs (III), and clonic (IV). When $ C_{EIN-PY} $ is significantly large, the local maximum and minimum values are nearly identical, causing the system to revert back to the saturated state (I). In addition, as shown in Figure 3(b), the oscillation frequency ranges from 2 to 4 Hz, which is consistent with the dominant frequency of clonic oscillations and absence seizures.

The single-parameter bifurcation of the system concerning $ C_{EIN-PY} $ is shown in Figure 4. The stable and unstable fixed points are represented by the red solid and black dashed lines, respectively, while the stable and unstable limit cycles are depicted by the green and dark blue lines, respectively. In detail, for smaller $ C_{EIN-PY} $ (i.e., $ C_{EIN-PY} $ $ < \approx $ 0.07543), the system reaches a stable fixed point, which corresponds to the saturated state shown in Figure 2(a). As $ C_{EIN-PY} $ increases, the fold limit cycle bifurcation ($ \rm LPC_1 $) gives rise to a bistable region comprising a stable limit cycle and a stable fixed point (i.e., 0.07543 $ < \approx $ $ C_{EIN-PY} $ $ < \approx $ 0.14929). Further increase in $ C_{EIN-PY} $ leads to the reoccurrence of the fold limit cycle bifurcation ($ \rm LPC_2 $), resulting in the appearance of another stable and unstable limit cycles. This transforms the system into a tristable region, which consists of one stable fixed point and two stable limit cycles (i.e., 0.14929 $ < \approx $ $ C_{EIN-PY} $ $ < \approx $ 0.15875). Due to the appearance of $ \rm LPC_3 $ at $ C_{EIN-PY} $ $ \approx $ 0.15875, the stable limit cycle generated by $ \rm LPC_1 $ disappears, causing the system to enter a bistable region composed of a stable fixed point and a stable limit cycle generated by $ \rm LPC_2 $, within the range of 0.15875 $ < \approx $ $ C_{EIN-PY} $ $ < \approx $ 0.20743. Subsequently, at the point $ C_{EIN-PY} $ $ \approx $ 0.20743, the appearance of the Hopf bifurcation ($ \rm HB_1 $) causes the stable equilibrium point to become unstable, and the system enters a single stable state region composed of stable limit cycles. Until another Hopf bifurcation ($ \rm HB_2 $) occurs, creating a new unstable limit cycle, the unstable fixed point turns into a stable fixed point, and the system restores to the bistable state (i.e., 0.4008 $ < \approx $ $ C_{EIN-PY} $ $ < \approx $ 0.44182). At last, the fold limit cycle bifurcation ($ \rm LPC_4 $) occurs at $ C_{EIN-PY} $ $ \approx $ 0.44182, leading to the disappearance of the stable and unstable limit cycles and the system enters a monostable state consisting of stable fixed points.

Figure 4.  Dynamical bifurcation diagram corresponding to Figure 3(a). Numbers 1–3 respectively represent the monostable, bistable, and tristable regions. Red solid and black dashed lines represent stable and unstable fixed points; green and dark blue lines are stable and unstable limit cycles. The pink and light blue marking points indicate fold limit cycle bifurcations ($ LPC_{i}, i = 1, 2, 3, 4 $) and Hopf bifurcations ($ HB_{j}, j = 1, 2 $), respectively. Parameters $ C_{IN-PY} $ = 1.5, $ C_{TC-PY} $ = 1.

DownLoad: Full-Size Img PowerPoint

3.2.   Epilepsy activities modulated by the inhibitory projection of IN on PY

In this subsection, we mainly explore the effect of the inhibitory projection from IN to PY on seizures. As shown in Figure 5(a), the system goes through saturated state, SWDs, 2-SWDs, saturated state, and tonic.

Figure 5.  (a) The extreme diagram induced by the projection from IN to PY, $ C_{IN-PY} $, red and blue dots represent the local maximum and minimum values, respectively. Regions are divided into four: (I) saturated state, (III) SWDs, (II) 2-SWDs and (V) tonic. (b) Dominant frequency corresponding to (a), the frequency of absence seizures is in the range of 2–4 Hz and the frequency of tonic seizures is above 25 Hz.

DownLoad: Full-Size Img PowerPoint

The single-parameter bifurcation of the system with respect to $ C_{IN-PY} $ is shown in Figure 6. As $ C_{IN-PY} $ varies in the range of [1, 1.67871], showing a stable fixed point of the system corresponding to the saturated state. Afterward, as $ C_{IN-PY} $ continuously varies from 1.67871 to 1.69792, the fold limit cycle bifurcation ($ \rm LPC_1 $) occurs, resulting in the appearance of one stable and one unstable limit cycle. Hence, this region becomes bistable, with a stable fixed point and a stable limit cycle. Next, the Hopf bifurcation ($ \rm HB_1 $) appears at $ C_{IN-PY} $ $ \approx $ 1.69792, where the unstable limit cycle disappears and the fixed point loses its stability, and the region [1.69792, 1.78611] becomes monostable. Then, the system enters into a bistable region when 1.78611 $ \leq $ $ C_{IN-PY} $ $ \leq $ 1.83806. The reason why the unstable fixed point changes into a stable fixed point and a new unstable limit cycle creates is the appearance of Hopf bifurcation ($ \rm HB_2 $), at $ C_{IN-PY} $ $ \approx $ 1.78611. Subsequently, because the fold limit cycle bifurcation ($ \rm LPC_2 $) arises at $ C_{IN-PY} $ $ \approx $ 1.83806, the stable and unstable limit cycle disappears, and the system turns into a monostable region [1.83806, 2.35184]. Finally, the system remains in a monostable region [2.35184, 3], and $ \rm HB_3 $ appears at $ C_{IN-PY} $ $ \approx $ 2.35184, which changes the fixed point to unstable and generates a new stable limit cycle.

Figure 6.  Dynamical bifurcation diagram corresponding to Figure 5(a). Numbers 1 and 2 respectively represent the monostable and bistable regions. Red solid and black dashed lines represent stable and unstable fixed points; green and dark blue lines are stable and unstable limit cycles. The pink and light blue marking points indicate fold limit cycle bifurcations ($ LPC_{i}, i = 1, 2 $) and Hopf bifurcations ($ HB_{j}, j = 1, 2, 3 $), respectively. Parameters $ C_{EIN-PY} $ = 0.8, $ C_{TC-PY} $ = 1.

DownLoad: Full-Size Img PowerPoint

3.3.   Epilepsy activities modulated by the excitatory projection of TC on PY

The excitatory effect of TC on PY in the regulation of seizure discharges is shown in Figure 7(a). The system transits from saturated state to clonic oscillation and then returns to saturated resting state. Furthermore, it can be inferred from Figure 7(b) that the clonic oscillation state falls within the frequency range of 2–4 Hz.

Figure 7.  (a) The extreme diagram induced by the projection from TC to PY, $ C_{TC-PY} $, red and blue dots represent the local maximum and minimum values, respectively. Regions are divided into two: (I) saturated state and (IV) clonic. (b) Dominant frequency corresponding to (a), the frequency of clonic seizures is in the range of 2–4 Hz.

DownLoad: Full-Size Img PowerPoint

The single-parameter bifurcation of the system concerning $ C_{TC-PY} $ is shown in Figure 8. Initially, as $ C_{TC-PY} $ increases from 0 to 0.3028, the system exhibits saturated state characterized by stable fixed points, resulting in a monostable state. However, as $ C_{TC-PY} $ continues to increase, the system reaches the Hopf bifurcation point ($ \rm HB_1 $), causing the fixed point to lose stability and a new stable limit cycle to emerge. The system then remains in a monostable state with a stable limit cycle. Eventually, at the occurrence of another Hopf bifurcation point ($ \rm HB_2 $), the stable limit cycle disappears, and the fixed point becomes stable again, leading to the system remaining in a monostable state.

Figure 8.  Dynamical bifurcation diagram corresponding to Figure 7(a). The Number 1 represents the monostable region. Red solid and black dashed lines represent stable and unstable fixed points; green line is stable limit cycle. The light blue marking points indicate Hopf bifurcations ($ HB_{j}, j = 1, 2 $). Parameters $ C_{EIN-PY} $ = 0.8, $ C_{IN-PY} $ = 1.5.

DownLoad: Full-Size Img PowerPoint

3.4.   Two-parameter bifurcation analysis

In the previous subsection, we separately analyzed the effects of the three pathways EIN-PY, IN-PY, and TC-PY on seizures using a single-parameter bifurcation. To further analyze the co-regulatory effects of the above three pathways on seizures, we perform two-parameter bifurcation analysis. Firstly, we investigate the state transitions and corresponding dynamical mechanisms when the parameters $ C_{EIN-PY} $ and $ C_{TC-PY} $ are varied in the range of [0, 0.8] $ \times $ [0, 1]. Figure 9(a) provides a detailed view of the state regions under the panel ($ C_{EIN-PY} $, $ C_{TC-PY} $), which are classified into the saturated region (I) and oscillation regions (II–V). When $ C_{EIN-PY} $ and $ C_{TC-PY} $ are small, tonic occurs in the region (V). From a vertical perspective, when $ C_{EIN-PY} $ is small, an increase in $ C_{TC-PY} $ only leads to a transition from tonic to saturated state. When $ C_{EIN-PY} $ increases slightly, the system develops 2-SWDs. The area of SWDs is larger when $ C_{EIN-PY} $ lies near 0.4. As $ C_{EIN-PY} $ increases further, the area of SWDs gradually shrinks. These observations indicate that the modulation of epileptic seizures is influenced by both $ C_{EIN-PY} $ and $ C_{TC-PY} $. Figure 9(b) displays the bifurcation of the system for $ C_{EIN-PY} $ and $ C_{TC-PY} $, where the Hopf bifurcation curve and the fold limit cycle bifurcation curve are represented by the blue and black lines, respectively. These curves illustrate how the occurrence of SWDs and state transitions induced by $ C_{EIN-PY} $ and $ C_{TC-PY} $. When one of the parameters is fixed, it is easy to obtain the transition mechanism for the other parameter and to determine the exact type and position of the bifurcation. Due to the complexity of the system solution, the panel is divided into numerous multistable regions. The two narrow regions between the blue and black curves, i.e., the regions between curves 2 and 3 and between curves 4 and 5 are known as the bistable regions of the system. These regions allow for the coexistence of stable and unstable oscillations. In addition, between the two black curves, i.e., between curves 2 and 5, SWDs and clonic oscillations can occur, corresponding to regions (II–IV) in Figure 9(a). There is also the region enclosed by the Hopf bifurcation curve (i.e., curve 1) and the coordinate axis where epileptic tonic oscillations can occur, corresponding to the region (V) in Figure 9(a). In other areas, stable fixed points indicating saturated firing of PY can appear (shown in region (I) in Figure 9(a)).

Figure 9.  (a) State evolution diagram under the combined influence of connectivity strengths $ C_{EIN-PY} $ and $ C_{TC-PY} $. Seizure-free (I), 2-SWDs (II), SWDs (III), Clonic (IV), Tonic (V) are indicated by text labels; (b) Two-parameter bifurcation diagram by varying ($ C_{EIN-PY} $, $ C_{TC-PY} $) in regimes [0, 0.8] $ \times $ [0, 1]. The blue and black lines are the Hopf bifurcation curve and the fold limit cycle bifurcation curve, respectively. MS represents monostable region, BS represents bistable region, TS represents tristable region. Parameter $ C_{IN-PY} $ = 1.5.

DownLoad: Full-Size Img PowerPoint

The synaptic connection strengths $ C_{EIN-PY} $ and $ C_{IN-PY} $ jointly regulate epileptic discharges, and the state transitions and corresponding dynamical mechanisms are investigated when $ C_{EIN-PY} $ and $ C_{IN-PY} $ are varied in the range of [0, 0.8] $ \times $ [1, 3]. The panel ($ C_{EIN-PY} $, $ C_{IN-PY} $) in Figure 10(a) shows a detailed view of state regions, which are categorized into four different regions: the background saturated region (I), tonic region (V), 2-SWDs region (II), and SWDs region (III).

Figure 10.  (a) State evolution diagram under the combined influence of connectivity strengths $ C_{EIN-PY} $ and $ C_{IN-PY} $. Seizure-free (I), 2-SWDs (II), SWDs (III), Tonic (V) are indicated by text labels; (b) Two-parameter bifurcation diagram by varying ($ C_{EIN-PY} $, $ C_{IN-PY} $) in regimes [0, 0.8] $ \times $ [1, 3]. The blue and black lines are the Hopf bifurcation curve and the fold limit cycle bifurcation curve, respectively. MS represents monostable region, BS represents bistable region. Parameter $ C_{TC-PY} $ = 1.

DownLoad: Full-Size Img PowerPoint

Specifically, from a horizontal perspective, when the value of $ C_{IN-PY} $ is small or when the value of $ C_{IN-PY} $ is large (above 2.4), the system will only produce saturated discharge state or tonic state regardless of any further variation in $ C_{EIN-PY} $. From a vertical perspective, by fixing an arbitrary value of $ C_{EIN-PY} $ so that $ C_{IN-PY} $ increases from 1 to 3, the system goes through saturated state, SWDs, 2-SWDs, saturated state, and tonic. The two-parameter bifurcation of the system for $ C_{EIN-PY} $ and $ C_{IN-PY} $ is shown in Figure 10(b). The regions between curves 2 and 4 and between curves 5 and 6 are the bistable regions of the system. 2-SWDs and SWDs can occur between the fold limit cycle bifurcation curves, i.e., between curves 2 and 6, corresponding to regions (II and III) of Figure 10(a). The region enclosed by the Hopf bifurcation curve (i.e., curve 1) and the coordinate axis, where tonic oscillations can occur, corresponds to the region (V) of Figure 10(a).

The synaptic connection strengths $ C_{IN-PY} $ and $ C_{TC-PY} $ jointly regulate epileptic discharges. The state transitions and bifurcation mechanisms are investigated when $ C_{IN-PY} $ and $ C_{TC-PY} $ are varied in the range of [1, 3] $ \times $ [0, 1]. The state transitions of the system under the combined influence of parameters $ C_{IN-PY} $ and $ C_{TC-PY} $ are shown in Figure 11(a). When the values of $ C_{IN-PY} $ and $ C_{TC-PY} $ are small, the system will only produce saturated discharge state. As $ C_{IN-PY} $ and $ C_{TC-PY} $ increase, the system undergoes a transition from the background saturated state to the SWDs. When the value of $ C_{IN-PY} $ is large, the system will only produce tonic, regardless of the value of $ C_{TC-PY} $. The two-parameter bifurcation of the system for $ C_{IN-PY} $ and $ C_{TC-PY} $ is shown in Figure 11(b), the regions between curves 1 and 2 and between curves 3 and 4 are the bistable regions of the system. The region between the two fold limit cycle bifurcation curves (i.e., the two black curves 1 and 4), where 2-SWDs and SWDs occur, corresponds to regions (II and III) of Figure 11(a). The region enclosed by the Hopf bifurcation curve (i.e., curve 5) and the coordinate axes, undergoes tonic oscillations, corresponding to the region (V) of Figure 11(a).

Figure 11.  (a) State evolution diagram under the combined influence of connectivity strengths $ C_{IN-PY} $ and $ C_{TC-PY} $. Seizure-free (I), 2-SWDs (II), SWDs (III), Tonic (V) are indicated by text labels; (b) Two-parameter bifurcation diagram by varying ($ C_{IN-PY} $, $ C_{TC-PY} $) in regimes [1, 3] $ \times $ [0, 1]. The blue and black lines are the Hopf bifurcation curve and the fold limit cycle bifurcation curve, respectively. MS represents monostable region, BS represents bistable region. Parameter $ C_{EIN-PY} $ = 0.8.

DownLoad: Full-Size Img PowerPoint

4.   Conclusions

The excitatory pyramidal neuronal population plays a key role in seizures. In this paper, we simplify the model based on Yan et al. ^[17,18] and exploit the thalamocortical neural field model to investigate the effects of excitability of EIN-PY and TC-PY and inhibition of IN-PY on epileptic seizures.

Firstly, by changing the excitatory effect of EIN on PY, the system can go through saturated state, 2-SWDs, SWDs, clonic, and saturated state in sequence. In terms of dynamics, the emergence of fold limit cycle bifurcation and Hopf bifurcation is due to the increase in connection strength, which places the system in the multistable region. Secondly, when the inhibitory effect of IN on PY is changed, the system generates a transition from saturated discharge state to tonic oscillation state in addition to producing SWDs and 2-SWDs. This phenomenon arises because the Hopf bifurcation that changes the fixed point from stable to unstable and creates a new stable limit cycle. Then, when the excitatory effect of TC on PY is changed, the Hopf bifurcation occurs, which puts the system to produce a state of clonic oscillation and eventually transition to saturation. Finally, using two-parameter bifurcation analysis, we find that a wider variety of seizure states could be induced by the combined effect of EIN and TC on PY. At the same time, the stable parameters of tonic-clonic oscillations, SWDs, 2-SWDs, and saturated discharges can be obtained in the two-dimensional plane, and the reasonable transition path between tonic-clonic and absence seizures can be identified.

The presence of multiple stable regions implies that external disturbances could trigger shifts between the fixed point and limit cycle. Therefore, appropriate electrical, electromagnetic, and optogenetic stimulation may suppress seizures. This also provides new directions for future research.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant Nos. 12002001 and 12272002) and Cultivation Plan for "Yujie'' Team of North China University of Technology (No. 107051360022XN725). The authors would like to thank all the reviewers for their valuable comments.

Conflict of interest

The authors declare there is no conflicts of interest.