Seizure control in a neural mass model by an active disturbance rejection approach

Abstract

A closed-loop neuromodulation automatically adjusts stimuli according to brain response in real time. It is viewed as a promising way to control medically intractable epilepsy. A suitable closed-loop modulation strategy, which is robust enough to unknown nonlinearities, dynamics, and disturbances, is in great need in the clinic. For the specialization of epilepsy, the Jansen’s neural mass model is utilized to simulate the undesired high amplitudes epileptic activities, and active disturbance rejection control is designed to suppress the high amplitudes of epileptiform discharges. With the help of active disturbance rejection control, closed-loop roots of the system are far from the imaginary axis. Time domain response shows that active disturbance rejection control is able to control seizure no matter whether disturbances exist or not. At the same time, frequency domain response presents that enough stability margins and a broader range of tunable controller parameters can be obtained. Stable regions have also been presented to provide guidance to choose the parameters of active disturbance rejection control. Numerical results show that, compared with proportional-integral control, more accurate modulation with less energy can be achieved by active disturbance rejection control. It confirms that the active disturbance rejection control-based neuromodulation solution is able to achieve a desired performance. It is a promising closed-loop neuromodulation strategy in seizure control.

Introduction

Epilepsy is a typical chronic neurological disorder resulting from the development of synchronous firing in a massive group of neurons. It afflicts nearly 1% of the people worldwide. The newest research reports illuminate that almost 1/3 patients never achieved seizure freedom by antiepileptic drug therapy or even experience unexpected side effects. ¹ For those drug-resistant patients, only a small fraction qualifies for surgical resection of the epileptogenic zone. Patients, who are not suitable for surgery and who are still suffering seizures after surgeries, attract much attention not only from the physicians and surgeons but also from the social societies and governments. Since epilepsy imposes significant psychological, physical, and financial burdens on those patients and their families. In order to address it, numerous new and effective treatments have been proposed. In the last decades, neuromodulation technologies, such as vagus nerve stimulation, deep brain stimulation, and responsive neurostimulation, have been utilized to alleviate the intractable epilepsy. ² It has been confirmed that neuromodulation is one of the most impressive achievements in epilepsy treatments. Additionally, invasive, especially, noninvasive neurostimulation techniques will become more effective and popular in the next couple of years. ³

For the stimulation of nerve tissues, there are two modes. One is the open-loop way, that is, stimulus signals are generated from the device and delivered to the brain. Such mode provides an invariant train of stimulatory pulses in a fixed frequency. It is independent of ongoing and varying neuronal activities. In contrast to the scheduled stimulation, the other is the closed-loop mode. It delivers stimulatory pulses upon the detection of seizure activities, that is, it generates the stimulatory pulses according to the real-time neural activities. Such a system is a ‘clever’ system which is flexible and individually adaptable to patients’ specific ictal pattern. However, by comparison with the open-loop mode, it is technically more challenging.

For the closed-loop mode, the challenge is the ‘on-demand’ stimulation strategy. It is known that system control theory is just the method generating the control signals based on the real-time system outputs. From this point of view, system control theory is the proper approach to address such a challenge. In addition, based on the closed-loop system, high efficacy of neurostimulation, a lower number of stimuli, and fewer side effects can be obtained. ⁴ Therefore, the closed-loop mode may be a promising method to suppress epileptic seizures. However, for improvements in seizure suppression, a ‘clever’ stimulation strategy, that is, an advanced control algorithm, which is able to automatically and adaptively generate stimuli with appropriate amplitudes and frequencies, is in great need. In order to acquire improvements, much effort has been made. ^{1,5 –11} Proportional–integral (PI) control ⁵ and proportional–derivative (PD) control ⁶ have been designed to suppress epileptic activities simulated by a neural mass model (NMM). Relationship between parameters of PI/PD control and parameters of an NMM has been analyzed to provide insights into the mechanism of suppressing epileptic activities. Considering nonlinearities and uncertainties of a neural system, a fuzzy inference system is taken to acquire the parameters of proportional–integral–derivative (PID) control, that is, fuzzy PID control is designed to restrain the epileptiform spikes. ⁷ For minimizing control energy and estimating system states from noisy measurements, unscented Kalman filter (UKF) ^8,9 has been employed, and the UKF control strategy has been proposed for the NMM and the Wilson–Cowan equations. Based on the Pinsky–Rinzel model, input–output linearization (IOL) approach is designed to terminate epilepsy. ¹⁰ Reported approaches are mostly designed to suppress the epileptic activities simulated by those computational epileptic models. For an epileptic animal model, electroencephalogram (EEG) features are used to detect seizures, and stimulation for seizure elimination is applied in a predetermined mode with a fixed frequency and amplitude. ^{11 –14} It does greatly reduce the rate of ictal events, ⁴ however, such mode is not able to adjust stimulations in response to the specific electrical activities. Therefore, it cannot provide a personalized solution.

If a stimulation protocol is adaptive, and the closed-loop control algorithm is able to adjust its amplitudes and frequencies based on the real-time electrical activities, better control of seizure activities is able to achieve. Proportional feedback simulation has been tried in rats, ¹⁵ mean reduction in amplitude variance is significant. It shows that proportional feedback stimulation is a promising way to suppress seizure activities, and more complicated control algorithms for generating feedback stimulation may provide further improvements in seizure suppression. ^4,15 However, for improving the seizure suppression, larger feedback gains have been tested and it has been found that higher gains are unable to maintain the stability of a system. As a result, the animals used were killed. Obviously, stimulation parameters are determined by physicians or surgeons empirically, and it cannot guarantee the safety of animals or even the patients. It may be effective for some patients, but it still needs some improvements. Firstly, if stable range of a closed-loop neuromodulation system can be provided in advance, one can guarantee animals’ or even the patients’ safety. Additionally, in a stable range, one may get a group of proper stimulation parameters so as to get better modulation. However, such requirements cannot be satisfied in animal experiments or even clinic trials on epileptic patients. Solid guidance for both closed-loop neuromodulation and determining a group of better stimulation parameters can be obtained by numerical simulations with a computational model of epilepsy.

In this article, a computational model is also utilized, and closed-loop neuromodulation by an advanced control algorithm is designed. Jansen’s NMM is characterized by interactions of the interlinked excitatory and inhibitory feedback loops. It is able to successfully simulate epileptic activities that are similar to those observed in physiological experiments. Therefore, it has been commonly utilized in most of the studies. ^{5,6,16 –23} PI/PD control ^5,6 is able to get its control signals when tracking errors exist, however, it is a passive control approach. Such an approach depends on tracking errors to generate the control signals, it is susceptive to nonlinearities and uncertainties, and it needs more energy to modulate epileptic activities. Fuzzy PID ⁷ is robust, but fuzzy rules depend largely on experiences. Additionally, fuzzy rules lower system performance more or less. For the complexities and specialties of epilepsy, it is difficult to set a faithful model. Additionally, uncertainties and nonlinearities are ubiquitous in neuromodulation systems. Therefore, rather than PID, UKF, or IOL approaches, an advanced control algorithm, which is robust enough to those undesired factors, is able to deliver stimulation signals adaptively to real-time seizure activities, reduce the average number of stimulations and the adverse effects, and prolong the battery life, is in great need. On the other hand, closed-loop neurostimulation holds the greatest therapeutic potential, ²⁴ but there still exist numerous obstacles in realizing the closed-loop stimulations. A very conspicuous one is the contradiction between the computational complexity of a closed-loop neuromodulation algorithm and its realization. Current systems have limited processing abilities. It is not practical to run computationally complex algorithms smoothly, although those algorithms are promising. In other words, they are too complex to be realized in the real-time modulation. ²⁴

From the analysis above, a clear picture can be drawn. Closed-loop neuromodulation attracts much attention and it is increasingly utilized in the United States, ²⁴ however, neuromodulation with a fixed frequency and amplitude or realized via simple closed-loop control approaches need to be improved. Additionally, neuromodulation realized by complex algorithms is not practical for clinical applications. Therefore, a control algorithm, which can achieve the desired performance and relatively simple to be realized on a chipset, is of great importance.

In recent years, disturbance estimation and rejection approaches attract much more attention. For example, disturbance observer-based sliding model control has been designed for a class of underactuated robotic system ²⁵ and mobile wheeled inverted pendulum systems. ²⁶ A flexible spacecraft ²⁷ and polymer electrolyte membrane fuel cell systems ²⁸ have also employed disturbance observer-based control techniques. Active disturbance rejection control (ADRC), a special kind of disturbance estimation and rejection approaches proposed by Han in the 1990s, is a promising way to deal with nonlinearities and uncertainties. By virtue of an extended state observer (ESO), one can get system output, the derivatives of the system output, and the generalized disturbance. Then, a control signal can be constructed. Undesired factors, which affect system performance, are estimated and cancelled before they corrupt the system output. Therefore, ADRC is an active control approach and it is capable of providing proper neuromodulation stimuli according to real-time EEG signals. Moreover, it can be seen later that ADRC is simple enough to be realized. Therefore, in recent years, ADRC has been accepted and achieve satisfactory performance in numerous fields. ^{29 –32} From the engineering point of view, ADRC is the very suitable closed-loop neuromodulation solution for its practical design and satisfied performance. Main contributions are

ADRC has been firstly introduced in suppressing epileptic activities,

closed-loop performance has been confirmed via both time domain and frequency domain when excitatory and inhibitory gains of an NMM vary in a certain range,

stable region has been presented to provide a solid guidance to determine the parameters of ADRC in seizure control, and

advantages of ADRC over PI in seizure control have been confirmed from the time domain indexes, including root mean square error (RMSE), mean absolute error (MAE), integral of time-multiplied absolute-value of error (ITAE), and the energy expenditure in the neuromodulation.

In this article, the underlying reasons for epilepsy are analyzed by frequency analysis of the NMM. Then, ADRC is taken to suppress epileptiform activities generated by an NMM. Parameters and performance are discussed. The following parts are organized like this. The NMM and its performance are presented in the second section. ADRC is introduced in the third section. Seizure control by ADRC, root loci, stable regions, frequency domain performance, and numerical simulations are presented in the fourth section. Finally, concluding remarks are drawn in the fifth section.

The neural mass model and its performance

The neural mass model

Epileptiform activities, resulting from the lack of balance between the excitation and inhibition, can be described by an NMM. Its structure can be illustrated in Figure 1. ^6,16

Figure 1.

Scheme of the NMM. NMM: neural mass model.

As shown in Figure 1, an NMM model is composed of three coupling subpopulations, including the main subpopulation (in the middle of the structure), the excitatory feedback subpopulation (on the top of the structure), and the inhibitory feedback subpopulation (at the bottom of the structure). The main subpopulation comprises principal cells that receive feedback signals from both the excitatory and the inhibitory subpopulations.

The excitatory and inhibitory subpopulations transform the average presynaptic firing rates into the average postsynaptic membrane potentials. Each of them can be represented by following impulse responses, ^6,16 that is

$\left\{\begin{array}{l}{h}_e(t)=H_{e}t{e}^{-t}{}^{/{\tau }_e}/{\tau }_e\\ h_i(t)=H_{i}t{e}^{-t}{}^{/{\tau }_i}/{\tau }_i\end{array}\right.$

1

where H_e and H_i are excitatory and inhibitory average synaptic gains, and ${\tau }_e$ and ${\tau }_i$ are lumped representations of the membrane time constants and all other spatially distributed delays in the dendritic network. Parameters can be chosen to simulate the loss of homeostatic control generated by the lack of balance between the excitation and inhibition.

Based on impulse responses given in equation (1), dynamics of the excitatory and inhibitory subpopulations can be obtained and they can be described by two second-order ordinary differential equations

$\left\{\begin{array}{l}{\ddot{y}}_e(t)+2\frac{{\dot{y}}_e(t)}{{\tau }_e}+\frac{y_e(t)}{{\tau }_e^2}=\frac{H_e{i}_e(t)}{{\tau }_e}\\ {\ddot{y}}_i(t)+2\frac{{\dot{y}}_i(t)}{{\tau }_i}+\frac{y_i(t)}{{\tau }_i^2}=\frac{H_i{i}_i(t)}{{\tau }_i}\end{array}\right.$

2

where $i_e(t),i_i(t)$ and $y_e(t),y_i(t)$ are input and output signals of the excitatory and inhibitory subpopulations, respectively.

Nonlinear static function S presented in Figure 1 transforms the average membrane potential into the average firing rate, and it can be depicted by a sigmoid function

$S(v)=\frac{2e_0}{1+e^{r(v_0-v)}}$

3

where e ₀ determines the maximum firing rate, v ₀ is the postsynaptic potential when 50% firing rate is achieved, r is the steepness of the sigmoid function, and v is the average presynaptic membrane potential.

Interactions among main cells, excitatory and inhibitory interneurons are characterized by four connectivity constants C ₁ through C ₄. It includes the total number of synapses established by interneurons. An excitatory input is represented as p(t), which is utilized to simulate influences from neighboring areas and modeled by the white noise. ⁵ Output of an NMM, that is, y(t), is the average synaptic activities of the pyramidal cells, which can be interpreted as the EEG signals. ^5,6

Nominal values of those parameters can be found in Table 1. ^6,16

Table 1.

Interpretation and values of the parameters in an NMM.

Parameters	Description	Values
He	Average excitatory synaptic gain	3.25 mV
Hi	Average inhibitory synaptic gain	22 mV
τe	Average synaptic time constant for excitatory population	0.0108 s
τi	Average synaptic time constant for inhibitory population	0.02 s
C 1, C 2	Average number of synaptic contacts in the excitatory feedback loop	C 1 = 135; C 2 = 108
C 3, C 4	Average number of synaptic contacts in the inhibitory feedback loop	C 3 = 33.75; C 4 = 33.75
v 0, e 0, r	Parameters of the nonlinear S function	v 0 = 6 mV; e 0 = 2.5 s−1; r = 0.56 mV−1

NMM: neural mass model.

Characteristic of the NMM

In Figure 1, only sigmoid function S is nonlinear. In order to analyze the NMM by linear system theory, nonlinear function S is linearized. According to small perturbation around the equilibrium point, dynamics of S depends on its slope at the equilibrium point. It can be defined as the sigmoid ‘gain’ K _s

$K_s=\dot{S}(v){|}_{v=v_0}=\frac{2e_{0}r{e}^{(v_0-v)}}{{[1+e^{r(v_0-v)}]}^2}{|}_{v=v_0}=\frac{e_{0}r}{2}$

4

Remark 1

Linearization is utilized so as to present the root loci and frequency responses in the following sections. It does not affect the analysis of ADRC, since ADRC is less dependent on the faithful system model in its design and analysis. Additionally, the local linearization model also coincides with the properties of the nonlinear NMM model in terms of the influence of H_e s and H_i s.

Obviously, the local linearization approach is much simpler, and it is enough in the analysis. ^5,6 After the Laplace transformation, the transfer function of equation (2) can be obtained

$\left\{\begin{array}{l}{G}_e(s)=\frac{H_e{\tau }_e}{{({\tau }_{e}s+1)}^2}\\ G_i(s)=\frac{H_i{\tau }_i}{{({\tau }_{i}s+1)}^2}\end{array}\right.$

5

where $G_e(s)$ and $G_i(s)$ are transfer functions of $h_e(t)$ and $h_i(t)$ , respectively. Then, the NMM given in Figure 1 can be redrawn as in Figure 2.

Figure 2.

Scheme of the linearized NMM. NMM: neural mass model.

According to Figure 2, the transfer function of the NMM is

$G_{\text{NMM}}(s)=\frac{Y(s)}{P(s)}=\frac{G_e(s)}{1+K_s^2{G}_e(s)[C_3{C}_4{G}_i(s)-C_1{C}_2{G}_e(s)]}$

6

where Y(s) and P(s) are the Laplace transformation of y(t) and p(t), respectively.

Substitute $G_e(s)$ and $G_i(s)$ into equation (6), we have

$G_{\text{NMM}}(s)=\frac{N(s)}{D(s)}$

7

where

$\begin{array}{l}N(s)=H_e{\tau }_e[{\tau }_e{}^2{\tau }_i{}^2{s}^4+2({\tau }_e{}^2{\tau }_i+{\tau }_e{\tau }_i{}^2)s^3+({\tau }_e{}^2+4{\tau }_e{\tau }_i+{\tau }_i{}^2)s^2+2({\tau }_e+{\tau }_i)s+1]\\ D(s)={\tau }_e{}^4{\tau }_i{}^2{s}^6+(2{\tau }_e{}^4{\tau }_i+4{\tau }_e{}^3{\tau }_i{}^2)s^5+({\tau }_e{}^4+8{\tau }_e{}^3{\tau }_i+6{\tau }_e{}^2{\tau }_i{}^2)s^4\\ +(4{\tau }_e{}^3+12{\tau }_e{}^2{\tau }_i+4{\tau }_e{\tau }_i{}^2)s^3+(6{\tau }_e{}^2+8{\tau }_e{\tau }_i+{\tau }_i{}^2+K_s{}^2{H}_e{C}_3{C}_4{H}_i{\tau }_i{\tau }_e{}^3\\ -K_s{}^2{H}_e{}^2{\tau }_e{}^2{C}_1{C}_2{\tau }_i{}^2)s^2+(4{\tau }_e+2{\tau }_i+2K_s{}^2{H}_e{\tau }_e{}^2{C}_3{C}_4{H}_i{\tau }_i-2K_s{}^2{H}_e{}^2{\tau }_e{}^2{C}_1{C}_2{\tau }_i)s\\ +K_s{}^2{H}_e{\tau }_e{C}_3{C}_4{H}_i{\tau }_i-K_s{}^2{H}_e{}^2{\tau }_e{}^2{C}_1{C}_2+1\end{array}$

8

As we know, the unbalance between the excitation and inhibition results in epileptic activities. In the NMM, homeostasis can be achieved by regulating the excitation and inhibition. Abnormal values of average excitatory synaptic gain H_e or average inhibitory synaptic gain H_i generate the epilepsy. Here, by increasing H_e or decreasing H_i , we can see how H_e and H_i influence the dynamics of the NMM vividly from the root loci.

Note that the abnormal value of the excitatory gain is larger than the nominal value, and the abnormal value of the inhibitory gain is smaller than the nominal one. Therefore, with an attempt to see the influence of H_e and H_i clearly, firstly, let the excitatory gain H_e range from 3.25 to 10, and the inhibitory gain H_i still chooses the nominal value (i.e. $H_i=22$ ). Secondly, keep the excitatory gain H_e be nominal one (i.e. $H_e=3.25$ ), and the inhibitory gain H_i is changed from 0.5 to 21. Other parameters are taken to be the nominal values. The root loci of the NMM are shown in Figure 3.

Figure 3.

Root loci of the NMM. NMM: neural mass model.

Since the NMM is a six-order system, there are six branches. When H_e and H_i are abnormal, Figure 3(a) and (b) presents unstable branches. Obviously, abnormal excitatory gains or over excitations result in roots that are far from the imaginary axis in the right half s plane (see Figure 3(a)). It means larger excitation does lead to epileptic activities with higher probability. At the same time, decreasing inhibition also makes the roots of the NMM be still far from the imaginary axis (see Figure 3(b)). It also signifies that smaller inhibitory gains also result in epilepsy with higher probability.

For epilepsy, it results from the loss of homeostasis between the excitation and inhibition. Therefore, an approach, which is able to make the system be less sensitive to the excitatory gains and inhibitory gains, is in great need. Next, ADRC is introduced, and the neuromodulation will be checked.

Active disturbance rejection control

ADRC is a special approach, which is independent of concrete model of a system. It is able to deal with both tracking and disturbance rejection problem by its controller and ESO. For its satisfying performance, an increasing number of applications have been emerging. The basic structure of a closed-loop system by ADRC is presented in Figure 4.

Figure 4.

Structure of a closed-loop system by ADRC. ADRC: active disturbance rejection control.

In order to illustrate ADRC, following system is considered

$\ddot{y}=f(·)+bu$

9

Here, y is the system output, and $f(·)$ is the generalized disturbance, which depicts the dynamics that make the system be different from a pure series integrators. b is an unknown control coefficient. Control law u is designed as

$u\text{=}\frac{u_0-z_3}{b_0}=\frac{k_p(y_r-z_1)-k_d{z}_2-z_3}{b_0}$

10

z₁, z₂, and z₃ are outputs of the ESO, which estimate a system output, derivative of a system output, and the generalized disturbance, respectively. An ESO can be designed as

$\left\{\begin{array}{l}{\dot{z}}_1=z_2+{\beta }_1(y-z_1)\\ {\dot{z}}_2=z_3+{\beta }_2(y-z_1)+b_{0}u\\ {\dot{z}}_3={\beta }_3(y-z_1)\end{array}\right.$

11

where y_r is the set-point value; $k_p,k_d$ , and b₀ are tunable controller parameters; and ${\beta }_1,{\beta }_2,{\beta }_3$ are tunable observer parameters.

After substituting equation (10) into (9), we have

$\ddot{y}=\overline{f}(·)+k_p(y_r-z_1)-k_d{z}_2-z_3$

12

where $\overline{f}(·)\triangleq f(·)+(b-b_0)u$ .

Let estimation errors be $e_1=y-z_1,e_2=\dot{y}-z_2,e_3=\overline{f}-z_3$ , and tracking errors be $e=y_r-y$ , we have

$\ddot{y}=k_p(y_r-y)-k_d\dot{y}+k_p{e}_1+k_d{e}_2+e_3$

13

Obviously, when ESO performs well, that is, $e_1=y-z_1\approx 0,e_2=\dot{y}-z_2\approx 0,e_3=\overline{f}-z_3\approx 0$ , equation (13) can be rewritten as

$\ddot{y}\approx k_p(y_r-y)-k_d\dot{y}$

Accordingly, the desired closed-loop transfer function is

$G_{\text{cld}}(s)=\frac{k_p}{s^2+k_{d}s+k_p}$

14

By choosing proper k_p and k_d , desired closed-loop dynamics can be obtained. However, there are six parameters needed to be determined. To address it, let

$\begin{array}{l}{s}^3+{\beta }_1{s}^2+{\beta }_{2}s+{\beta }_3={(s+{\omega }_o)}^3\\ s^2+k_{d}s+k_p={(s+{\omega }_o)}^2\end{array}$

then the bandwidth-parameterization based tuning approach can be obtained, ³³ that is

${\beta }_1=\text{3}{\omega }_o,{\beta }_2=\text{3}{\omega }_o^2,{\beta }_3={\omega }_o^3,k_p={\omega }_c^2,k_d=2{\omega }_c$

15

Therefore, tunable parameters become ${\omega }_c,{\omega }_o$ , and b₀. From equations (10) and (15), we have

$u=\frac{{\omega }_c^2(y_r-z_1)-2{\omega }_c{z}_2-z_3}{b_0}$

16

Taking Laplace transformation for both sides of equation (11), we have

$\left\{\begin{array}{l}{z}_1(s)=\frac{3{\omega }_o{s}^2+3{\omega }_o^{2}s+{\omega }_o^3}{{(s+{\omega }_o)}^3}y(s)+\frac{b_{0}s}{{(s+{\omega }_o)}^3}u(s)\\ z_2(s)=\frac{(3{\omega }_o^{2}s+{\omega }_o^3)s}{{(s+{\omega }_o)}^3}y(s)+\frac{b_0(s+3{\omega }_o)s}{{(s+{\omega }_o)}^3}u(s)\\ z_3(s)=\frac{{\omega }_o^3{s}^2}{{(s+{\omega }_o)}^3}y(s)-\frac{b{\omega }_o^3}{{(s+{\omega }_o)}^3}u(s)\end{array}\right.$

17

Substituting equation (17) into (16), one has

$\begin{array}{l}u=\frac{1}{b_0}\frac{{(s+{\omega }_o)}^3}{{(s+{\omega }_o)}^3+2{\omega }_c{s}^2+({\omega }_c^2+6{\omega }_o{\omega }_c)s-{\omega }_o^3}[{\omega }_c^2{y}_r\\ -\frac{(3{\omega }_c^2{\omega }_o+6{\omega }_c{\omega }_o^2+{\omega }_o^3)s^2+(3{\omega }_c^2{\omega }_o^2+2{\omega }_c{\omega }_o^3)s+{\omega }_c^2{\omega }_o^3}{{(s+{\omega }_o)}^3}y]\end{array}$

18

Seizure control by active disturbance rejection approach

Control scheme

Seizure control by ADRC is designed. Based on equation (18), the equivalent structure of Figure 4 can be shown in Figure 5.

Figure 5.

Equivalent structure of a closed-loop system by ADRC. ADRC: active disturbance rejection control.

$G(s),H(s)$ given in Figure 5 can be depicted as

$\begin{array}{l}G(s)=\frac{1}{b_0}\frac{{(s+{\omega }_o)}^3}{{(s+{\omega }_o)}^3+2{\omega }_c{s}^2+({\omega }_c^2+6{\omega }_o{\omega }_c)s-{\omega }_o^3}\\ H(s)=\frac{(3{\omega }_c^2{\omega }_o+6{\omega }_c{\omega }_o^2+{\omega }_o^3)s^2+(3{\omega }_c^2{\omega }_o^2+2{\omega }_c{\omega }_o^3)s+{\omega }_c^2{\omega }_o^3}{{(s+{\omega }_o)}^3}\end{array}$

19

Here, seizure control is discussed. Epileptic activities, usually, are characterized by high-amplitude oscillations. ^{34 –36} To suppress the high-amplitude epileptic activities, similarly to literatures., ^1,5,6,8,23 the desired output is set to be zero, that is, y_r = 0. Then, the structure of seizure control by ADRC can be shown in Figure 6.

Figure 6.

Seizure control by ADRC. ADRC: active disturbance rejection control.

Remark 2

Here, ADRC is designed to overcome the influence of abnormal H_e s and H_i s. By an active control, the high amplitudes of epileptiform wave are expected to be suppressed. However, it is worth pointing out that the controlled EEG signals will not be zero, since EEG signals will not be zero even if H_e s and H_i s are normal. Setting value for seizure control is chosen to be zero so as to reduce the high amplitudes, and zero setting value can also be found in many other literatures (see literatures ^1,5,6,8,23 ).

Performance analysis

Root loci analysis

From Figure 6, the open-loop transfer function can be obtained.

$\begin{array}{l}{G}_{\text{op}}(s)=G(s)H(s)G_{\text{NMM}}(s)\\ =\frac{(3{\omega }_c^2{\omega }_o+6{\omega }_c{\omega }_o^2+{\omega }_o^3)s^2+(3{\omega }_c^2{\omega }_o^2+2{\omega }_c{\omega }_o^3)s+{\omega }_c^2{\omega }_o^3}{b_0[{(s+{\omega }_o)}^3+2{\omega }_c{s}^2+({\omega }_c^2+6{\omega }_o{\omega }_c)s-{\omega }_o^3]}{G}_{\text{NMM}}(s)\end{array}$

20

Obviously, if parameters of ADRC are fixed, root loci can be obtained. Here, let ${\omega }_c,{\omega }_o$ , and b ₀ be 500, 1300, and 500, respectively, root loci of the closed-loop system are given in Figure 7. Figure 7(a) and (c) presents the root loci when $H_i=22$ and $H_e\in [3.25,10]$ . From Figure 7(a), we can see that there is no branch in the right hand of s plane. It means that the closed-loop system is stable. Additionally, Figure 7(c) shows the details in an enlarged scale, it confirms that even if H_e is increasing, by ADRC, there still exists enough stability margins. Figure 7(b) and (d) gives out the root loci in the case of $H_e=3.25$ and $H_i\in [0.5,21]$ . Figure 7(b) shows all closed-loop roots are stable. In addition, Figure 7(d) describes the enlarged details of Figure 7(b). It presents that enough stability margins are guaranteed by ADRC no matter how H_i varies.

Figure 7.

Root loci of the closed-loop system by ADRC. ADRC: active disturbance rejection control.

Comparing with Figure 3, we can see that ADRC is able to correct the influence of the excitatory synaptic gain H_e and the inhibitory synaptic gain H_i . By ADRC, system dynamics will be more insensitive to H_e and H_i . It is the very aim of neuromodulation.

Stable region analysis

In order to determine the stable region of the closed-loop ADRC–NMM system in (b ₀, ω _c) parameter space, the Routh stability criterion is utilized. Based on the open-loop transfer function (20), we have the closed-loop characteristic equation

$\Delta (s)=1+G(s)H(s)G_{\text{NMM}}(s)=0$

21

equation (21) is a ninth-order equation. It is somewhat complex. Let ${\omega }_o=2.6{\omega }_c$ , however, the stable region can be obtained by Maple. Figure 8 gives out the green stable region when $H_e=4.5,H_i=22$ .

Figure 8.

Stable region of ADRC in suppressing epileptic activities generated by an NMM. ADRC: active disturbance rejection control; NMM: neural mass model.

For different abnormal H_e s and H_i s, stable regions have been discussed to check their influence. Let ${\omega }_o=2.6{\omega }_c,H_i=22$ , and $H_e=3.25$ , Figure 9 gives out several stable regions. By ADRC, the influence of the varied H_e s and H_i s is minimized. Similar to Figure 7(c) and (d), Figure 9(a) and (b) also confirms that the varied excitatory gains (H_e s) have a greater impact on the stability of the closed-loop system. The fact will also be verified in frequency analysis given in the section ‘Frequency analysis’.

Figure 9.

Stable regions of ADRC for different abnormal excitatory and inhibitory gains. ADRC: active disturbance rejection control.

Frequency analysis

From the sections ‘Root loci analysis’ and ‘Stable region analysis’, we can see that the closed-loop system is less sensitive to H_e and H_i by ADRC. In this section, an explanation for the root loci will be presented in the section ‘Root loci analysis’ from the point of frequency domain will be presented. Then, the influence of the tunable ADRC parameters will also be discussed.

Figure 10(a) and (c) shows that the closed-loop neuromodulation system is insensitive to the varying parameters H_e s. Gain margins and phase margins presented in Table 2 also confirm that enough stability margins can be guaranteed even if the excitatory synaptic gain H_e varies. On the other hand, when inhibitory synaptic gain $H_i\in [0.5,21]$ , Figure 10(b) and (d) shows that magnitude frequency curves and phase frequency curves are almost overlapped near the cutoff frequency and high-frequency band. Gain margins and phase margins listed in the last four columns of Table 2 do confirm the fact. Therefore, from Figure 10 and data given in Table 2, we can safely arrive at the conclusion that, with a group of fixed parameters (i.e. parameters of ADRC chosen in the section ‘Root loci analysis’), ADRC is robust enough even if H_e and H_i vary in a certain range. In addition, facts presented here also illustrate Figure 9 from the viewpoint of frequency domain.

Figure 10.

Frequency response for a group of fixed parameters of ADRC. ADRC: active disturbance rejection control.

Table 2.

Frequency indexes for a group of fixed ADRC parameters.

ωc = 500, ωo = 1300, b 0 = 500	Hi = 22				He = 3.25
	He = 3.25	He = 5	He = 7	He = 10	Hi = 0.5	Hi = 5	Hi = 10	Hi = 21
Gain margin	18.04	14.30	11.37	8.28	18.04	18.04	18.04	18.04
Phase margin	48.26	48.05	43.31	34.09	48.31	48.30	48.29	48.26

ADRC: active disturbance rejection control.

However, is just a group of parameters effective? Here, we also check the closed-loop neuromodulation performance when choosing different groups of parameters. Figure 11 gives out the frequency responses corresponding to groups of controller bandwidths and observer bandwidths. Data presented in Table 3 also confirm that not only a group of fixed parameters given in the section ‘Root loci analysis’ is effective in suppressing epilepsy, but also groups of tunable ADRC parameters take effect. Especially, gain margins and phase margins are very close to each other. It signifies that ADRC is capable of taking effect in a wider range. It does help ADRC to apply in clinical practice.

Figure 11.

Frequency response for groups of parameters of ADRC. ADRC: active disturbance rejection control.

Table 3.

Frequency indexes for groups of ADRC parameters.

He = 7, Hi = 22, b 0 = 500	ωo = 1200				ωc = 500
	ωc = 430	ωc = 600	ωc = 770	ωc = 940	ωo = 980	ωo = 1100	ωo = 1220	ωo = 1340
Gain margin	11.62	10.90	10.49	10.26	11.08	11.18	11.29	11.41
Phase margin	44.91	41.51	39.54	38.35	43.47	43.31	43.28	43.34

ADRC: active disturbance rejection control.

Steady-state performance analysis

In this section, steady-state error of the closed-loop system is discussed. According to Figures 5 and 6, the tracking error transfer function from the input signal to tracking error can be derived, that is

$G_e(s)=\frac{E(s)}{R(s)}=\frac{1}{1+G(s)G_{\text{NMM}}(s)H(s)}$

22

where $G_e(s)$ is the transfer function of tracking error, and E(s) and R(s) are Laplace transformations of the tracking error e(t) and set-point value r(t), respectively.

When the set-point value $r(t)=c$ (c is a constant), we have $R(s)=\frac{c}{s}$ . Then

$\begin{array}{l}E(s)=R(s)G_e(s)=\frac{R(s)}{1+G(s)G_{\text{NMM}}(s)H(s)}\\ =\frac{c}{s}\cdot \frac{1}{1+G(s)G_{\text{NMM}}(s)H(s)}\end{array}$

23

According to final value theorem, the steady-state error is

$\begin{array}{l}{e}_{ss}=\underset{t\to \infty }{\text{lim}}e(t)=\underset{s\to 0}{\text{lim}}sE(s)\\ \begin{array}{c}\begin{array}{c}\end{array}\end{array}\begin{array}{c}\end{array}=\underset{s\to 0}{\text{lim}}s\cdot \frac{c}{s}\cdot \frac{1}{1+G(s)G_{\text{NMM}}(s)H(s)}\\ \begin{array}{c}\begin{array}{c}\end{array}\end{array}\begin{array}{c}\end{array}=\underset{s\to 0}{\text{lim}}\frac{c}{1+\overline{G}(s)}\\ \begin{array}{c}\begin{array}{c}\end{array}\end{array}\begin{array}{c}\end{array}=c\cdot {\overline{e}}_{ss}\end{array}$

24

where

$\begin{array}{l}{\overline{e}}_{ss}=\underset{s\to 0}{\text{lim}}\frac{1}{1+\overline{G}(s)},\\ \overline{G}(s)=G(s)G_{\text{NMM}}(s)H(s)\\ =\frac{(3{\omega }_c{}^2{\omega }_o+6{\omega }_c{\omega }_o{}^2+{\omega }_o{}^3)s^2+(3{\omega }_c{}^2{\omega }_o{}^2+2{\omega }_c{\omega }_o{}^3)s+{\omega }_c{}^2{\omega }_o{}^3}{{(s+{\omega }_o)}^3+2{\omega }_c{s}^2+({\omega }_c{}^2+6{\omega }_c{\omega }_o)s-{\omega }_o{}^3}\cdot \\ \frac{H_e{\tau }_e{({\tau }_{e}s+1)}^2{({\tau }_{i}s+1)}^2}{{({\tau }_{e}s+1)}^4{({\tau }_{i}s+1)}^2+K_s{}^2{H}_e{\tau }_e[C_3{C}_4{H}_i{\tau }_i{({\tau }_{e}s+1)}^2-C_1{C}_2{H}_e{\tau }_e{({\tau }_{i}s+1)}^2]}\end{array}$

Since ${\overline{e}}_{ss}=\underset{s\to 0}{\text{lim}}\frac{1}{1+\overline{G}(s)}=0$ , then the steady-state error $e_{ss}=0$ . It signifies that system output equals to the desired out in steady state no matter what set-point value c is. Therefore, in this article, to suppress the epileptic activities, generally let $r(t)=0$ , and the steady-state error still be zero. In other words, a neuromodulation with zero steady-state error can be obtained. However, it should be noted that noise signal exists in the NMM, therefore, small control error does exist in the closed-loop neuromodulation. It will be also illustrated clearly in the section ‘Simulations’.

Simulations

In this section, time domain responses are presented to illustrate that ADRC is able to suppress epileptic activities effectively. A group of hyper-excitatory parameters (i.e. H_e = 7, H_i = 22) is employed to simulate the epileptic behaviors. ADRC is designed to modulate the undesired activities. PI, a commonly used controller in control engineering, is also taken to make a comparison. The control action is activated at 8 s, and simulations last for 16 s. Two scenarios are considered. The first is the nominal case, and the second is the case when sinusoidal disturbance exists. Control structure is shown in Figure 12.

Figure 12.

Seizure control structure.

In both scenarios, the tunable parameters of PI and ADRC are selected from Table 4.

Table 4.

Tunable parameters of PI and ADRC.

	ωc /Kp	ωo /Ki	b 0
PI	310	2	–
ADRC	320	420	800

ADRC: active disturbance rejection control; PI: proportional-integral.

Scenario I

Here, a nominal case, that is, d = 0 in Figure 12, is considered. Both PI and ADRC are designed to control seizure. Figure 13 shows the comparisons of time domain responses and control signals between PI and ADRC.

Figure 13.

Time responses and control inputs of PI and ADRC. PI: proportional-integral; ADRC: active disturbance rejection control.

Comparisons of the indexes are listed in Table 5.

Table 5.

Comparisons of indexes between PI and ADRC.^a

	RMSE	MAE	ITAE	Energy
PI	0.2163	0.1502	28.7132	4.7220
ADRC	0.0582	0.0097	1.6993	4.4849
Improvement	73.09%	93.54%	94.08%	5.02%

ADRC: active disturbance rejection control; PI: proportional-integral; RMSE: root mean square error; MAE: mean absolute error; ITAE: integral of time-multiplied absolute-value of error.

^a $\text{RMSE}=\sqrt{\frac{1}{k}{\displaystyle \sum _{i=1}^k{e}_i^2}},\text{MAE}=\frac{1}{k}{\displaystyle \sum _{i=1}^k|e_i|},\text{ITAE}=T{\displaystyle \sum _{i=1}^k{t}_i|e_i|},E=\frac{T}{k}{\displaystyle \sum _{i=1}^k{u}_i^2}$ , e_i is the tracking error at the i-th time, and T is the sample time.

Figure 13(a) shows that ADRC is able to control seizure with much shorter transient time. In other words, EEG signals modulated by ADRC are better than the ones modulated by PI. Moreover, Figure 13 signifies that, with less oscillations of the control signal, better modulation can be obtained by ADRC. It means that, in seizure control, less energy is needed by ADRC. RMSE, MAE, ITAE, and control energy listed in Table 5 also show that, with less energy, ADRC is able to get much smaller RMSE, MAE, and ITAE. It signifies that ADRC can control seizure more effectively.

Additionally, it is worth pointing out that the energy expenditure of the stimulation signal directly affects the life of the stimulation battery. Therefore, energy expenditure is a very important factor, and it is also considered in the comparison. Both Figure 13 and Table 5 tell us that, in a nominal scenario, ADRC is better than PI in suppressing epileptic activities.

Scenario II

With an attempt to check the robustness of the neuromodulation, sinusoidal signal, that is, d = 120sin(πt), is introduced at the 12th second to simulate the external disturbance. Both PI and ADRC are designed, and the tunable parameters are also selected from Table 4. Figure 14 shows the time domain responses and control signals of PI and ADRC, respectively.

Figure 14.

Time responses and control inputs of PI and ADRC (sinusoidal disturbance exists). PI: proportional-integral; ADRC: active disturbance rejection control.

For the sinusoidal disturbance given in Figure 15, ADRC is capable of estimating and cancelling it effectively. Therefore, from Figure 14(a), we can see clearly that the time response of ADRC is much better than the one of PI. RMSE, MAE, and ITAE values of PI and ADRC given in Table 6 also confirm the fact mentioned above.

Figure 15.

Sinusoidal disturbance.

Table 6.

Comparisons of the indexes between ADRC and PI (sinusoidal disturbance exists).

	RMSE	MAE	ITAE	Energy
PI	0.2593	0.1587	30.1558	6.6909
ADRC	0.0616	0.0156	3.0274	6.2884
Improvement	76.24%	90.17%	89.96%	6.02%

ADRC: active disturbance rejection control; PI: proportional-integral; RMSE: root mean square error; MAE: mean absolute error; ITAE: integral of time-multiplied absolute-value of error.

Additionally, from Figure 14(b) and the last column data shown in Table 6, it is obvious that the desired performance can also be achieved by ADRC with less control energy. It is the very target that we want in the neuromodulation of epilepsy.

So far, the results of both scenarios tell us that ADRC provides much better performance in seizure control. Therefore, it may be a more promising approach to suppress epileptic activities.

Discussion

An NMM provides a straightforward way to describe the activities of the populations of neurons. It cannot describe exactly how neurons interact, however, it is suitable to depict the role of exhibition and inhibition, and the balance/unbalance between them from the macroscopic view. ¹⁹ Based on such a computational model, which is also employed in most studies, ^{5,6,16 –23} an active disturbance rejection modulation approach is discussed. Here, in order to verify the robustness of the closed-loop system, excitatory gains are varied within ±5% of its nominal value. Both PI and ADRC are taken the same parameters listed in Table 4. Monte Carlo tests have been performed, and the results are presented in Figure 16.

Figure 16.

Monte Carlo test results in the presence of H_e varies.

From Figure 16, one can see that, even if H_e varies in a certain range, RMSE values of ADRC are much less than the ones of PI.

Considering that an epileptic model is somewhat conservative, stable regions and frequency responses are also discussed when excitatory and inhibitory gains vary. System frequency domain performance is verified when different groups of parameters of ADRC are taken. Additionally, the time domain responses, in the presence and absence of external disturbance, are also discussed. Both the frequency responses and the time domain responses confirm that ADRC is robust enough to the nonlinearities, uncertainties, and external disturbances. In other words, there are enough stability margins for ADRC to deal with the conservative of a computational model. As a guidance for clinical applications, it may be a suitable closed-loop neuromodulation approach in the near future, if it is confirmed by animal experiments and clinical trials.

Conclusion

In this article, a closed-loop neuromodulation approach is discussed. Jansen’s NMM is taken to be the test platform, and the active disturbance rejection approach is employed to suppress the epileptic activities. Both time domain responses and frequency domain responses confirm that ADRC is capable of dealing with undesired factors, such as nonlinearities, uncertainties, and external disturbances, so as to guarantee the satisfied modulation. It is worth pointing out that ADRC is independent of a concrete model. Thus, it can be naturally developed to other neural models, animal models, or even clinical trials. From this point of view, the results presented in this article are able to provide a helpful reference, guidance, or even solution to the closed-loop neuromodulation in suppressing epilepsy. The future work should be online seizure control to standard animal epileptic models by ADRC or even applying ADRC to alleviate epilepsy patients as a long-term goal.