The preview control of a corticothalamic model with disturbance

1.   Introduction

Epilepsy is a neurological disorder characterised by abnormal brain rhythms recorded at the macroscopic level of EEG electrodes typically during seizures ^[1]. The electrographic epileptiform phenomena in the case of absence seizures, myoclonic seizures and complex partial seizures are often comprised of periodic spike and wave discharges (SWD) in the EEG, which is chacterized by a fast spike followed by a slow wave ^[2].

Epileptic seizure is fundamentally a dynamic disease ^{[3, 4]} which occurs in an intact physiological system. This system operates in a range of control parameters that lead to abnormal dynamics ^{[5, 6, 7]}. For the detection of seizures, scholars are also exploring new markers ^[8]. The further clinical challenge is to take into account the evolving dynamics of epilepsy with the goal of developing effective treatment strategies. There is a considerable literature on the possible mechanism of generation of SWD which is ascribed to the most commonly accepted 'cortico-reticular' theory ^[9], i.e., SWD bursts are generated by an interplay between thalamus and cortex. Thus, the computational challenge is to understand the relationship between the structure of the nervous system and the dynamic abnormalities generated by epileptic neural populations ^{[10, 11]}.

Investigations showed that a brief sensory or electrical stimulus given at seizure onset can sometimes incur ^{[12, 13]}, abate ^[14] or abort ^{[15, 16, 17]} the seizure events of epilepsy. The mantra of dynamic systems theory is that qualitative changes in dynamics, referred to as bifurcations, occur as an important parameter(s) crosses stability boundaries ^{[18, 19]}. In fact, epileptic phenomena are generally modeled by a bistability of silence and seizure-like bursting. It is demonstrated in a bistable model that single pulse perturbations can both induce and abate abnormal epileptiform activity ^{[20, 21]}. This approach is to apply single pulse perturbations in state space beyond the manifold which separates the seizure and non-seizure attractor ^{[22, 23]}. One of the assumptions of our study is that the background state coexists with the SWD limit cycle in the state space ^[23].

While there is obvious appeal to single pulse stimulation, there are many difficulties with that approach, especially in stochastic systems where repeated success can be troublesome ^[22]. In addition, stimulation involves many control parameters such as stimulation cycle, pulse width, intensity and stimulation direction, which are difficult both theoretically and clinically. In contrast, nonlinear feedback controller and observer are designed to suppress the epileptic seizures ^{[24, 25]}. Neural Mass Model (NMM) is a biologically inspired model able to reproduce signals observed in spontaneous electroencephalograms (and evoked potentials), and it allows one to apply those controllers.

The control of a system with epileptic SWD oscillations is highly nontrivial since the system is nonlinear ^[26]. Broadly speaking, optimal control, according to the control theory, is a mathematical framework that allows for the systematic selection of time-varying inputs to effectively drive a dynamical system in a desired way ^{[22, 27]}. Optimal control theory has been utilized not just in seizures ^[27], but also in the studies of other disease models, including the SIR Model of infectious illnesses ^{[28, 29]} and the complication Type $2$ diabetes ^[30]. However, for epileptic seizures, the nonlinear controller designed in the previous literature can effectively control it ^{[20, 27]}, they do not consider how to effectively eliminate or suppress random disturbances. From this perspective, we try to utilize the theory of preview control to attenuate or eliminate the effects of random disturbances on the epileptic system. Preview control can make full use of the known future information to reduce the static error and to improve the performance of inhibiting transient responses ^{[31, 32, 33, 34]}. Optimizing the utilization of the known future reference signals or disturbance signals to improve the control performance of epileptic seizure system continues to be a significant area of research.

Indeed, in addition to impulse stimulation, the random disturbance (real brain environmental noise or other external stimulation interference) has been theoretically demonstrated to induce epileptic seizures and further propagation due to the mechanisms of bistability or excitability ^[4]. That is, when the controlled system receives the disturbance again, a transient oscillation may be induced, which as the disturbance input may also induce the seizure of the other adjacent coupled nodes. This might be an issue that almost all controller designs ignore.

We conduct the first study attempting to use preview controllers to control epileptic seizures. We explore whether the preview controller can stabilize the controlled system (the current system) after a second disturbance. Therefore, our contributions and innovations are as follows:

$\bullet$ In this study, based on a computational model of epileptic spike-wave dynamics with disturbances, we aim to design a preview controller that is capable of suppressing epileptic seizures after receiving impulsive and noisy disturbances. We compare the control effect of the preview controllers and non-preview controllers, and the results show that the preview controllers have better effect.

$\bullet$ In particular, we consider single-input and multi-input cooperative control strategies to explore the feasibility and effectiveness of clinical multi-target stimulation in epilepsy regulation from a computational perspective. By discussing the control strategies of various nuclei, we aim to explore the strategies that can make the system oscillate less or even not oscillate after secondary disturbance. This research serves as a point of reference for clinical epilepsy regulation.

$\bullet$ In addition, we propose a cost performance function for evaluating the energy consumption of preview control and non-preview control under different input control strategies. The lower performance consumption can extend the operating time of controllers. Comparing the consumption of different strategies, we discuss controller strategies with better performance and lower consumption.

2.   Dynamics of the corticothalamic model

2.1.   The corticothalamic model generating generalized spike and wave discharges (SWD) of idiopathic generalized epilepsy (IGE)

Considerable literatures showed that SWD are generated by an interplay between the thalamus and cortex ^{[35, 36, 37]} by which normal thalamic discharges are sent to a slightly hyperexcitable cortex, which responds with spike and wave activity. Mathematical models of neural field dynamics provide a safer way to explore the effect of brain stimulations than in vivo experimentation. Several models of SWD have been proposed in literature to study mechanisms of SWD seizure genesis and dynamics ^{[38, 39, 40]}.

Neural population models describe the macroscopic neural activity that can be clinically recorded by an electroencephalogram (EEG) (see Figure 1(a)). Since most patient data is collected by EEG, which operates at the macroscopic scale, similar to the previous work, in this work we propose a clinically relevant thalamocortical circuit neural population model with random disturbance.

The original model from the work ^[27] and used in this paper is given as follows:

As seen from Figure 2, the cortical subsystem is composed of excitatory pyramidal ($PY$) and inhibitory interneuron ($IN$) populations. The thalamic subsystem includes variables representing populations of thalmocortical relay cells ($TC$) and neurons located in the reticular nucleus ($RE$).

All populations are interconnected in agreement with experimentally known connections using the connectivity parameters $C_{i} {(i = 1, 2, ...9)}$. $h_{py}, h_{in}, h_{tc}, h_{re}$ are input parameters, $\tau_{1}, \tau_{2}, \tau_{3}, \tau_{4}$ are timescale parameters, and $C_{1} = 1.8, C_{2} = 4, C_{3} = 1.5, C_{4} = 0.2, C_{5} = 10.5, C_{6} = 0.6, C_{7} = 3, C_{8} = 3, C_{9} = 1, h_{py} = -0.35, h_{in} = -3.4, h_{tc} = -2, h_{re} = -5, \tau_{1} = 26, \tau_{2} = 1.25\times26, \tau_{3} = 0.1\times26, \tau_{4} = 0.1\times26.$ $S{(x)}$ is the sigmoid function:

$L{(x)}$ is the linear function:

and $\epsilon = 250, 000, a = 2.8, b = 0.5.$

2.2.   Bistability dynamics

Figure 1(a) shows that under the parameters given above the system is at background resting state. However, when applied a pulse stimulus on the system at $10$ s, the system transits from a resting state to the SWD oscillations (see Figure 1(b)). We begin with the simplest of our scenarios. Figure 1(b) shows the maxima and minima of the model output for different values of the parameter $h_{tc}$. If $h_{tc}$ take smaller negative values ($h_{tc} < \approx-2$, left side of Figure 1(a)), there is only one stable equilibrium solution, all simulations converge to the steady state (stable focus). For less negative values ($-2 < \approx h_{tc} < -1.5$, shaded area of Figure 1(b)), a bistable region exists between the stable focus and the SWD oscillations. This arises following a fold of cycles bifurcation at $h_{tc}\approx-2$. The schematic diagram of bistability can be seen from Figure 1(c). The purple and yellow regions are the attractors of the stable focus and SWD, respectively. The stable focus can be considered analogous to resting state background EEG, and the high amplitude oscillatory attractor to be the seizure state. Figure 1(a) illustrates that stimulus can drive the system out of the attractor of stable focus and into the attractor of SWD. In the bistable region, a separating manifold (separatrix) exists between the two states in four dimensional state space. This manifold is highly complex in structure ^[22]. Transitions between non-seizure and seizure states can occur when a stimulus beyond the separatrix occurs. Beyond the disappearance of the stable focus (due to a subcritical Hopf bifurcation) at $h_{tc} >$ -1.5, monostable SWD and slow waves exist (right hand side of Figure 1(b)).

3.   A network model with disturbance and the design of a preview control method

In this section, we propose a corticothalamic neural network model with disturbance and simultaneously design the corresponding preview control method which can suppress the epileptic seizure SWD dynamics. Figure 3 gives the configuration of serve systems with optimal preview control.

3.1.   Control model with disturbance

Let

The time series of model output is defined as the mean value of $PY$ and $IN$ in variables:

where $C = \big(0.5 \quad 0.5 \quad 0 \quad 0 \big)$.

Then, based on the $(2.1)$, the CT model with both disturbance and control can be rewritten as

where

The control input of the system is

and

is the control matrix of $u(t)$. Considering different control strategies inputs, the number of elements $u_i$ in $u(t)$ is adjusted with the choice of strategy, while the number of $B_i$ in $B_0$ is adjusted to be consistent with the number of $u_i$. $B_i$ is a $4\times 1$ column vector.

The exogenous disturbance to the system is denoted by

And the exogenous disturbance matrix is taken as

The values in $D_0$ represent the influence degree of disturbance on this variable. The larger the value, the greater the influence of disturbance. This paper chooses the value that has a significant effect of simulation.

The nonlinear term of the $(2.1)$ can be written as

3.2.   System discretization

First, the discretization of system $(3.3)$ is carried out. When the sampling period is $\delta$ = 0.001 $ms$, the approximate discrete-time model is of the following form:

where

We take the reference signal $r{(k)}$ as the step signal in order to compare the action of the system when it is added to the controller, and the step value of $r{(k)}$ is considered the equilibrium value of the system output stabilization. Furthermore, the external disturbance $d{(k)}$ as the rectangular signal or random gaussian white noise. Since the preview controller can use known future information, the basic assumptions about the reference signal and the disturbance signal are required as follows:

Assumption 1 The preview length of the reference signal is $M_{r}$, which means that at time $k$, $r{(k)}$, $r{(k+1)}$, $r{(k+2)}$, ..., $r{(k+M_{r})}$ are available. The future values of the reference signal after time $k+M_{r}$ are constant with $r{(k+M_{r})}$, namely, $r{(k+M_{r}+j)} = r{(k+M_{r})}$, $j = 1, 2, 3, ...$.

Assumption 2 The preview length of the exogenous disturbance is $M_{d}$, i.e., at time $k$, the values $d{(k)}$, $d{(k+1)}$, $d{(k+2)}$, $...$, $d{(k+M_{d})}$ are available. After time $k+M_{d}$, the future values of the exogenous disturbance are equal to $d{(k+M_{d})}$, namely, $d{(k+M_{d}+j)} = d{(k+M_{d})}$, $j = 1, 2, 3, ...$.

The error signal $e{(k)}$ is defined as the difference between the output and reference signal, i.e.,

The purpose of the preview controller designed is that the output $y{(k)}$ can track the reference signal $r{(k)}$, asymptotically. In particular, we also aim to observe the effects of preview controllers and non-preview controllers when attacked by the random disturbances.

3.3.   Derivation of the augmented error system

We construct an augmented error system to transform the tracking problem of the original system into the regulation problem of the augmented error system and use the linear matrix inequality ($LMI$) technique to design a preview controller.

First, we introduce the first-order backward difference operator $\Delta$:

Applying the operator to both sides of $(3.10)$ leads to

where $\Delta f_{k} = f{(x{(k)})}-f{(x{(k-1)})}$ is the difference of the nonlinearity. Then, applying the same operator on the error signal

Furthermore, combining $(3.14)$ and $(3.15)$, we get

Combining the first equation of $(3.14)$ and $(3.16)$ yields

where

To introduce the preview information on reference signal and disturbance signal, we define new vectors:

From assumptions about the reference signal and the exogenous disturbance, it is easily seen that

where

Considering $(3.17)$ and $(3.21)$, we can obtain the augmented error system:

where

3.4.   Design of the preview controller

For $(3.23)$, we design a state feedback

where $K$ is a matrix to be determined. Substitute the controller into $(3.23)$, we get

The following lemma is needed to analyze conditions for asymptotic stability of the closed-loop system $(3.27)$.

Lemma ^[41] Supposing the Assumptions 1 and 2 hold, if there exist a positive definite matrix $P > 0$ and matrices $N$, $M$, $R$, and constant $\mu$ such that

where $\bar{F} = \begin{pmatrix} 0&I & 0 & 0 \end{pmatrix}$ represents the relationship between $\Delta x{(k)}$ and $\bar{x}{(k)}$, and $\gamma$ is the Lipschitz constant, i.e., for any $x, x' \in\textbf{R}^{n}$, the inequation

holds, then the closed loop system $(3.23)$ is asymptotically stable, where the state feedback gain matrix $K = RM^{-1}$.

By utilizing the MATLAB toolbox, the linear matrix inequality ($LMI$) $(3.28)$ in Lemma can be resolved. According to $f(x(t))$ in $(3.11)$, we can calculate the Lipschitz constant $\gamma$ of $f(x(t))$. Subsequently, the feasible solution $P$, $M$, $N$, $R$ can be obtained, which ensures the validity of the $LMI$ $(3.28)$. If $(3.28)$ is feasible, then the state feedback controller is $K = RM^{-1}$ that ensures the asymptotic stability of the closed-loop system (3.27). $K$ is divided as follows in accordance with the division of state variables in $(3.18)$, $(3.20)$, and $(3.24)$, i.e.,

After substituting $K$ into $(3.26)$, we can get

According to the definition of $\Delta u{(k)}$, the preview controller of system $(3.10)$ can be taken as the following theorem.

Theorem If there exist a positive definite symmetric matrix $M$ and a matrix $N$ such that $LMI$ $(3.28)$ holds, and $LMI$ $(3.28)$ has a feasible solution, then the preview controller of system $(3.10)$ is

Note that it is assumed that when $i < 0$, $x{(i)} = 0$, $u{(i)} = 0$, $r{(i)} = 0$. The coefficients in $u{(k)}$ are shown in $(3.30)$. Under this controller, the output of the system $(3.10)$ can asymptotically track target signal.

4.   Simulation results

Epileptic phenomena are generally modeled by a bistability of silence and seizure-like bursting. The stimulation disturbance is often able to induce epileptic seizures which can also be affected by random factors such as noise in the seizure process. Our aim is to design a controller that can control the system to normal state. In particular, it can maintain system ability to prevent secondary seizures when the system is disturbed again.

Observed from the analysis above, a preview controller is obtained based on the corticothalamic network dynamic model capable of generating the spike and wave discharges (SWD) of idiopathic generalized epilepsy (IGE). In fact, the control problem of SWD is essentially the tracking problem of the closed-loop system for the normal state.

We use parameters which place the model in bistability state. Thus, the appropriate stimulation disturbance can induce the SWD seizures. Considering that SWD control of the normal state is essentially a problem of tracking the normal state, we take the reference signal as the following step function:

We have assumed that the initial value of the state is $x(0) = (0.1724 \quad 0.1787 \quad -0.0818 \quad 0.2775)^T$. We select the initial values and then test that they are in the domain of attraction of the system stability. Moreover, we focus on the oscillatory case of the system and the case after being disturbed, which by its very nature is not affected by the initial values. Therefore, we do not provide an in-depth discussion of this detail. According to the system output value when the system is stable under this initial value, the tracking signal $r(k)$ takes on the value of 0.1755.

The exogenous disturbance $d_1(t) = d_2(t) = d_3(t) = d_4(t)$ is as follows

We take values of impulses and white noise for disturbances and verify that they cause seizures without the controller.The simulation results are shown in Figure 4.

We obtain a time length of $2300\; ms$ in Figure 4, during which the controller does not function. The results in Figure 4 demonstrate that our choice of interference is informative. To examine the difference before and after the controller takes effect, we will add the controller for the time period $t > 2300\; ms$ in following simulations.

Additionally, starting from the physical and actual meanings, we compared the effects of single-input controller and multiple-input controller on the model outputs while tracking the reference signal and receiving random disturbance. On this basis, we give some definitions of the coefficient matrix $B$ of the controller $u(t)$ in system $(3.10)$.

We take the dimension of $B$ as $4\times i$ (i = 1, 2, 3, 4), where $i$ indicates the number of input channels, and we assume that each input channel can only affect one nucleus in the model (i.e., $PY$, $IN$, $TC$, $RE$). The element of $B$ are 0 or 1, then if $B_{ji} = 1$ ($i, j$ = 1, 2, 3, 4), it means that the $i^{th}$ input channel has an effect on the $j^{th}$ nucleus. The number of element $1$ is equal to $i$.

Moreover, in order to observe the effect of preview references and disturbances on the tracking performance, we choose a preview length with good results through many experimental simulations. Then, two cases are considered, separately: ${(i)}\; M_{r} = 0\; ms$, $M_{d} = 0\; ms$ (no preview); ${(ii)}\; M_{r} = 3\; ms$, $M_{d} = 3\; ms$.

4.1.   Preview controller can completely abate the SWD discharges with effectively eliminating disturbance

From $t = 2300\; ms$ to $t = 5000\; ms$, the controller is applied to the system. We divide the simulation results into $4$ groups according to the dimension of $B$ and compare the effect of the controller with or without preview.

In the period $0 \le t < 500\; ms$, the system is in the non-epileptic seizure state. Then at time $t = 500\; ms$, the external disturbance $d(t) = 0.1$ leads to the seizure of the system, which is represented by the production of spike-wave discharges. After $t = 2300\; ms$, the controller starts to work and suppress seizures (i.e., track the reference signal). At this period, we apply a impulse disturbance $d(t) = 0.1$ and $d(t) = -0.1$ to the system. Additionally, random factors such as noise can also induce seizures (Figure 4(a)). For example, input from subcortical to cortical. This hypothesis is supported by the dynamics of the thalamic-cortical circuit and mechanical computational models of SWD seizure duration statistics obtained from rat and human models. We simulate the noise in the time period $3700\; ms \leq t \leq 4700\; ms$. The random disturbance is $d(k) = normrnd(0, 0.02, 1, 1001)$.

As shown in the Figure 5 about controllers with only one input channel, although seizures both can be suppressed by controller with or without preview, the preview controller can effectively eliminate or reduce the effect of secondary disturbance on the system.

The local zoomed-in graph in Figure 5(a) shows that the preview controller can use the information from the future reference signal with a time length $M_r = 3\; ms$ to control the system outputs to react earlier. The preview controllers can make the system shift more gently than the non-preview controller which has a more drastic change in the system output. This may cause patient discomfort in clinical applications, which is another reason why we considered applying the preview controller to epilepsy seizures. Furthermore, we also utilize the future information of the disturbance signal with a length of $M_d = 3\; ms$. The effect graph is similar to the local zoomed-in graph in Figure 5. In order to make the images more concise, we do not draw local zoomed-in diagrams for the other simulation result diagrams.

We also show the results of the other strategies to compare the effects of preview and non-preview controllers in the following Figures 6–8.

The results show that the controllers (without preview and with preview) can effectively terminate the epileptic seizure at $t = 2300\; ms$, and when the system receives secondary disturbances (impulses and white noise), the preview controller(red line in the figures) can make the system oscillate less or even no oscillations. Comparing controller effects in different strategies, we can find that the effect of the controllers containing $IN$ channel (the controllers without $PY$ channel) are worse than the effect of other controllers with $PY$ channel. We discuss $PY$ and $IN$ in terms of their physiological roles.

$PY$ expresses excitatory neurons and $IN$ indicates inhibitory neurons. The purpose of the controllers we designed is to regulate neurons. Based on the different physiological roles of $PY$ and $IN$ and the result of simulations, we speculate that $PY$ has a greater role in influencing epileptic seizures combined with the knowledge we already know. Therefore, the controller is better able to keep the system in a stable state when it functions on the $PY$ channel.

4.2.   Preview control can produces a lower performance consumption

Practically speaking, we would like the controllers to be effective in attenuating or terminating seizures while generating lower energy consumption, which would increase the feasibility of clinical applications. Therefore, we introduce a performance index function to evaluate the energy consumption of both preview controller and non-preview controller under different input control strategies. With minor modifications, we integrate our needs into the formula of the performance function $J$ in ^[42]. We only consider the tracking performance of the system, i.e., the tracking error $e(k)$ and the consumption of input $u(t)$ are both as small as possible. The performance index function $J$ for error system $(3.10)$ is defined as

where $Q_e$ and $Q_u$ are identity matrices of appropriate dimensions.

First, we compare the controller input values for different strategies. For inputs with many input channels, we compute the modulus length of $u(t)$, which is defined as

The input variation images for Figures 5–8 are provided in Figures 9–12.

We can tentatively draw a conclusion by comparing the controller input values for the different strategies described above. Depending on the number of input channels required, we can appropriately avoid $IN$ and $RE$ input channels, as they generate larger input values and thus higher consumption.

Then, we combine the partial statistics of the controller input values and the values of the performance function $J$ to filter the controller strategies for low energy consumption. We display the data in Table 1 (without preview) and Table 2 (with preview). The statistics in the table are derived from he time period when the controllers function.

At time $t = 2300\; ms$, the controllers adjust and control the pathological signal to reach the normal reference state, which lead to the input increase greatly. $Max1$ indicates the maximum value of the input when the controllers start to work. However, we consider the preview length of the reference signal is $M_{r}$, the preview controller will use the future information with $M_{r}$ length to adjust input value. Then, $Max1$ of preview controller is higher than the $Max1$ without preview. $Max2$ is the maximum value of input after reach the normal reference state for the first time. $Min$ is the minimum value of input from $t = 2300$ to $t = 5000$ and $Average$ represents the average value during this period.

According to the data in the tables, we can find that the performance consumption $J$ of the preview controller is significantly lower than that of the non-preview controller. Combining the analysis of the controller effects in Section 4.1 with the analysis of the controller input values and performance consumption, we can choose the preview controller strategy with $PY$ and $IN$ input channels. The performance consumption of this strategy is not the lowest, but it is the lowest compared to other strategies with great system output effects. Additionally, based on the conclusion of having large input values for the $IN$ and $RE$ channels, this strategy does not include both $IN$ and $RE$ channels. Moreover, in future practical clinical applications, the selected controller strategy must be modified to suit the situations of the individual patient. For example, is there a need to decrease the quantity of input channels, taking into account the potential effects of different channel numbers on the patient? Or is a lower system performance consumption needed to restore equilibrium?

5.   Conclusions

In this paper, a preview controller is designed based on the methods of discretization, augmented error system and linear matrix inequality (LMI), using, as example, the corticothalamic network dynamic model capable of generating the generalized spike and wave discharges (SWD) of idiopathic generalized epilepsy (IGE).By comparing the effect of the action of the controllers designed and the performance function $J$, we obtained the following conclusions:

1) Compared to non-preview controllers, the preview controller enables the system to restore the equilibrium state, i.e., terminate the seizure, more quickly and gently; similarly, the preview controller enables the system to maintain smaller oscillations after the system receives a secondary disturbance. Combined with Figure 2, the preview controller can make the output $y(t)$ of CT system Ⅰ has a lower impact on CT system Ⅱ.

2) We discuss the controller input designed for different strategies. According to the output $(3.2)$ of system, we analyze the advantages and disadvantages of the choice of strategies based on $PY$ and $IN$ neuronal nuclei from the perspectives of keeping system stability and performance consumption. Considering the possible negative effects of many input channels on the patient, we think that the preview controller with $PY$ and $IN$ input channels can provide both better control and lower consumption for clinical applications. The low consumption can reduce the loss of the controller in order to prolong the working time.

3) According to the data in Tables 1 and 2, the preview controller can produce overall lower performance consumption $J$. However, since the preview controller can utilize future information to make adjustments to the system at the current moment, it will have larger output values at individual moments, which is likely to cause physical discomfort to the patient. Therefore the preview length should not be selected too large and adjusted to a suitable length according to the actual needs.

To our knowledge, this is the first work investigating the preview control problem of epileptic seizures. This may be helpful to design clinically the robust and reliable seizure modulator.

Use of AI tools declaration

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

Acknowledgments

This research was supported by the National Natural Science Foundation of China (Grants Nos. 12072021, 12372061 and 12332004).

Conflict of interest

The authors declare that there are no conflicts of interest.