Closed-loop controller based on reference signal tracking for absence seizures

Absent epilepsy is a kind of refractory epilepsy, which is characterized by 2–4 Hz spike and wave discharges (SWDs) in electroencephalogram. Open-loop deep brain stimulation (DBS) targeting the thalamic reticular nucleus (TRN) is an effective method to treat absent epilepsy by eliminating SWDs in the brain. Compared with open-loop DBS, closed-loop DBS has been recognized by researchers for its advantages of significantly inhibiting seizures and having fewer side effects. Since traditional trial-and-error methods for adjusting closed-loop controller parameters are too dependent on the experience of doctors, in this paper we designed two proportional integral (PI) controllers based on the basal ganglia-cortical-thalamic model, whose PI parameters are calculated from the stability of the system. The two PI controllers can automatically adjust the frequency and amplitude of DBS respectively according to the change of the firing rate detected by substantia nigra pars reticulata (SNr). The parameters of the PI controller are calculated based on the Routh-Hurwitz stability criterion of a linear system which transformed by the original system using controlled auto-regressive (CAR) model and recursive least squares (RLS) method. Numerical simulation results show that both PI controllers significantly destroy the SWDs of the cerebral cortex and restore it to the other two normal discharge modes according to the different target firing rate, which supplies a promising brain stimulation strategy.

Absence epilepsy has been widely studied for its unique pathogenesis 1,2 , including the main frequency of 2-4 Hz spike and wave discharges (SWDs) in electroencephalogram (EEG) 3,4 , sudden and temporary loss of consciousness 5 , and widespread seizures in the adolescent years 6,7 . Previous studies have suggested that the thalamus plays a crucial role in epileptic seizures and that abnormal feedback from the thalamus to the cerebral cortex is responsible for the appearance of SWDs 8,9 . Subsequently, electrophysiological records in the cerebral cortex and thalamus of epileptic patients also demonstrated the corticothalamic system was associated with SWDs 10,11 , and the corresponding corticothalamic model was established 12,13 . The basal ganglia, as the information processing unit of the brain, directly or indirectly participates in the information transmission between the cortex and the thalamus 14,15 . Clinical studies have shown that damage to basal ganglia can lead to various brain disorders, including epilepsy 16 , cognitive impairment, and Parkinson's disease 17,18 . The basal ganglia are composed of subthalamic nucleus (STN), globus pallidus internal (GPi), globus pallidus external (GPe), and substantia nigra pars reticulata (SNr), of which the most important output nucleus is SNr neurons. Moreover, experimental studies on SNr in rodent models suggested that SNr plays an important role in the control of absence epilepsy 19,20 . Therefore, given that underlying mechanism of SNr in the brain is not clear, further research on SNr is necessary.
Slow oscillations of 0.1-0.2 Hz in the human brain are traveling waves that can trigger the thalamic spindles and periodically sweep across the cerebral cortex 21 . Spindles, which are associated with slow oscillations, can be observed in the electroencephalogram (EEG) in deep sleep 22 . Furthermore, experimental records showed that sleep spindles disappear when the transmission between thalamus reticular nucleus (TRN) and specific relay nucleus (SRN) is disrupted 23,24 . Therefore, researchers believe that TRN can act as a pacemaker of the spindles, and that the slow oscillation of different brain regions can also be regulated by TRN 25,26 . In recent years, more and more biophysical models have been established to study TRN. Fan et al. achieved the interconversion between spindles and SWDs in the cerebral cortex by applying deep brain stimulation (DBS) stimulation to TRN 27 . Wang et al. investigated the mechanism of external stimulation on TRN in the cortical thalamic model 28 . These studies have stimulated our enthusiasm to study the underlying mechanism of DBS stimulation on TRN.
There are many ways to treat epilepsy and the traditional methods are drug therapy or surgical removal 29,30 . Although these methods have a certain control effect on the absence epilepsy, the risks of surgery and the side  55,56 . Based on this model, we can successfully achieve the transformation of spike and wave discharges (SWDs) to spindle oscillations in the cerebral cortex. As shown in Fig. 1, the whole network model can be divided into three categories according to color, namely cerebral cortex (orange), thalamus (blue) and basal ganglia (green), each of which can be divided into different sub-modules. For example, the cerebral cortex is made up of two orange rectangles, called inhibitory interneurons (II) and excitatory pyramidal neurons (EPN). The different submodules are connected by projections: excitatory projections mediated by glutamate (purple arrow solid line), inhibitory projections mediated by GABA A (blue arrow solid line), and inhibitory projections mediated by GABA B (blue arrow dot line).
The biophysical model of absence epilepsy is shown as follows 55 : (1) where α is synaptodendritic decay time constant, and β is synaptodendritic rise time constant. γ e denotes cortical damping rate, and i ∈ [EPN, II, TRN, SRN, D 1 , D 2 , SNr, GPe, STN] represents one of nine groups of neurons. r indicates the spatial position, σ represents the threshold variability of firing rate, θ i represents the mean firing threshold, and Q max i indicates the maximum firing rate. φ e denotes cortical excitatory axonal field. τ and φ n are represent the GABA B delay and the constant nonspecific subthalamic input onto SRN. v i,j indicates coupling strength from i to j,i, j ∈ [EPN, II, TRN, SRN, D 1 , D 2 , SNr, GPe, STN] . The other parameters are shown in Table 1 55 .
Establish stimulus-response relation. The design diagram of the closed-loop DBS controller is illustrated in Fig. 2A. We selected the firing rate of SNr as the reference signal r(m) . Then, we input the difference between r(m) and expected mean firing rate r rs (m) into the proportional integral (PI) controller to obtain the DBS stimulus parameter s(m) . Finally, we obtained DBS stimulation for the control of absence epilepsy. Since absence epilepsy is highly nonlinear, direct use of PI controller is not appropriate. Figure 2B shows the linear system of a closed-loop DBS controller, which is used to describe the linear relationship between DBS stimulus parameter and the mean firing rate of SNr. The CAR model is calculated as follows 57 : where z, r(m) , and s(m) are the lag operator, output signal, and input signal, respectively. n a ,ε(m) , and n b are the order of output, error, and the order of input, respectively.
As we know, the CAR model describes a linear relationship between stimulus and responses and relies heavily on input and output data. Therefore, to extract more information from the relationship between DBS parameters and mean firing rate of SNr, we let DBS parameters change randomly in a certain range. For example, Fig. 3A 1 + a 1 z −1 + a 2 z −2 + · · · + a n a z −n a r(m) = b 0 + b 1 z −1 + b 2 z −2 + · · · b n b z −n b s(m)+ε(m), www.nature.com/scientificreports/ respectively show the discharge of DBS current during 4-6 s, and we change the frequency and amplitude of current respectively every 0.2 s. In addition, we also consider the influence of time window bin on the control accuracy of closed-loop DBS. As shown in Fig. 3C, we plotted the box diagram of mean firing rate of SNr as a function of the time window bin. The simulation results show that when the time window bin is 0.2 s, the mean firing rate of SNr measured each time is scattered, which means that the time window bin of 0.2 s is sensitive to the changes of the system. The recursive least squares (RLS) method is used to estimate CAR model parameters a 1 , a 2 , . . . , a n a and b 1 , b 2 , . . . , b n b . The least square form of the CAR model is shown as follows 57 : T is composed of past input and output data and current input data. θ = [a 1 , a 2 , · · · , a na , b 0 , b 1 , · · · , b nb ] T can be estimated by the following equation: Figure 4 shows the performance of CAR model tracking dynamically changing firing rate of SNr when n a = 2 and n b = 2 . We periodically change the DBS parameters from 1 s and start tracking the firing rate of SNr using the CAR model. When the input of CAR model is the frequency and amplitude of DBS current respectively, the corresponding performance of CAR model in model training is illustrated in Fig. 4A1,B1 respectively. Figure 4A2,B2 show the parameter estimation process of CAR model and the fitting of test data is shown in Fig. 4A3,B3 respectively. As we can see in Fig. 4, the CAR model performs well in training and testing, even  www.nature.com/scientificreports/ though the order of the CAR model is relatively low. We should be more interested in the stability than the fitting accuracy of CAR model. Therefore, the CAR model in this paper is expressed as follows: Design of closed-loop controller. The system replacing the absence epilepsy model with CAR model is shown in Fig. 2B and the structure of PI controller is as follows: where k p and k i are computed by the Routh-Hurwitz stability criterion, which is introduced by the following equations.
The forward transfer function of the system in Fig. 2B is as follows: The closed-loop transfer function is The characteristic equation of this system is where: Then, multiplying both sides of equation by (w − 1) 3 , we get where: Based on the Routh-Hurwitz stability criterion, the stability of this system is equivalent to n i k p , k i > 0(i = 0, 1, 2, 3) , and the corresponding subset of feasible solutions is k p and k i > 0.
Traditional closed-loop DBS designs mostly control the amplitude of DBS to adjust the intensity of DBS 58 . However, from the perspective of signal analysis, frequency modulation has stronger anti-interference ability than amplitude modulation. To sum up, we designed two types of closed-loop DBS controllers, one of which is based on frequency modulation (BoFM). We limit the frequency to 20-100 Hz, and the BoFM controller is shown as follows: The other is a closed-loop DBS controller based on amplitude modulation (BoAM). We limit the amplitude to 1-4 mA, and the BoAM controller is shown as follows: .
(17) (w − 1) 3 D(w) = n 3 w 3 + n 2 w 2 + n 1 w + n 0 = 0, (MathWorks, USA) simulation environment. The differential equation was solved by the standard fourth order Runge-Kutta method, with the temporal resolution of numerical integration is 0.05 ms. In the process of parameter training of CAR model, the optimal parameters can be obtained by running for sufficient time (> 15 s). The time series φ e of cortical excitatory axon field was used to describe the cortical macro dynamics, the transient waveform which was unstable at the beginning was discarded to ensure that the time series used for analysis was obtained under the premise that the system was stable. The fast Fourier transform is used to estimate the dominant frequency of time series φ e .

Dynamic change of firing rate in substantia nigra reticulum induced by open-loop deep brain stimulation.
Our model successfully demonstrated the dynamic transition from spike and wave discharges  Correspondingly, the main frequency of cortical electrical activity also decreases to 0 Hz in Fig. 5D. As we can see, Fig. 6 clearly shows the four different firing states of the cerebral cortex and the corresponding firing rate of the substantia nigra pars reticulata (SNr). Figure 6A2-D2 respectively show the time series of cerebral cortex discharges in saturation state, SWDs state, simple oscillation state and low firing state, and the corresponding SNr firing rates are shown in Fig. 6A1-D1. Compared with the more stable firing rates such as saturation state and low firing state, the mean firing rates of SWDs oscillating state and simple oscillating state are constantly changing, which brings us difficulties in designing the closed-loop controller with mean firing rate as reference signal in the next section. For simplicity, we calculated the mean firing rate of four firing activity and find that the mean firing rate of SNr is arranged from high to low in the order of saturation state, SWDs oscillation state, simple oscillation state, and low firing state. Therefore, we chose the average value of firing rates of simple oscillation and low oscillation respectively as the expected firing rates of the closed-loop DBS controller, which are r rs = 27.7884 Hz and r rs = 28.4144 Hz, respectively.
It has been known from Fig. 5C that DBS stimulation current can inhibit the firing activity of cerebral cortex, therefore we also want to research whether DBS stimulation current can also control the firing rate of SNr. For this purpose, we investigated the effects of frequency and amplitude of DBS stimulation on MFR of SNr. As shown in Fig. 7A1,A2,B1,B2, when the frequency of DBS stimulation increases from 20 to 40 Hz, the cortical firing activity changes from SWDs oscillation state to simple oscillation state, and the amplitude of the firing rate oscillation of SNr decreases. At the same time, we found a similar situation when the amplitude of DBS stimulus increased. As shown in Fig. 8A1,B1,A2,B2, as the amplitude of DBS stimulation increases from 1 to 3 mA, the amplitude of the firing rate of SNr decreases significantly, and the cortical firing activity changes from SWDs oscillation state to simple oscillation state. In addition, an interesting phenomenon is that we detected simple oscillations of 60 Hz with very low amplitude in the low firing state suppressed by DBS stimulation in Figs. 7C1,C2, 8C1,C2. This may be due to the introduction of artificial high-frequency electrical stimulation in the absence epilepsy model.
To test whether the above results can be generalized to a larger parameter range, we extended the frequency and amplitude ranges to [20,100] Hz and [1,4] [1,4] mA respectively. As shown in Fig. 9A1, the mean firing rate of SNr increases to a peak at Freq2 and then decreases slowly with increasing frequency. Combined with Fig. 9A1,A2, it is found that the mean firing rate of SNr when the firing activity is in simple oscillation state is higher than that when the firing activity is in SWDs oscillation state. The same situation can be observed in Fig. 9B1, B2. Therefore, it can be concluded that when DBS stimulation intensity increases, including the  www.nature.com/scientificreports/ frequency or amplitude of DBS stimulation, the mean firing rate of SNR increases to a peak at the transition point from simple oscillation state to low discharge state, and then slowly decreases.
Efficacy of closed-loop control strategy based on frequency modulation. After the analysis in the previous section, we known that the increase of DBS stimulation intensity will decrease the amplitude of SNr firing rate, which constitutes the theoretical basis for our design of closed-loop DBS controller. In this section, we will focus on the control effect of absence epilepsy under the BoFM strategy. As shown in Fig. 10, we designed a closed-loop controller based on BoFM control strategy for the purpose of controlling the state of cortical discharge activity. Meanwhile, considering the influence of expected firing rate on DBS control results, we selected the mean firing rate of SNr under the condition of cortical low firing state as the expected firing   Fig. 10A,C, respectively. Moreover, we also selected the mean firing rate of SNr under the condition of cortical simple oscillation state as the expected firing rate to study the control effects of cortical activity. As shown in Fig. 11B, the cortical firing activity realized the transition from SWDs oscillation state to simple oscillation state at T1 = 1.1772s. The variations of SNr firing rate and closed-loop DBS current are shown in Fig. 11A,C, respectively. Compared with Fig. 10B, the cortical firing activity in Fig. 11B did not switch to a low firing state but presented a simple oscillation state. The reason may be that we increase the expected firing rate of the closed-loop controller. Compared with Fig. 11C and Fig. 10C, The DBS current in Fig. 11C is sparser than that in Fig. 10C, which means that the power consumption of the closed-loop DBS current is inversely proportional to the value of the expected firing rate. In addition, an interesting finding is that the response time (T1 = 1.1772s) of simple oscillation in Fig. 11B is smaller than the response time (T1 = 1.20825 s) of simple oscillation in Fig. 10B.
Efficacy of closed-loop control strategy based on amplitude modulation. In the previous section, we successfully implemented the closed-loop control under the BoFM strategy, and we found that under the BoFM strategy, the control effect induced by different expected rate of fire is different. Considering the influence of DBS amplitude on SNr firing rate, we will focus on the control effect of absence epilepsy under the BoAM strategy in this section.
Like Fig. 10B, when the expected firing rate is 27.7884 Hz, we realize the transition from SWDs oscillation state to low firing state under the BoAM strategy. As shown in Fig. 12B, the simple oscillation state appeared at time T1 = 1.19885 s, and the low firing state appeared at time T2 = 4.8032 s. The variations of SNr firing rate and closed-loop DBS current are shown in Fig. 12A,C, respectively. However, the response time T2 of low discharge state under BoAM strategy is much larger than that under BoFM strategy. Moreover, the amplitude of the closed-loop control BDS current under the BoAM strategy rises from 1 mA to nearly 4 mA and gradually remains stable in Fig. 12C. We also analyze the control of cortical firing activity when the expected firing rate is 28.4144 Hz, which is the mean firing rate of SNr when the cortex is in simple oscillation state. As shown in Fig. 13B, when the expected firing rate is 28.4144 Hz, the cortical firing activity changes from SWDs oscillation state to simple oscillation state at time T1 = 1.19885 s. The variations of SNr firing rate and closed-loop DBS current are shown in Fig. 13A,C, respectively. The response times of simple oscillations are similar in Figs. 12B and 13B, which is different from we concluded in the previous section. In addition, we find that the amplitude of closed-loop DBS current in Fig. 13C is chaotic, which is different from Fig. 12C. However, like the conclusion in the previous section, the DBS current in Fig. 13C is sparser than that in Fig. 12C, which means that the DBS current in Fig. 13C consumes less power than that in Fig. 12C.

Conclusion and discussion
In this paper, we successfully achieved the conversion of spike and wave discharges (SWDs) oscillating state to other firing state in the cerebral cortex using a closed-loop deep brain stimulation (DBS) controller based on a model of absence epilepsy involving cerebral cortex, thalamus, and basal ganglia. The mean firing rate of substantia nigra pars reticulata (SNr) is selected as a reference signal, and the cortical firing activity is selected as a biomarker reflecting the epileptic state. The controlled auto-regressive (CAR) model and the Routh-Hurwitz According to the numerical results, the intensity of the open-loop DBS current, including the amplitude or frequency of the DBS current, will reduce the amplitude of the SNr firing rate. With the increase of current intensity, the mean firing rate of SNr increase firstly and then decrease, and the firing activity of cerebral cortex change from SWDs oscillation state to simple oscillation state and finally kick to low firing state. Therefore, we designed closed-loop controllers based on BoFM strategy and BoAM strategy, respectively. The closed-loop controllers with two control strategies can achieve the desired control effect under different expected firing rates. By adjusting the expected firing rate from 27.7884 Hz to 28.4144 Hz, the final firing state under the closed-loop control is transformed from low firing state to simple oscillation state. Meanwhile, the power consumption of closed-loop DBS current decreases with the increase of expected firing rate. In addition, an interesting phenomenon is that the response time of low firing state under BoAM strategy is much larger than that under BoFM strategy when the expected firing rate is 27.7884 Hz, we hope that the research can provide reference and help for the treatment and prevention of epilepsy patients.  www.nature.com/scientificreports/ The PI parameters of the closed-loop DBS controller proposed in this paper are calculated by using the stability of the system rather than by traditional trial-and-error adjustment 44 . Therefore, the challenge to apply the proposed closed-loop DBS method to clinical or experimental applications is to identify the relationship between stimulus intensity and reference signal. In our study, the discharge rate of SNr was selected as the reference signal, and the frequency and amplitude of DBS were selected as the stimulus intensity. Although the proposed closed-loop control algorithm can eliminate the SWDs oscillation in the cerebral cortex, some limitations of the study cannot be ignored. First, the PI controller can track the mean firing rate of SNr well only when the change of frequency and amplitude is less than 1 Hz. In other words, a faster change in the target mean firing rate than 1 Hz will lead to an increase in tracking error, which may be caused by unmodeled dynamics. Therefore, to improve the tracking performance of the dynamic reference signal, an adaptive controller is designed to adjust the parameters of PI controller to adapt to the dynamic change of mean firing rate. Second, although TRN has been shown to be effective as a stimulus target for the treatment of absent epilepsy 27,28 , the centromedian thalamic nucleus 59 , subthalamic nucleus 60,61 , and anterior nucleus 62,63 may also be potential stimulation targets as they are associated with other types of refractory epilepsy. Finally, the basal ganglia-cortical-thalamic network is highly non-linear, so it is necessary to use a non-linear controller in the following work.

Data availability
The data used to support the findings of this study are included within the article.