Suppressing epileptic activity in a neural mass model using a closed-loop proportional-integral controller

Closed-loop control is a promising deep brain stimulation (DBS) strategy that could be used to suppress high-amplitude epileptic activity. However, there are currently no analytical approaches to determine the stimulation parameters for effective and safe treatment protocols. Proportional-integral (PI) control is the most extensively used closed-loop control scheme in the field of control engineering because of its simple implementation and perfect performance. In this study, we took Jansen’s neural mass model (NMM) as a test bed to develop a PI-type closed-loop controller for suppressing epileptic activity. A graphical stability analysis method was employed to determine the stabilizing region of the PI controller in the control parameter space, which provided a theoretical guideline for the choice of the PI control parameters. Furthermore, we established the relationship between the parameters of the PI controller and the parameters of the NMM in the form of a stabilizing region, which provided insights into the mechanisms that may suppress epileptic activity in the NMM. The simulation results demonstrated the validity and effectiveness of the proposed closed-loop PI control scheme.


Methods
Model. The schematic and block diagram of the NMM is shown in Fig. 1A,B 44 , respectively. The NMM was composed of the following three interacting subpopulations: the main subpopulation (middle part), the excitatory (top part) and inhibitory (bottom part) feedback subpopulations. C 1 , C 2 , C 3 and C 4 are connectivity constants that represent interactions between the subpopulations and characterize the average numbers of synaptic contacts; p(t) is the input of the NMM and modelled by Gaussian noise; and y(t) is the output of the NMM, corresponding to the average synaptic activity of the pyramidal cells, that is interpreted as an EEG signal.
Each subpopulation of the NMM was composed of an excitatory, h e (t), or inhibitory synaptic dynamic function, h i (t), and a sigmoid static function, S(v). The synaptic functions, h e (t) and h i (t), transform the average pre-synaptic firing rates into average post-synaptic membrane potentials and are defined as S v e e ( ) 2 / [1 ] r v v 0 ( ) 0

46
, where 2e 0 is the maximum firing rate, v 0 is the post-synaptic potential corresponding to a firing rate of e 0 , and r describes the steepness of the sigmoid function.
The parameter values of the NMM are listed in Table 1 44 .
Graphical stability analysis method. A graphical stability analysis method was employed to determine the stabilizing region of PI controller for suppressing epileptic activity in the NMM. The key point of designing the PI controller was to choose the control parameter values such that homeostasis of the NMM was maintained. To achieve this objective, we wanted to determine the stabilizing region of the PI controller.
If it is assumed that the characteristic polynomial of one control system is denoted as Δ (s), then all of the coefficients of Δ (s) are real; therefore, the characteristic roots of Δ (s) must be complex conjugates. According to linear stability theory 47 , the stability boundary of the control system is defined by Δ (jω) = 0. Hence, it is sufficient to consider the following two cases: ω = 0 and ω ∈ (0, ∞ ). The stability boundary of the control system was defined as follows 48 : where Ω 1 and Ω 2 represent the control parameter sets of the employed controller to be determined for the two cases ω = 0 and ω ∈ (0, ∞ ), respectively. In this study, Equation (1) was utilized to determine the stabilizing region of the PI controller in the control parameters space.

Results and Discussion
Control Scheme. The control scheme was developed based on the fact that epileptic activity can be characterized as high-amplitude limit cycle oscillation born in Hopf bifurcation 4,28,31,33,36,41,42,[49][50][51][52][53] , which indicates that the fixed point of the NMM lost its stability. Thus, we aimed to design a PI controller to stabilize the unstable, fixed point to prevent the generation of the Hopf bifurcation and to further suppress high-amplitude epileptic activity.
Epilepsy is thought to be caused by an imbalance between excitation and inhibition resulted from hyper-excitation or low inhibition 32 . Thus, the design objective was to determine the PI controller parameter values that stabilized the unstable NMM caused by abnormally large excitatory or small inhibitory parameter values 32,42 such that homeostasis of the NMM was maintained and the high-amplitude epileptic activity was suppressed. To do that, we employed a graphical stability analysis method that determined the stabilizing region of the PI controller in the parameter space.
The proposed PI control scheme for suppressing epileptic activity is shown in Fig. 2A, where u(t) is the output of PI controller, corresponding to the electrical stimulation signals, and y(t) represents the local field potential of a neural mass. The local electrical field in the neural mass was measured by the recording electrode and feeds back to the PI controller via the stimulating electrode. By using the NMM to model activity of the neural mass in Fig. 2A, we derived the block diagram of the proposed control scheme, shown in Fig. 2B, where y(t) is the output of the NMM and u(t) is the output of the PI controller.
To take advantage of the graphical stability analysis method used to design the PI controller, we replaced the nonlinear sigmoid function, S(v), with its linear approximation around the equilibrium point, v = v 0 , as follows: According to Fig. 2B, we derived the transfer function of the NMM as follows: are the Laplace transform of h e (t) and h i (t) in the NMM, respectively. Figure 2B can be simplified as a standard form of typical control systems, shown in Fig. 2C, where r(t) is the desired output of the NMM, e(t) is the error signal of the closed-loop control system and the system is defined as e(t) = r(t) − y(t). For simplicity, in the remainder of the article, we will use PI-NMM to represent the closed-loop control system composed of the PI controller and the NMM.
In addition to epileptic activity, the NMM may feature multiple parallel states, such as low-amplitude oscillation and fluctuation output at a fixed point 33  only when the amplitude of the output of the NMM exceeded a predetermined threshold. Correspondingly, the control signal u(t) is defined as where K p and K i are the proportional and integral coefficients of the PI controller, respectively. A 0 is the predetermined amplitude threshold beyond which stimulation is applied. By conducting a Laplace transform of Equation (3), we can obtain the transfer function of the PI controller as follows: where E(s) is the Laplace transform of the error signal, e(t).
According to feedback control principles 47 , the control objective of a feedback control system is to make the error signal e(t) → 0. This indicates that the output of the PI-NMM control system will approach the desired output under the control of the PI controller, i.e., y(t) → r(t). In this study, our goal in implementing the PI controller was to suppress the high-amplitude epileptic activity in the NMM; thus, the desired output of the PI-NMM control system was set to zero, i.e., r(t) = 0.
Theoretical results of the stability analysis. To determine the stabilizing region of the PI controller in the (K p , K i ) parameter space in which the PI-NMM control system was stable, we used the graphical stability analysis method to determine the stability boundary of the PI controller in the K p , K i plane.
According to Fig. 2C, we can derive the characteristic equation of the PI-NMM control system as

pi NMM
where Δ (s) is the characteristic polynomial of the PI-NMM control system. Defining s = jω and substituting it into Equation (5), one can obtain pi NMM Equation (6) defines the stability boundary of the PI-NMM control system. According to the graphical stability analysis method shown in Equation (1), we can discuss the two cases ω = 0 and ω > 0.
Case 1: ω = 0. Substituting Equations (2) and (4) into Equation (5), we can derive the characteristic equation of the PI-NMM control system as follows Let s = jω. Thus, Equation (7) can be rewritten as For ω = 0, Equation (8) can be simplified as Thus, we can obtain the stability boundary for the case ω = 0, i.e., H 1 , as follows: This equation is a line in the (K p , K i ) plane that is shown in Fig. 3 (the purple line). Next, according to Equation (7), we can further derive the characteristic equation of the PI-NMM control system as follows: Note that the coefficient of s 7 in the characteristic polynomial, i.e., τ τ e i 4 2 , is positive. According to the Routh stability criterion 47 , if the system is stable, then all of the coefficients of the characteristic equation should have the same sign. The constant term should also be positive, i.e., K i H e τ e > 0. Thus, it follows that i Therefore, we can obtain the following lemma 1. Lemma 1 For the case ω = 0, the region above the K i axis in the (K p , K i ) plane is the stabilizing region of the PI controller.
Case 2: ω = 0. Supposing that δ R_NMM (ω) and δ I_NMM (ω) are the real and the imaginary components of G NMM (jω), respectively, then Equation (6) can be rewritten as  (11) can be split into the following two parts where δ R (ω) and δ I (ω) are the real and imaginary components of Δ (jω), respectively. Note that both δ R (ω) and δ I (ω) are dependent on K p and K i . Thus, the stability of the characteristic Equation (12) can be investigated in the parameter space (K p , K i ).
According to Equation (12), we can further derive where is the complex modulus of G NMM (jω).  4.5 mV). The arrows represent the direction of the curve (K p (ω), K i (ω)) in which ω increases. The light blue curve is defined by Equation (13), and the purple line is defined by Equation (9). The parameter space to the right of the curve and above the line is the stabilizing region of the PI controller.
According to Equation (1), the stability boundary is composed of the stability boundary line defined by Equation (9) and the stability boundary curve defined by Equation (13). As illustrated in Fig. 3, the stability line and the stability boundary curve divide the (K p , K i ) plane into two regions, which are denoted R 1 and R 2 , respectively. In the following section, we investigate which is the stabilizing region of the PI controller.
Next, we introduce the following proposition 56 . If one travelled along the curve defined by H 2 in the direction of increasing ω, then the right side is the stabilizing parameter region where det J < 0 and the left is where det J > 0. Here, J is the Jacobian matrix defined as R p R i I p I i From Equation (12), we obtain Substituting the above four equations into Equation (14), we can obtain Therefore, we derive the following lemma 2. Lemma 2 For the case ω > 0, following the curve (K p (ω), K(ω)) in the direction of increasing ω, the parameter space to the right of the curve defined by Equation (13) is the stabilizing region of the PI controller.
Main theoretical results. Combining Lemma 1 with Lemma 2 yields the following theorem.
Theorem 1 For the PI-NMM control system depicted in Fig. 2, the exact stabilizing parameter region of the PI controller is the parameter space above the line defined by Equation (9) and to the right of the curve defined by Equation (13) if it is followed in the direction of increasing ω.
Thus, according to Theorem 1, the region denoted by R 1 in Fig. 3 was identified as the stabilizing region of parameters (K p , K i ), and R 2 was the destabilizing parameter region.

Effect of excitatory and inhibitory NMM parameters on the stabilizing region of the PI controller.
Epileptic activity is caused by the imbalance of excitation and inhibition in the NMM, which is caused by extremely large excitatory, H e , or small inhibitory parameters, H i , respectively. Therefore, we discussed the effect of the two parameters on the stabilizing region of the PI controller. Figures 4 and 5 illustrate the stabilizing regions of the NMM PI controller for different abnormal values of H e and Hi, respectively. Figures 4B and 5B show the stabilizing region of the PI controller for small values of the integral coefficient, K i . The results demonstrate that the stabilizing region of the PI controller becomes smaller with increasing H e (hyper-excitation) and decreasing H i (low inhibition).
Theoretical results of the steady-state performance analysis. According to Fig. 2C, we can derive the closed-loop transfer function of the PI-NMM from the input, r(t), to the error signal, e(t), as where E(s) and R(s) are the Laplace transforms of e(t) and r(t), respectively. According to Equation (16), one can obtain In this study, the desired output of the PI-NMM was set to a constant value, r(t) = a, and its Laplace transform is derived as follows: Substituting Equations (16) and (18) into Equation (17), we further obtain which shows that the output, y(t), of the NMM is equal to the desired output, r(t), at steady state.
In this study, we set r(t) = a = 0; in this case, Equation (21) still holds. This result indicates that the PI controller can achieve an error-free control performance at steady-state and, thus, successfully suppresses epileptic activity in the NMM. It should be noted that the input of the NMM was noise signal; therefore, the PI-NMM control system could not be operated at steady-state 57 and a small control error existed, as shown in Figs. 6 and 7.
Simulation results. Simulations were conducted to illustrate the efficiency of the proposed PI control scheme to suppress epileptic activity in the NMM. We simulated the following two cases: hyper-excitation (H e = 7.0 mV) and low inhibition (H i = 17.0). According to Figs. 4 and 5, we determined the parameter values of the PI controller to be K p = 310 and K i = 2 for the hyper-excitation case and K p = 90 and K i = 2 for the low-inhibition case, respectively. The simulation results are illustrated in Figs. 6 and 7. The results showed that the output of the NMM without the PI controller was high-amplitude epileptic activity, which became low-amplitude activity under the control of the PI controller. Thus, the high-amplitude epileptic activities were successfully suppressed by the designed PI controller.

Limitations.
To determine the stabilizing region of the PI controller, we used the graphical stability analysis method and conducted a linearized approximation of the NMM sigmoid function, which is extensively used in previous NMM studies 27,46,[58][59][60] . It should be noted that the linearized approximation may provide a conservative estimate of the PI controller stabilizing region. In the future, a non-linear method, such as the bifurcation analysis method 28,42,43,45 , should be employed to determine a more specific stabilizing region for the proposed PI-based control scheme.

Conclusions
In the present study, we used Jansen's neural mass model as a test bed to develop a systematic design approach to determine the control parameters of a PI controller. It should be noted that the proposed design method of the PI controller was independent of a specific model. Thus, it can be extended to other neural models 37,40,45,46,[50][51][52][53][54][55] , such as the epileptic activity model developed by Wendling et al. 40 , and other low-order closed-loop controllers [24][25][26][27][28][29][30][31] , such as the proportional controller and the proportional-derivative controller. Here, we took a parametric approach to seizure behaviour to demonstrate how to use the PI controller to suppress epileptic activity. However, there are other factors that cause epilepsy 1 in addition to the excitatory and inhibitory parameters, such as the stimulus 43,45 . The design method and process were not dependent on the factors that cause epilepsy; thus, the proposed design method can be extended to other cases involving epileptic activity that is characterized as high-amplitude oscillations 1,4,28,40,41,45,49,52 .
A graphical stability analysis method was utilized to determine the stability region of the PI controller in the parameter space. This provided a region of the PI control parameters that would suppress epileptic activity in the NMM. The proposed method ensured that the design of the controller was analytical, enabled theoretical analysis and revealed cause and effect relationships in a theoretical manner. This allowed us to explore the relationship  between control parameters and model parameters that induced epileptic activity, which helped us understand the mechanism that suppresses neural diseases through closed-loop neurological electrical stimulation. In future work, we should attempt on-line seizure suppression by applying the proposed control scheme to standard animal models of epilepsy, with the long-term goal of applying it to human patients.