Dynamics of a coupled epileptic network with time delay

1.   Introduction

Epilepsy is a chronic brain disease, which is characterized by the occurrence of spontaneous seizures ^[1]. The affected neurons produce synchronized abnormal discharges during seizure ^[2]. There are about 50 million patients with epilepsy globally, and there are two million new patients each year. The focal onset seizure is the most common type for patient with epilepsy. The onset starts locally and propagates to normal brain tissues ^[3,4]. There are interactions between pathological and normal brain regions ^[5]. The epilepsy is now considered as a brain network disease. The epileptic discharges are with complex temporal and spatial characteristics ^[6,7]. How the seizure onselectroencephalogram (EEG)ets and propagates to neighbor regions is still in an infancy ^[8]. Some early studies proposed mathematical models to simulate the seizure onset, propagation and termination, which are critical to drug treatment, surgical intervention and neuromodulation of epilepsy ^[9,10].

There are different approaches to model the epileptic behaviors of the brain. Wendling $et\ al$. introduced a computational macroscopic model of electroencephalogram (EEG) activity that included a physiologically relevant fast inhibitory feedback loop and the transition from interictal to fast ictal activity was explained by the impairment of dendritic inhibition in the model ^[11,12]. Kalitzin $et\ al$. proposed an analytical model as examples of transitions between various dynamical states ^[13]. The model was then used to study the dynamics of convulsive seizure termination and postictal generalized EEG suppression ^[14]. Jirsa $et\ al$. proposed a generic model called Epileptor. They showed that the onset and offset of ictal-like discharges were well-defined mathematical events: a saddle-node and homoclinic bifurcation, respectively ^[15].

Furthermore, the coupled Epileptor network can be used to predict the spatiotemporal diversity of seizure propagation and termination in human focal epilepsy ^[16]. Time lags between epilepsy network nodes were considered as an important fact linked to seizure onset ^[17]. Bandt $et\ al$. studied connectivity strength and time lags between different network nodes using resting-state functional Magnetic Resonance Imaging (fMRI). They showed that there were decreased time lags within the seizure onset node and globally increased time lags throughout all regions of the brain not involved in seizure onset or propagation ^[18].

In this study (see Figure 1), we investigated that the dynamical behaviors of epilepsy network based on simplified Epileptor model with time delay, and how the variation of time delay between the network nodes induce seizure-like event.

2.   Local stability and Hopf bifurcation

The Epileptor model was introduced to represent complex behavior of the brain for patient with epilepsy. Taking advantage of time scale separation and focusing on the slower time scale, two dimensional Epileptor was proposed as follows ^[16]:

where $i\in\left\{1, 2\right\}$; characteristic time scale is fixed at $\tau_0 = 2857$, and resting-state current $I_{1, i} = 3.1.$ $x_{1, i}$ is the state variables on the slower time scale account for spike and wave events observed in electrographic seizure recordings. $z$ is the state variable which guides the neural population through the seizures including seizure onset and offset. $K_{ij}$ is the connection strength between Epileptor $i$ and Epileptor $j$ as given by the connectivity matrix with $i, j = 1, 2$. The parameter $x_{0, i}$ describes the excitability degree of each Epileptor.

Some studies have shown that the time lags between different brain regions will change in patients with epilepsy. Proix $et \ al$ found delays of recruitment can be either highly variable or stereotypical by computing for several seizures the mean and SD of the delays of recruitment and the values of excitability and coupling parameters ^[8]. Zhang $et\ al$ found that time delay is related to coupled seizure network, which is unavoidably encountered due to finite speed, reaction time of signal transmission over the couplings, can lead to instabilities, bursting transitions, bifurcations of periodic ^[19]. In this paper, we consider the following epileptic model with two time delays between coupled nodes:

The nontrivial equilibrium point of model system (2.1) is $P(x_{11}^*, z_1^*, x_{12}^ *, z_2^*),$ where

and $x_{11}^*$ satisfies the following equation:

where

The characteristic equation of model (2.1) at the nontrivial equilibrium point $P$ is given by:

which can be reduced into the following form

where

When $\tau_1 = \tau_2 = 0$, the characteristic equation (2.2) becomes

According to Routh-Hurwitz criteria, the positive equilibrium of model (2.1) in the absence of time delay is locally asymptotically stable when the following conditions are satisfied:

Next, we will analyze the effects of the two time-delays on the stability of the nontrival equilibrium point $P$. We substitute $i\omega$ into equation(2.2) and separate the real and imaginary parts, then we obtain the following equations:

from which it follows that

If $A^2_4-A^2_5 < 0,$ the Equation (2.5) has a positive real root $\omega_{0}^*$. We define $\tau = \tau_{1}+\tau_{2}$, then we get

According to the Butler's Lemma, we can conclude that the equilibrium point $P$ is stable for $\tau < \tau^*_{0}.$ Now, we discuss whether Hopf bifurcation occurs at $P$ of model (2.1) when $\tau$ increases through $\tau^*_{0}$.

Based on Equation (2.2), we can further obtain

From Equation (2.6), it is clear that the sign can be determined by

where $f(x) = x^4+(A^2_1-2A_2)x^3+(2A_4+A^2_2-2A_1A_3)x^2+(A^2_3-2A_2A_4)x+A^2_4-A^2_5$.

Theorem 1 Assume that the conditions in (2.3), $A^2_4-A^2_5 < 0$ and $f^{'}(\omega^{*2}_{0}) \neq 0$ hold. The equilibrium point $P$ of model (2.1) is locally asymptotically stable for $\tau \in [0, \tau^*_{0})$, and Hopf bifurcation occurs at $P$ when $\tau = \tau^*_{0} \ $, where

3.   Stability of bifurcated periodic solutions

Generally speaking, the time delays between two network nodes are always similar, so we consider a special case $\tau_1 = \tau_2$. In what follows, we will investigate the direction of Hopf bifurcation, stability and period of the periodic solution bifurcating from the equilibrium point $P$. Following the ideas of Hassard et al. ^[20], we derive the explicit formulae for determining these properties of Hopf bifurcation at the critical value $\tau^*_{s} = \frac{\tau^*_{0}}{2}$ by employing the normal form method and center manifold theorem. Let

and $u_1(t) = x_{11}(\tau_st), u_2(t) = z_1(\tau_st), u_3(t) = x_{12}(\tau_st), u_4(t) = z_2(\tau_st), \tau_s = \tau^*_{s}+\mu, \mu\in R.$ For system(2.1), $\mu = 0$ is the Hopf bifurcation value. In the fixed phase space $C = C([-1, 0], R^4)$, system (2.1) is transformed into a FDE as

where $u(t) = (u_1(t), u_2(t), u_3(t), u_4(t))^T\in R^4$, and $L_\mu :$ $C \rightarrow R, F:R\times C\rightarrow R$ are given respectively by

and

where $\phi(\theta) = (\phi_1 (\theta), \phi_2 (\theta), \phi_3 (\theta), \phi_4 (\theta))^T\in C.$

According to the Riesz representation theorem, there exists a $4 \times 4$ matrix $\eta(\theta, \mu)$ of bounded variation for $\theta \in [-1, 0]$, which follows that

In fact, we can choose

For $\phi \in C^1([-1, 0], R^4)$, we define

and

Then system (3.1) is equivalent to

where $u_t(\theta) = u(t+\theta)$, for $\theta\in[-1, 0].$

For $\psi\in C^1([0, 1], (R^4)^*),$ we define

and the bilinear form

where $\eta(\theta) = \eta(\theta, 0).$ Thus, $A = A(0)$ and $A^*$ are adjoint operators.

Owning to $\pm i\omega^*_{0}\tau^*_{s}$ are eigenvalues both of $A$ and $A^*$, we can compute the eigenvector of $A$ corresponding to $i\omega^*_{0}\tau^*_{s}$ and $A^*$ corresponding to $-i\omega^*_{0}\tau^*_{s}.$ Assume that the eigenvector of $A(0)$ corresponding to $i\omega^*_{0}\tau^*_{s}$ is $q(\theta) = (1, \Delta_1, \Delta_2, \Delta_3)^T e^{i\omega^*_{0}\tau^*_{s}\theta}$, then $A q(\theta) = i\omega^*_{0}\tau^*_{s}q(\theta).$ The definition of $A(0)$ and $\eta(\theta, \mu)$ yields

Thus, we can get

Similarly, the eigenvector of $A^*$ corresponding to $-i\omega^*_{0}\tau^*_{s}$ is $q^*(s) = D(1, \Delta^*_1, \Delta^*_2, \Delta^*_3) e^{i\omega^*_{0}\tau^*_{s}s}$, then we get

From (3.4), we get

Thus, we get

Therefore, it is clear that both $\langle q^*(s), q(\theta) \rangle = 1$ and $\langle q^*(s), \overline {q}(\theta) \rangle = 0$ are satisfied. Let $u_t$ be the solution of (3.3) when $\mu = 0.$ Define

On the center manifold $C_0$, we get

and

where $z$ and $\overline{z}$ are local coordinates for center manifold $C_0$ in the direction of $q^*$ and $\overline{q^*}.$ Note that $W$ is real if $u_t$ is real. We only consider real solutions. For solution $u_t \in C_0$ of (3.3), due to $\mu = 0$, we have

We rewrite this equation as

where

From (3.5) and (3.6), we have

where $q(\theta) = (1, \Delta_1, \Delta_2, \Delta_3)^T e^{i\omega^*_{0}\tau^*_{s}\theta}$, then

Comparing the coefficients with (3.7), we can yield the following important parameters:

Since $g_{21}$ depends on the coefficient $W_{20}(\theta)$ and $W_{11}(\theta),$ we need to compute them. From (3.3) and (3.5), we can get

On the other hand, on the center manifold $C_0$ near the origin, we can obtain

By comparing the coefficients of $\frac{z^2}{2}$ and $z \overline{z}$, we know that

Similarly,

We know that for $\theta \in [-1, 0),$

Substituting $q(\theta) = q(0)e^{iw^*_{0}\tau^*_s\theta}$ into(3.10), then we have

where $E_1 = (E^{(1)}_1, E^{(2)}_1, E^{(3)}_1, E^{(4)}_1)\in R^4$ is a constant vector.

Similarly, we obtain

where $E_2 = (E^{(1)}_2, E^{(2)}_2, E^{(3)}_2, E^{(4)}_2)\in R^4$ is also a constant vector.

Based on (3.10), (3.11), (3.12) and (3.13), we can obtain

where

and

where

Thus, $W_{20}(\theta)$ and $W_{11}(\theta)$ can be determined and $g_{20}$, $g_{11}$, $g_{02}$, $g_{21}$ can be expressed. Hence, we can judge property of Hopf bifurcation by three parameters $\mu_2$, $\beta_2$, $T_2$ and compute the following values:

From the conclusion of Hassard et al. ^[20], we have the conclusion

Theorem 2 $\mu_2$ determines the direction of the Hopf bifurcation, $\beta_2\ $determines the stability of the bifurcating periodic solution, and $T_2$ determines the period of the bifurcating periodic solution. Moreover, if $\mu_2 > 0$ $(\mu_2 < 0),$ the Hopf bifurcation is supercritical (subcritical); if $\beta_2 > 0$ $(\beta_2 < 0),$ the bifurcating periodic solution is stable (unstable); if $T_2 > 0$ ($T_2 < 0$), the period increases (decreases).

4.   Numerical simulations

We use some numerical simulation to illustrate the bifurcation analysis. Firstly, we consider the following Epileptor model:

where $x_{01}$ and $I_{11}$ are regarded as the two parameters. We plot a critical curve with respect to these two parameters in Figure 2, which gives the stable region. Moreover, fixing

we consider a coupled Epileptor model with the following form

The pathological changing of brain tissues always increase or decrease the time lags of signal transmission bidirectionally at the same time. Therefore, we consider a special scenario $\tau_1 = \tau_2$. We can compute the nontrivial equilibrium point $P\quad (0.0114, 4.0997, 0.3640, 3.7870)$ of model (4.1). In the absence of delay, $\tau_1 = \tau_2 = 0$, the dynamic behaviors of model (4.1) are described in Figure 3, which indicates that $P$ is locally asymptotically stable. Using Matlab software, we can find the unique positive solution $\omega^*_{0} = 0.3$ in (2.5). Then, we can obtain $\tau_s^* = 2.6$, $c_1(0) = -3.69-1.36i$, $\mu_2 = 8.62$, $\beta_2 = -0.85$, $T_2 = -5.41$. Based on the Theorem 2, model (4.1) undergoes a supercritical Hopf bifurcation at the nontrivial equilibrium point $P$ and the bifurcating periodic solution exists for $\tau_s$ slightly larger than $\tau_s^*$ and the bifurcated periodic solution is unstable, which can be seen in Figure 4. We also plot Figure 5 to present phase map between $x_{11}$ and its delay $x_{11}^d$.Moreover, we take $\tau_0$ as the parameter and exhibit a critical curve in the plan $\tau_0-\tau_s$. From Figure 6, we can see that time delay decreases with the increasing of characteristic time scale $\tau_0$.We also consider other cases such as $\tau_2 = 0.5\tau_1$, $\tau_2 = 2\tau_1$ and $\tau_2 = 4\tau_1$. The periodic solution still exist which can be seen in Figure 7.

5.   Discussions

Our theoretical analysis and numerical simulation both show that the increase of time delay between nodes change the network to SLE state through Hopf bifurcation. The bifurcation value $\tau_s$ also relies on characteristic time scale of brain tissue $\tau_0$. The increasing of $\tau_0$ will increase the threshold of bifurcation value $\tau_s$. Time delay between nodes of epileptic network is important for understanding the seizure. The seizure onsets are mainly due to the changes of brain tissues, which alter the topology of epileptic network. Early studies using diffusion tensor imaging (DTI) indicated the decrease in fractional anisotropy (FA) of brain regions of patients with epilepsy, which further changed the time delay of signal transmission between brain regions ^[21]. Study using fMRI showed global increased time lags through all regions of the brain not involved in seizure onset and propagation ^[18]. Our study show that the stable network may change to SLE state by introducing time delays between nodes, which indicates that time delay is one of the important reasons for seizure onset. Moreover, our analysis shows that the epileptic network of two coupled nodes goes to the SLE state through Hopf bifurcation.

Our results support that the brain region with increased time delay with other regions can be selected for epilepsy treatment. Noninvasive neuromodulation is a potential measure to treat drug-resistant epilepsy, which use electric or magnetic field to stimulate specific brain region in order to eliminate seizure, such as repetitive Transcranial Magnetic Stimulation (rTMS), transcranial Alternating Current Stimulation (tACS), etc. Where and when to make the stimulation are the key factors for efficacy of treatment. Our theoretical study shows that the changing in time delay between brain regions is a good indicator for forecasting seizure onset. Scalp EEG is a convenient and noninvasive tool of measuring electrophysiological activities of brain. There are also some measures, such as phase lag index, to quantify time lags of EEG signal transmission between different recording sites ^[22]. Time lags of nodes of epileptic network can be used to select target for epilepsy treatment.

Conflict of interest

The authors declared that there is no conflicts of interest in this paper.