Epilepsy is the third most common neurological disorder right after Alzheimer's disease and stroke, affecting approximately 50 million people around the world. Epilepsy is characterized by seizure-like high amplitude oscillations of global brain activities affecting various mechanisms on microscopic and macroscopic scales. It is a complex set of disorders that can involve many areas of the cortex as well as underlying deep-brain systems.^[1] At the system level, epilepsy is a dynamical disorder of the brain, which makes it particularly suitable to be studied from the perspective of computational modeling and dynamical-system theory.^[2]

In this context, computational models have become attractive due to their ability to model complex neurological phenomena with relative ease by modeling the effects of various neurological factors.^[2] Over the past few years, substantial progress has been made in modeling epilepsy at levels ranging from microscopic to macroscopic.^{[3, 4]} At the macroscopic level, the neural population model (NPM) is developed to relate macroscopic phenomena, at the level of the electroencephalogram (EEG), with oscillations in smaller neuronal populations. The neural population model is based on a biologically plausible parameterization of the dynamic behavior of the layered neocortex, which starts from the fact that the EEG is a reflection of ensemble dynamics arising from interconnected populations of pyramidal cells and interneurons.^[5] The NPM was initially proposed to study the origin of the alpha rhythm, ^[6] and subsequently improved and extended to describe more general cortical electrical activities.^[7] Recently, NPMs have been extensively used to generate EEG, magnetoencephalography (MEG), and functional magnetic resonance imaging (FMRI) signals.^[8] Specifically, an NPM has been successfully used to generate an epileptic behavior similar to that observed experimentally, as well as the transition between normal and epileptic activities.^[3] Therefore, NPMs are useful test-beds to develop seizure control strategies.^[9]

The most commonly used treatment for epilepsy is administration of antiepileptic medication, usually anti-consultants. While frequently successful, such drugs can have significant systemic and central nervous system side effects. An alternative, on which we focus in this report, is the suppression of epileptiform activities through electrical brain stimulation, which seems to have minimal or no advserse side effects.^[10] Furthermore, epilepsy is drug-resistant in a substantial fraction of patients, and the affectivity of these drugs frequently diminishes over time even in patients who first responded favorably to anti-epileptic medication. Minimally invasive transcranial application of direct current pulses has been demonstrated to be an effective anti-consultant measure in a rodent model of epilepsy.^[11] This is the case even though the applied current is very simple and open loop (trains of rising bipolar pulses through a unilateral epicranial electrode).  The adaptive (closed-loop) control is likely to be a much more powerful stimulation control strategy and it is therefore the focus of current experimental and theoretical researches.^[12] An optimal control theory was employed to design an event-based, minimum-energy, desynchronizing control stimulus for a network of pathologically synchronized, heterogeneously coupled neurons.^[12] One possibility was to use a radial basis function neural network responsive neuromodulator to control seizure-like events in a computational model of epilepsy.^[13] Also, a closed-loop controller was developed based on a recursively identified autoregressive model of the relationship between stimulation input and local field potential output.^[14] Although excellent performance has been achieved by these approaches, these intelligent or optimal control algorithms are somewhat complex. This suggests that simple control algorithms should be preferred in closed-loop seizure control.

Closed-loop control is a promising method for suppressing epilepsy.^[15] A much used closed-loop control scheme is a proportional-integral-derivative (PID) controller, which gains its popularity in engineering because of its robust performance and simplicity.^[16] A PID type closed-loop control strategy was introduced to control epileptic seizures by disrupting the synchronization patterns observed prior to the seizures.^[9] However, the selection of the PID control parameters in the previous work was determined by trial and error. This makes the outcome strongly dependent on the designer's experience. Furthermore, the work is time-consuming and efficiency is not guaranteed, which makes the selection of control parameters quite challenging.

Given the non-stationary of neuronal signals over long time periods, the integral part of a PID controller likely contributes little functionality. Proportional-derivative (PD) control is a special case of PID, which provides a fast response to the control error and increases the system's stability.^[16] Therefore, in this study, PD control is employed to control the epileptiform activity in an NPM. Since the epileptiform activity is caused by unbalanced excitation and inhibition, ^[17] intuitively the control parameters should be related to, and dependent on, the excitatory and inhibitory parameters. How to build a quantitative relationship between the PD control parameters and the model's excitatory and inhibitory parameters is, however, not known. We therefore apply the control theory to establish the above-mentioned relationship in the form of stabilizing regions, which can provide a guideline for the choice of control parameters. Our strategy is to start with a neural population model that shows the epileptic behavior due to its (abnormal) excitatory and inhibitory parameters. The epileptiform activities generated by the model are then suppressed by subjecting it to stimulation applied by a closed-loop PD controller. With the stabilizing region of the PD controller available for a given set of model parameters, we can avoid the time-consuming stability check for each set of PD controller parameters in the search process and thereby save the controller tuning time.

This report is organized as follows. After this introduction section, the neural population model is introduced and its transfer function representation is derived in Section 2. The stability of the NPM is analyzed based on the roots locus method in Section 3. In Section 4, a PD type closed-loop control strategy is proposed to suppress epileptiform activities in the NPM, and the stabilizing region of the PD controller is determined using a graphical approach. Simulation examples are given to demonstrate the performance of the proposed control scheme in Section 5. Finally, concluding remarks are given in Section 6.

Epilepsy is a dynamical disease of the brain^[2] caused by lack of balance between excitation and inhibition.^[16] In this section, we consider the establishment and disturbances of neurophysiologtical homeostasis in an NPM to study the dynamical mechanisms underlining epileptiform activity in the NPM. This will be helpful to understand the mechanisms of seizure control strategies developed in the next section.

Homeostasis is a fundamental mechanism that stabilizes and maintains “ normal” activity levels in neural (and technical) systems.^[20] In the NPM, a homeostatic equilibrium is achieved by the balance between excitation and inhibition. Excessive parameter drift of the NPM (an abnormal increase of excitatory parameters H_e or an abnormal decrease of inhibitory parameters H_i) can cause unbalance between excitation and inhibition, resulting in the generation of epileptiform activity by destabilizing the NPM. Thus, stability analysis of the NPM will provide insight into the dynamical mechanisms that generate epileptiform activities.

According to the linear control theory, the stability of a control system can be characterized completely by the location of its characteristic roots.^[21] Thus, we plot the root loci of the characteristic equation of the NPM and use them to analyze the influence of the model parameters H_e and H_i on the stability of the NPM.

From Eq.  (5), the characteristic equation of the NPM is

To analyze the influence of the parameters H_e and H_i on the stability of the NPM, we vary them independently while keeping all other model parameters at their standard values, as shown in Table  1. Therefore, the characteristic roots of Eq.  (6) depend on H_e and H_i, i.e., s_i = f(H_e, H_i), where s_i is the i-th root of Eq.  (6), i = 1, 2, … , 6.

Fig.  3.

	Figure Option ViewDownloadNew Window
	Fig.  3. Root loci of the neural population model. Only parameters (a) He and (b) Hi are varied, the other model parameters stay at the standard values given in Table  1.

According to Eq.  (6), we plot the root loci of the characteristic equation of the NPM, which are illustrated in Figs.  3(a) and 3(b). The characteristic equation (6) is sixth order, so there are six blue branches. The model parameters are given in Table  1 except for H_e which is varied in Fig.  3(a) and H_i which is varied in Fig.  3(b). The variation range of each parameter is specified in the title of the corresponding subplot. Stars show the locations of the roots of the characteristic equation of the NPM, while diamonds and squares show the root locations when the corresponding parameter is set at its minimum and maximum values, respectively.

This procedure provides considerable insight into how the parameters influence the system's stability. Abnormal increase of excitation (i.e., hyperexcitation) as well as decrease of inhibition (i.e., low inhibition) results in roots located in the right half of the s plane, with positive real parts. These roots correspond to parameters H_e and H_i that lead to an unstable behavior of the NPM, which results in the generation of the epileptiform activities by the NPM.

As described in the previous section, the unbalanced excitation and inhibition makes the NPM unstable and leads to epileptiform activities with high amplitudes.

In this section, we consider the neural population model as a controlled plant from the viewpoint of control theory. We construct a PD controller to stabilize the unstable NPM such that homeostasis is maintained and epileptiform activities are suppressed. The basic requirement for the PD controller parameters is to make the closed-loop system stable. A graphical approach is convenient for determining the stabilizing region in the parameter space for tuning parameters of the PD controller.

4.1. Control scheme

The PD based NPM (PD-NPM) closed-loop control system is illustrated in Fig.  4, where G_NPM (s) is the transfer function of the NPM, and G_pd(s) is a PD controller with an ideal transfer function

Furthermore, u(t) is the control signal, p(t) the excitatory input of the NPM from neighboring areas, and y(t) is its output.

Fig.  4.

	Figure Option ViewDownloadNew Window
	Fig.  4. Schematic diagram of the PD-NPM control system.

4.2. Synthesis of stabilizing the PD controller

The objective is to determine the region of the PD controller in the (K_p, K_d) plane in which the PD-NPM control system is stable. We achieve this by using a graphical approach to determine the stability boundary of the PD controller in the (K_p, K_d) plane.

From Fig.  4, we see that the characteristic equation of the closed-loop PD-NPM system is

where Δ (s) is the characteristic polynomial of the PD-NPM control system.

According to the linear control theory, the stability of a control system can be determined according to the location of the characteristic roots. A feedback system is stable if all the characteristic roots lie in the open left-hand s plane of the imaginary axis, and unstable if any characteristic root lies in the open right-hand side.

Defining s = jω , and substituting it in Eq.  (8), we obtain

which determines the stability boundary of the PD-NPM control system.

As the coefficients of the characteristic polynomial (9) are real, the complex roots must be complex conjugate. Hence, it is sufficient to consider the case ω ∈ [0, ∞ ). The stability boundary is then given as follows:^[20]

where

In the following, we discuss the two cases, ω = 0 and ω > 0, to determine the boundary of the stabilizing region of the PD controller in the (K_p, K_d) plane.

4.2.1. Case 1: ω = 0

Substituting Eqs.  (5) and (7) in Eq.  (8), we can derive the characteristic equation of the PD-NPM control system as follows:

For ω = 0, i.e. s = 0, equation (11) becomes

Thus we can derive the stability boundary for the case, i.e., H₁, as follows:

This equation describes a line parallel to the axis through the point in the (K_p, K_d) plane, which is shown in Fig.  5 (red line).

Next, from Eq.  (11), we note that the coefficient of s⁶ in the characteristic polynomial, i.e., , is obviously positive. According to the Routh stability criterion, the constant term should be positive too, i.e.,

that is

Therefore, we obtain the following lemma.

Lemma 1 For the case ω = 0, along the direction of the arrow shown in Fig.  5, the region to the left of the line defined in Eq.  (12) in the (K_p, K_d) plane is the stabilizing parameter sets of the PD controller.

4.2.2. Case 2: ω > 0

From Eq.  (9), we obtain

where δ _{R_NPM} (ω ) and δ _{I_NPM} (ω ) are the real and the imaginary components of G_NPM (jω ), respectively. Partitioning the corresponding Δ (jω ) into real and imaginary components yields

where

Here δ _R (ω ) and δ _I (ω ) are the real and the imaginary components of Δ (jω ), respectively.

Note that both δ _R (ω ) and δ _I (ω ) are dependent on K_p and K_d. Thus, the stability of the characteristic equation  (9) can be investigated in the parameter space (K_p, K_d) as follows.

From Eqs.  (15) and (16), we obtain

Solving the above two equations simultaneously, we obtain

where is the complex modulus of G_NPM (jω ).

Equation  (17) defines the critical curve H₂ in the (K_p, K_d) plane for the case ω > 0. We sketch the stability boundary curve (K_p (ω ), K_d (ω )) for ω > 0 in Fig.  5 (blue curve).

We thus find that the stability boundary is composed of the stability boundary line defined in Eq.  (12) and the stability boundary curve defined in Eq.  (17). As illustrated in Fig.  5, the stability line and the stability boundary curve divide the (K_p, K_d) plane into three disjoint regions, denoted by R₁, R₂, and R₃, where R₂ is the stabilizing parameter regions, and R₂ and R₃ are the destabilizing parameter regions. In the following, we will investigate how to characterize the system behavior in these regions.

Fig.  5.

	Figure Option ViewDownloadNew Window
	Fig.  5. The stability boundary of the PD controller for suppressing epileptiform activities in PD-NPM with He = 4.5  mV. Arrows represent the direction of the curve (Kp (ω ), Kd (ω )) in which ω increases.

To begin with, we introduce the following proposition.^[22] Proceeding along the curve in the direction of increasing ω , the right side of the curve (K_p (ω ), K_d (ω )) is the stabilizing parameter region, where det J < 0; while to its left, det J > 0. Here J is the Jacobin matrix defined as

From Eq.  (17), we obtain

Thus we conclude that det J = ω | G_NPM (jω )| ². Since ω > 0 and| G_NPM (jω ) | ² > 0, it follows that

Therefore, we obtain the following lemma.

Lemma 2 For the case ω > 0, following the curve (K_p (ω ), K_d (ω )) in the direction of increasing ω , the parameter space to the left of the curve defined by Eq.  (17) is the stabilizing parameter region of the PD controller.

Combining Lemma 1 with Lemma 2 yields the following theorem.

Theorem 1 For the PD-NPM control system depicted in Fig.  4, following the curve (K_p (ω ), K_d (ω )) defined in Eq.  (17) in the direction of increasing ω (ω > 0), the portion of the parameter space to the left of the curve defined by Eq.  (17) and the line defined by Eq.  (12) is the exact stabilizing parameter region of the PD controller.

Thus, according to Theorem 1, the region denoted by R₁ in Fig.  5 can be identified as the stabilizing region of parameters (K_p, K_d).

4.3. Influence of model parameters on the stabilizing region of the PD controller

The epileptiform activities in the NPM with abnormal values of excitatory and inhibitory parameters are the result of unbalanced excitation and inhibition. Therefore, in this subsection, we will analyze how the parameters H_e and H_i influence the stabilizing region of the PD controller.

Fig.  6.

	Figure Option ViewDownloadNew Window
	Fig.  6. The influence of excitatory and inhibitory parameters He and Hi on the stabilizing regions of the PD controller. In panel (a), He is varied as shown and Hi is fixed at 22  mV; in panel (b), Hi is varied and He is fixed at 3.25  mV.

We consider two cases with abnormal values of H_e and H_i . The stabilizing regions of the PD controller for the NPM with different abnormal values of the parameters H_e and H_i are illustrated in Figs.  6(a) and6(b), respectively. The results indicate that increasing H_e or decreasing H_i makes the stabilizing region of the PD controller smaller. This agrees with the intuitive notion that a stronger control signal is needed to stabilize the NPM if the latter has stronger excitation and/or weaker inhibition.

We conduct simulations to demonstrate the feasibility and efficiency of the proposed PD based control scheme. We use the NPM with the standard parameters from Table 1 but with the following two variations. In the first, we increase the excitation to H_e = 6.5  mV. In the second, we increase the excitation to a moderate one (H_e = 3.95  mV) and also decrease the inhibition to H_i = 18  mV. As shown in Fig.  3, in each of these cases, the NPM is linear unstable. The simulation results of the model dynamics are illustrated for these cases in Figs.  7(a) and 8(a), respectively. The output function y(t) of the NPM, which is in agreement with the previous work, ^[8] is interpreted as an approximation of the EEG signal. Without control applied (Figs.  7(a) and 8(a)), the EEG signal consists of spike-like high amplitude oscillations with obvious epileptiform characteristics.

Fig.  7.

	Figure Option ViewDownloadNew Window
	Fig.  7. The output of NPM without (a) and with (b) the PD controller for the first case (He = 6.5). Note the difference in scales.

In order to suppress the seizure-like behavior, we apply a PD controller as described with control parameters inside the stabilizing region illustrated in Fig.  6. Specifically, the control parameters are selected as K_p = 230, K_d = 0.2 for the first case (Fig.  7) and K_p = 120, K_d = 0.2 for the second one (Fig.  8). The simulation results of the controlled NPM are shown in Figs.  7(b) and 8(b). The output is now a low amplitude, noise-like activity similar to that seen during the normal (non-epileptic) behavior.

Fig.  8.

	Figure Option ViewDownloadNew Window
	Fig.  8. The output of NPM without (a) and with (b) the PD controller for the second case (Hi = 18, He = 3.95). Note the difference in scales.

We use a commonly used neural population model to study the generation of epileptic behavior. Parameters of the model are chosen such that seizure-like activities are obtained. We then consider this system from the point of control theory as one whose behavior needs to be steered towards a normal (non-pathological) state. Our goal is to use a simple, proportional-derivative (PD) control scheme, both for robustness and ease of parameter determination. We show that the parameters of the PD controller can be determined from graphically-defined stability regions which are defined in terms of the model parameters. This establishes a quantitative relationship between the control parameters and the model's excitatory and inhibitory parameters, thus providing a guideline for the choice of the control parameters to suppress epilepsy in the model. Simulation results suggest the feasibility of the PD based control scheme. A direct current stimulation may be effective for neurological disorders other than epilepsy, ^[22] and the closed loop strategies developed here may prove to be much more powerful than the typically very simple stimulation protocols that were conventionally applied. The linear approximation of the sigmoid function in our model will generate more unstable output than the full (nonlinear) model. Thus, determining the stability region of the PD controller based on the stability analysis of the linearized model provides a conservative estimate. In future work, we will use a bifurcation analysis method to obtain more specific results. An open question is how to determine the control parameters from measurements in a biological system, i.e., when no computational model with known parameters is available. We also propose to apply the proposed PD controller to standard animal models of epilepsy with the goal of on-line seizure suppression, with the obvious long-term goal of application to human patients.