Suppression of seizure in childhood absence epilepsy using robust control of deep brain stimulation: a simulation study

Deep brain stimulation (DBS) is a promising technique to relieve the symptoms in patients with intractable seizures. Although the DBS therapy for seizure suppression dates back more than 40 years, determining stimulation parameters is a significant challenge to the success of this technique. One solution to this challenge with application in a real DBS system is to design a closed-loop control system to regulate the stimulation intensity using computational models of epilepsy automatically. The main goal of the current study is to develop a robust control technique based on adaptive fuzzy terminal sliding mode control (AFTSMC) for eliminating the oscillatory spiking behavior in childhood absence epilepsy (CAE) dynamical model consisting of cortical, thalamic relay, and reticular nuclei neurons. To this end, the membrane voltage dynamics of the three coupled neurons are considered as a three-input three-output nonlinear state delay system. A fuzzy logic system is developed to estimate the unknown nonlinear dynamics of the current and delayed states of the model embedded in the control input. Chattering-free control input (continuous DBS pulses) without any singularity problem is the superiority of the proposed control method. To guarantee the bounded stability of the closed-loop system in a finite time, the upper bounds of the external disturbance and minimum estimation errors are updated online with adaptive laws without any offline tuning phase. Simulation results are provided to show the robustness of AFTSMC in the presence of uncertainty and external disturbances.

Suppression of seizure in childhood absence epilepsy using robust control of deep brain stimulation: a simulation study Ehsan Rouhani 1* , Ehsan Jafari 2 & Amir Akhavan 1 Deep brain stimulation (DBS) is a promising technique to relieve the symptoms in patients with intractable seizures. Although the DBS therapy for seizure suppression dates back more than 40 years, determining stimulation parameters is a significant challenge to the success of this technique. One solution to this challenge with application in a real DBS system is to design a closed-loop control system to regulate the stimulation intensity using computational models of epilepsy automatically. The main goal of the current study is to develop a robust control technique based on adaptive fuzzy terminal sliding mode control (AFTSMC) for eliminating the oscillatory spiking behavior in childhood absence epilepsy (CAE) dynamical model consisting of cortical, thalamic relay, and reticular nuclei neurons. To this end, the membrane voltage dynamics of the three coupled neurons are considered as a three-input three-output nonlinear state delay system. A fuzzy logic system is developed to estimate the unknown nonlinear dynamics of the current and delayed states of the model embedded in the control input. Chattering-free control input (continuous DBS pulses) without any singularity problem is the superiority of the proposed control method. To guarantee the bounded stability of the closedloop system in a finite time, the upper bounds of the external disturbance and minimum estimation errors are updated online with adaptive laws without any offline tuning phase. Simulation results are provided to show the robustness of AFTSMC in the presence of uncertainty and external disturbances.
Epilepsy refers to a set of chronic neurological diseases typified by repetitive (at least two) spontaneous seizure attacks due to abnormal bursts of electrical activity in the brain that can bring about temporary brain dysfunction 1 . According to the regional emergence of excessive electrical activity in the brain, the International League Against Epilepsy (ILAE) classified seizures as focal (partial), generalized, and unknown 2,3 . Focal seizures are established in a local area of the brain. They can be further classified based on the patient's awareness during a seizure attack as aware (simple) and impaired awareness (complex). In generalized seizures, both hemispheres are involved at the onset of the seizure, and awareness is impaired in most cases. The regional onset is unspecified in unknown seizures, but other manifestations are known 3 . These classes can also be subdivided into the motor or non-motor subsets regarding movement manifestations [2][3][4] . Childhood absence epilepsy (CAE) is a subset of genetic/idiopathic generalized epilepsy characterized by frequent daily absence (non-motor) seizures in children between 2 to 13 years old. The typical features of CAE include sudden cessation of activity and awareness (e.g., staring blankly), short time duration episodes (approximately 9 s), and 2.5-3.5 Hz bilateral, synchronous spikeand-wave discharges (SWDs) in the recorded electroencephalogram (EEG) signal during seizure episodes in the brain 5,6 . Elucidating the pathogenesis of CAE, as well as developing efficacious remedies, are the fundamental challenges that should be tackled by researchers.
Various studies on patients with absence seizures and animal models acknowledged that the SWDs in electrophysiological recordings mainly originated from the aberrant neuronal interplay between the cerebral cortex and thalamus that constitutes the corticothalamic network. In fact, the reticular thalamic nucleus, thalamic relay nuclei, and cortical neurons form a reciprocal neuronal circuit demonstrating a significant contribution in the absence of seizure onset [7][8][9][10][11][12][13] . In accordance with experimental results, numerous computational corticothalamic models have been developed to discover the underlying mechanisms of transition from regular brain activity to SWDs [14][15][16][17][18][19][20][21][22][23][24][25] . To replicate the spontaneous SWDs in EEG recordings of rats with the absence of epilepsy, Suffczynski dynamics and parameter uncertainties due to the highly nonlinear and time-varying nature of the neurological systems, external disturbance rejection, and stable function of the controller in real applications have not been entirely addressed.
Apart from simulation studies, a handful of clinical DBS systems to suppress epilepsy have been developed in the last two decades. Responsive Neurostimulation System (RNS, NeuroPace) is the first FDA-approved commercially available closed-loop device that provides responsive brain stimulation with respect to the electrocorticogram (ECoG) recordings through one or two depth and/or subdural cortical strip leads that are implanted at the seizure focus of focal onset seizure patients. This device demonstrated favorable outcomes on 11 out of 27 patients with epilepsy [65][66][67] . In 68 , an implantable concurrent sensing and stimulation DBS device that relies on LFP and ECoG sensors was developed and validated in an ovine model of epilepsy. A support vector machine (SVM) classification algorithm was recruited to separate the biomarkers from the stimulation artifact. DiLorenzo et al. 69 devised a seizure advisory system in connection with an implantable 16-channels intracranial monitoring to predict seizures. The preclinical tests were performed in a canine model and then implantation of 15 patients to evaluate the method.
Sliding mode control (SMC) is a powerful nonlinear robust control strategy featuring outstanding performance to compensate for the parameter uncertainty and bounded external disturbances in highly nonlinear systems 70,71 . To achieve seizure abatement, a classical SMC was integrated with an RBFNN approximator to propel the EEG waves recorded from cortical areas with undesired dynamics to track the normal background waves 49 . In this work, the DBS-corticothalamic system was formulated to a canonical structure, and the system's dynamics were approximated by RBFNN, while disturbance rejection and model uncertainties were handled by the classical SMC. The main disadvantage of the proposed controller is the asymptotic stability of the SMC due to the linear structure of the sliding surface. To tackle the stability problem, terminal sliding mode (TSM) control that contains fractional order of sliding surfaces was developed to ensure the finite-time convergence of the system trajectories to the origin 70,72 . However, TSM cannot address the chattering phenomena 70 of the control input, which endangers the system's safety in real applications. Moreover, it suffers from the singularity problem that may lead to infinite control input to guarantee the ideal TSM motion [73][74][75] . Qian et al. 51 presented a finite-time fractional order SMC in combination with RBFNN to suppress epilepsy seizures in a well-established thalamocortical NMM. To design and stability analysis of the closed-loop system, an appropriate coordinate transformation was introduced to represent a regular form of the NMM. The main bottleneck of the work is that to transform the original NMM to the regular form and design the controller signal, prior knowledge of the NMM dynamics is necessary. To overcome the aforementioned problems (asymptotic stability, chattering, and singularity of the controller and the knowledge of the plant dynamics), in this paper, a robust control technique based on adaptive fuzzy TSM control (AFTSMC) is developed for control of the oscillatory spiking pattern of the three-neuron pattern of CAE whereby the controller is not dependent on any knowledge of the dynamics of the system to be controlled. The main contributions of this work are as follows: • To implement the proposed MIMO control algorithm, the dynamics of the membrane voltage of three-neuron CAE motif 8,25 should be represented and modeled as a three-input three-output nonlinear state delay model, and three continuous DBS pulses (control input) are designed to force the CAE pattern to track the normal state. This representation is based on the decomposition of the nonlinear dynamics of the current and delayed states of the CAE system 8,25 . Indeed, the continuous DBS pulses are automatically generated based on the MIMO model of the CAE motif. To achieve this, the unknown nonlinear dynamics of the MIMO model are approximated with the fuzzy logic system (FLS). • A fast TSM-type reaching law is embedded in the continuous nonsingular control input to accelerate the speed of the motion to the sliding surface far away from the surface in finite time. The superiority of the continuous control law is that it's chattering-free. • On the sliding surface, the fractional-integral structure of the surface guaranteed the convergence of the tracking error to the neighborhood of zero in finite time. • The nonlinear dynamics of the current and delayed states of the model embedded in the control input are approximated with the FLS, and adaptive laws based on the terminal gradient descent algorithm are developed for online updating the weights of the estimator. • To guarantee the bounded stability of the closed-loop system, the upper bounds of the external disturbance and minimum estimation errors are updated online using terminal-based adaptive laws without any offline tuning phase.

Model and control problem
Model. CAE is a 2.5-3.5 Hz bilateral, synchronous sudden spike followed by the wave discharges (SWDs) observed in the EEG pattern of the children during seizure episodes. SWD patterns with similar characteristics in human subjects were also reported in the experimental works conducted on animal homologs. Electrophysiological recordings of the spike-and-wave seizures in animal models showed the aberrant neuronal interplay between the cerebral cortex and the two main thalamic cells involved thalamocortical and the inhibitory neurons of the thalamic reticular nucleus. These experiments suggested that the connections of the GABA A and GABA B receptors can play an important role in the activation of spike-and-wave seizures. It has been proven that the seizure appeared in the mouse homolog by reduction of the conductance of GABA A receptors. Furthermore, both the thalamus and cortex are essential to the absence of seizures. As each cell type of the thalamus and cortex neurons contain the voltage-dependent currents necessary to define their properties, the dynamics of the electrical activity of the cortical (CT), thalamic relay (TC), and reticular (RT) nuclei neurons are represented as a conductance-based Hodgkin-Huxley type neuron model. www.nature.com/scientificreports/ in Fig. 1, and the dynamics of the membrane voltage of the neurons are modeled with the following nonlinear first-order differential equations 8,25 : where in (1), C TC = C CT = C RT = 1 μF cm −2 is the capacitance of the membrane of CT, TC, and RT neurons and V TC , V CT , and V RT is the membrane voltage of the CT, TC, and RT neurons, respectively. The ionic membrane currents (μA cm −2 ) are denoted by I L , a leak current, I Na , a transient voltage-gated current of Na + ions, I K , a transient voltage-gated current of K + ions, I TS , a low-threshold current of Ca ++ ions for RT neuron, I M , a depolarization-activated current of K + ions, I T , a low-threshold current of Ca ++ ions for TC neuron, I h , a mixed Na + -K + current activated by hyperpolarization, I K2 , a slow current of K + ions, and, I TC ext , I CT ext , and I RT ext , external currents of TC, CT, and RT neurons, respectively due to the input from other brain regions. Without the coupling effect of synaptic currents, if the external current of the neurons is set at zero, the voltage of the three neurons is at rest, and for a large nonzero value of the external currents, the voltage of the neuron is denoted with periodic spiking activity. In the simulation, the parameters of the CT and RT external current are set at 0 μA cm −2 and I TC is assumed to be 5 or 6 μA cm −2 . Detailed dynamics of the ionic membrane currents are provided in Appendix A of the Supplementary Materials.
The synaptic currents ( I i syn , i ∈ {TC, CT, RT} ) are in the following form: where the general form of AMPA, GABA A , and GABA B currents with the value of their parameters are summarized in Table 1. In the current equations, the parameters E (mv) and g (mS cm −2 ) are the reversal potential, and maximal conductance, respectively, and s is the fraction of the channels opened with the presynaptic neuron. In the dynamic equation of s, the parameter T is the concentration of neurotransmitters released with the presynaptic neuron, and the parameters α (m/s m/M) and β (m/s) are two constants. Two time delays τ 1 and τ 2 are included in the model of neurotransmitter concentration as the conduction delay for the corticothalamic and thalamocortical delays, respectively. In the simulation, the thalamocortical delay is set to 2.8 ms and the corticothalamic delay is varied between 2 and 10 ms based on the dynamics of normal and CAE states. In the model of the synaptic current from the RT neuron to the TC neuron due to the (GABA B ), the parameter r GABA B is the model of the receptor of GABA B current. In the simulation, the values of the external current of the TC neuron and maximal conductance from the RT neuron to the TC neuron ( I RT GABA A ) are different values for modeling the normal state ( I TC ext = 5 µA cm −2 , g GABA A = 0.65 mS cm −2 ) and CAE state with an oscillatory spiking pattern ( I TC ext = 6 µA cm −2 , g GABA A = 0.32 mS cm −2 ) 25 .
(1)  (1) is used as a CAE model (virtual patient). To implement AFTSMC by delivering the continuous DBS stimulation (control input) to suppress the oscillatory pattern of the neurons, the dynamics of the membrane voltage of the three-neuron CAE motif in (1) should be represented as a canonical MIMO nonlinear model with the following nonlinear differential equation: ., 3 is a measurable membrane potential of TC, CT, and RT neurons and u = [u 1 , . . . , u m ] T is the control input (continuous DBS pulses). d(t) denotes unknown external disturbances and uncertainties of the system, which is bounded with a positive unknown value, i.e., d(t) < υ . The unknown continuous vector functions f (x, t), f τ (x τ , t) and G(x, t) are defined as where χ 0 is a positive real parameter and I m is an m × m identity matrix.

Assumption 2
The desired trajectory x d (t) is a continuous membrane potential of TC, CT, and RT neurons in a normal condition that is measurable and its first-order dynamic exists.

Assumption 3
The time delay τ is known and measurable.
(3)  www.nature.com/scientificreports/ The tracking error of the conductance-based model is defined as e i = x d i − x i . To design AFTSMC the nonsingular terminal sliding surface is considered with the following function:

RT synapses
where σ > 0 and 1 < η < 2 are positive design numbers. If the initial states of the system are outside the sliding surface, a fast TSM-type reaching law is designed to ensure the convergence of the system states to the sliding surface in finite time as follows: where the matrices K 1 = diag(k 11 , ..., k 1m ) > 0 m×m and K 2 = diag(k 21 , ..., k 2m ) > 0 m×m are design control parameters and 0 < ρ < 1 . According to Eq. (6), the first derivative of the sliding manifold is The control input is designed in the following form:

Lemma 1 If a Lyapunov description of finite-time stability is given as follows:
where α 1 , α 2 > 0 and 0 < < 1 , then, the settling time is calculated by the following equation 61 : Lemma 2 If a 1 , a 2 , ..., a n are all positive numbers, and 0 < p ≤ 2 , then the following inequality always maintains 61 : (3), the assumptions (2) and (3) are held and the nonlinear dynamics of the system with the external disturbances are assumed to be known. The integral-based terminal switching surface is considered as (6), the control input is designed by (9), and the finite-time reachability condition to the sliding surface is ensured with (7). If the dynamic states of the model with any initial conditions are outside the sliding surface, then the finite-time convergence of the model dynamics to the switching surface s i = 0 is guaranteed. On the switching surface, the tracking error decreases to zero in a finite time.

Theorem 1 Consider the dynamical MIMO system with delay in
Proof Consider the first dynamic of the tracking error with ė(t) =ẋ d (t) −ẋ(t) , the first derivative of the sliding surface vector using (8) is rewritten as By considering the Lyapunov function V = 0.5s T s and substituting (13) in its first dynamic, it yields where K 1 = ησ diag(|e| η−1 )K 1 and K 2 = ησ diag(|e| η−1 )K 2 are modified gain matrices. Using Lemma 2, we have where k 1 = min j k 1j > 0 and k 2 = min j k 2j > 0, j = 1, ..., m are minimum eigenvalues of matrices K 1 and K 2 , respectively, and 0.5 < (1 + ρ) 2 < 1 to satisfy the condition of Lemma 2. The finite-time convergence to the sliding surface s = 0 is guaranteed based on the criterion presented in Lemma 2 and using Eq. (11), the finite time for reaching the sliding surface is as follows: www.nature.com/scientificreports/ After the reaching of the states to the sliding surface s = 0 , according to the terminal sliding surface considered in (4), the convergence time required to achieve the global stable attractor of (6) ( e i = 0 ) for any initial condition x i (t r i ) is finite and is computed as follows: This completes the proof of Theorem 1.
Adaptive fuzzy terminal sliding mode control (AFTSMC). Since in the real DBS systems, the dynam- G(x, t) are unknown, the control input (9) cannot implement. In the current study, to eliminate the limitation, FLS is applied to approximate the unknown nonlinear dynamics of the system.
and G(x, t) , respectively, the control input (9) can be modified as follows: where ε 0 > 0 is a very small parameter. In (18), the regular form Ĝ (x, ψ t g ) −1 is used to overcome the singularity of Ĝ (x, ψ t g ) and thus the well definition of control input (18). A detailed description of FLS is provided in Appendix B of the Supplementary Materials. Let us define the minimum approximation errors of the fuzzy estimator as , and Ĝ * (x, ψ * g ) are the optimal estimated functions. We have used the following inequalities to limit the values of the minimum approximation errors: where ε f i , ε f τ i and ε g are unknown small values and will be estimated with adaptive rules. To online designation , the adaptive parameters of the fuzzy estimator in (S.24) and (S.25) should be updated online with the following adaptation rules: , and ζ g ij (x) denote fuzzy basis vector fixed with the designer (see equation S.23 in Appendix B of the Supplementary Materials). To guarantee the stability of the closed-loop system in the presence of approximation errors, external disturbances, and regularized inverse Ĝ (x, ψ t g ) , the corrected control signal u r (t) is added to the control input (18) as a robustifying term where u r (t) is designed as follows: where k 1 and k 2 denote the minimum eigenvalues of matrices K 1 and K 2 , respectively. (28) and (29), confirm that the region �s� ≤ δ = min(δ 1 , δ 2 ) will be achieved in finite time and then, the tracking error decreases finite time to a boundary layer The proof is given in Appendix C of the Supplementary Materials.

Remark 1
According to the bounds of (28) and (29), the matrices K 1 and K 2 should be large enough to decrease the boundary δ . This increases the amplitude of the control input (pulse width or pulse amplitude of a DBS pulse) which limits the implementation of the controller in a real FES system.

Remark 2
The term sig(s) ρ in the reachability law embedded in control input (18) and η in sliding surface (6) are considered a bridge between linear SMC ( ρ → 1, η → 1 ) and TSMC. These parameters should be selected appropriately to guarantee the finite-time convergence of the tracking error to the boundary layer and to achieve the control input without any singularity or chattering.

Remark 3
In this paper, to satisfy the sufficient controllability condition of (3) Figure 3 shows the typical results of the membrane voltage of TC, CT, and RT neurons in response to the external current of TC neuron for normal and CAE states with the different values of maximal conductance from RT neuron to TC neuron. If the external current of the TC neuron increased from 5 to 6 μA cm −2 with the decrement of the strength of the maximal conductance from the RT neuron to the TC neuron from 0.65 to 0.32 mS cm −2 , the rhythm of the neurons changes from bursting or single spiking to tonic spiking. Indeed, the complexity of the dynamics of the neurons is dependent on the values of the external TC current and the maximal conductance from the RT neuron to the TC neuron. The smaller selection of the maximal conductance, the faster rhythm of the neurons in the CAE state. To calculate the tracking error, the desired trajectory is the membrane voltage of the TC, CT, and RT neurons in the normal state (Fig. 3a). The goal of the AFTSMC is to automatically generate the continuous controlled current (DBS pulses) to return the oscillator spiking patterns of the three neurons to the normal bursting patterns in the presence of uncertainty and disturbance current.   Figure 5 shows the results of the continuous (chattering-free) control input without any excessive control effort. When the membrane voltage of the three neurons is at rest, the controller generates the DBS current with the minimum energy and only in response to the spiking pattern of the desired voltages increases automatically the amplitude of stimulation pulses to a sharp value to rapidly track the peak of the action potential of the desired trajectory. In contrast, due to the discontinuity of the classical AFSMC in the reaching law of the control input (results of Fig. 6) across the sliding surface, the high switching control activity is appeared in the control input and degrades the performance of the tracking and causes the problem in the implementation of the control input in real-time applications.
An important issue in controlling the spiking behavior of the CAE model is evaluating the effects of external disturbances on the tracking performance of AFTSMC. The sources of external disturbances may include the un-modeled complexity of the external ionic currents, the unknown synaptic currents from the other neurons in the network applied to the CT, TC, and RT neurons, etc. The proposed controller should robustly handle the effect of these external disturbances without degradation of the tracking performance. To evaluate the ability of the proposed AFTSMC to reject the external disturbance, a non-regular constant current (pulse amplitude,  Figure 7 shows the results of the spiking pattern of the membrane voltage of TC, CT, and RT neurons for the CAE state in the presence of the external disturbance current applied to the CT neuron. The results in comparison with the results of Fig. 3(b) indicate that as the random current pulses are applied to the TC neuron, the inter-spike-interval nature of the three neurons changes so that the spiking pattern occurs following the external disturbance current is applied while sometimes as the current pulses are suddenly exerted, the spiking behavior is not observed or is delayed. Figure 8 shows the typical results of the tracking under the disturbance current of the CT neuron with a pulse amplitude of 1.25 μA cm −2 . The RMSE for TC, CT, and RT neurons are 0.3080, 0.2694, and 1.6736 mv, respectively. The controller could effectively adjust the control signal of the TC neuron to reject the effect of the disturbance current applied to the CT neuron. Results of the mean RMSE (± one standard deviation) obtained over 15 trials of the simulation for control of the CAE state in the presence of non-regular external disturbance current applied to CT neuron using AFTSMC in comparison with the AFSMC and super-twisting sliding mode control (STSMC) 70 are depicted in Fig. 9. The initial conditions of the gating variables of ionic channels and membrane voltage of TC, CT, and RT neurons, the desired pattern of the membrane voltages, the initial weights of the FLS estimators, and the onset of the pulses in the external current applied to the TC neuron were set in random values and changed subsequently from one trial to another trial of the simulation. The mean RMSEs  To evaluate the robustness of the controller to handle the effect of the time-varying uncertainty of the system parameters, the values of maximal conductance and reverse potential of the ionic currents of three neurons were changed randomly over 0-50% range about the nominal values (last column of the Tables S1, S2, and S3 of Appendix A of the Supplementary Materials) over the total time of the simulation. The variation of each parameter was obtained by passing the uniform distribution sequence with an appropriate selection of its parameters (minimum and maximum) to the second-order low-pass filter with the natural frequency of 0.05 Hz and adding the constant nominal value of the parameter. Figure 10 shows the tracking results and control pulses generated with the proposed controller over 25% uncertainty. The RMSEs over one trial of the simulation are 0.2428, 0.3613, and 1.3995 mv, respectively for TC, CT, and RT neurons. The interesting observation is the automatic generation of positive and negative amplitudes of DBS pulses for tracking control of the RT neuron. In contrast, the amplitude of the DBS pulses generated with the controller in the results of Figs. 5 and 8 was only in the negative range. Figure 11 shows the results of the mean RMSE (± one standard deviation) obtained over 15 trials of the simulation without uncertainty in comparison with random variation of 15, 30, and 50% uncertainty of the maximal conductance and reverse potential of the ionic currents of three neurons in each trial of the simulation

Discussion and conclusions
This paper developed a three-input (continuous DBS pulses), three-output (membrane voltage of neuron) modelbased controller using robust AFTSMC to suppress the spiking pattern in a CAE dynamical model. The proposed model consists of cortical, thalamic relay, and reticular nuclei neurons. Although in the previous work 48 , a finitetime fractional order SMC in combination with RBFNN was used to control the epilepsy seizures in a dynamical SWD model, prior knowledge of the model dynamics was required to design the DBS pulses. In contrast, in the current paper, the unknown nonlinear dynamics of the current and delayed states of the model used in the control input were estimated with the adaptive fuzzy controller. A critical issue in controlling the nonlinear dynamical systems in the simulations of brain networks is evaluating the effects of external disturbances and uncertainties on the tracking performance of the closed-loop system. The sources of external disturbances and uncertainties may include the un-modeled complexity of the external ionic currents, the subject-to-subject variability of the parameters, the unknown dynamics (interaction between the nuclei, the unknown different threshold firing levels of neurons in the population), etc. In the current study, the upper bounds of the external disturbance and minimum estimation errors were updated online with adaptive laws without any offline tuning phase. Although increasing the effect of the external disturbance or uncertainty degraded the tracking performance, the results of Figs. 9 and 11 demonstrated the efficiency of the proposed controller in handling the effects of the external disturbance and uncertainty of the model parameters in comparison with the AFSMC and STSMC. Moreover, according to the results of Figs. 8 and 10, the control signal (pulse amplitude) required for the proposed AFTSMC increased sharply in a normal range only when the spiking pattern of the three neurons was generated. Nevertheless, the DBS pulses showed very stable values without any extra activity in comparison with the results of Fig. 6(b) using AFSMC when the pattern of the membrane voltage was in the rest state. But, designing the automatic control pulses using the canonical MIMO nonlinear model is a significant challenge due to the singularity of the www.nature.com/scientificreports/ estimated matrix gain of the control input during its inverse calculation. To address this limitation in the current study, the regularized form of the estimated matrix with a proper selection of the parameter ε 0 was used. Another challenge facing the MIMO structure of the controlled system is the adjustment of several control parameters to achieve high tracking accuracy and robustness in the presence of external disturbance and uncertainty. To reduce the complexity of the centralized MIMO closed-loop controller with high calculation and easily select the control parameters, a decentralized form of the controller may be an alternative practical approach for the simple implementation of the controller in real DBS systems. The superiority of the decentralized structure is that the MIMO dynamical model is converted to the separated single-input single-output (SISO) subsystems with each controller acts solely on its subsystem, and the interactions between the subsystems are approximated with the adaptive fuzzy controller embedded in each subsystem. The performance of the decentralized closed-loop DBS using AFTSMC should be further investigated and constitutes future research.
In the current study, the structure of the MIMO-controlled system was designed based on the highly nonlinear (conductance-based Hodgkin-Huxley type neuron model) CAE dynamical system, and the adaptive fuzzy controller was recruited to estimate the unknown dynamics of the plant. One important issue that should be considered in a real closed-loop DBS system is the speed of the control algorithm to meet the real-time application. Regarding the computational complexity of the proposed control algorithm, with the sampling period of 0.001 ms, when the algorithm was executed on the CPU the mean run time of the closed-loop system for the total time of the simulation (2000 ms) was about 25 min. It is worth noting that a CPU sequentially processes the data. But, in real-time applications of closed-loop DBS, the control algorithm should be implemented on an FPGA (Field-Programmable Gate Array) which is an integrated circuit with so many parallel hardware resources. Implementation of the control algorithm on the FPGA benefits from parallel data processing and this leads to a great reduction in the run time of the proposed procedure, which may be addressed in future studies.
In a real DBS system, the monophasic (anodic or cathodic) and biphasic (with and without inter-pulse delay between symmetric and asymmetric anodic and cathodic phases) waveforms 76 may be generally used to deliver the stimuli pulses. But, to satisfy neural tissue safety, the pulse trains of charge-balanced biphasic currents are delivered into the neural tissues 77,78 and the stimulation parameters (pulse amplitude, pulse width, and frequency) can be adjusted individually or simultaneously with the control algorithms. Short pulse width reduces the damage risk of the neural tissue, and the increment of the pulse amplitude depends on the spatial relationship between the tip of the electrode and the target neural fibers. The high frequency of the stimulation can directly affect the average power consumption of the DBS hardware and decreases its battery life. Another problem is the spread of the electric field due to the location and geometry of the tip of the electrode during selective stimulation of the special neural tissues. The admissible amplitude or width of the stimuli pulses should be appropriately adjusted to focally stimulate the target brain nucleus with the lowest effect on the neighboring non-target tissues. In the current study, the control gain of the model in (3) was considered as a 3 × 3 nonlinear matrix and 9 fuzzy estimators were used to estimate the effect of the three control inputs on the dynamics of each neuron. Thus, we believe that the proposed control scheme can automatically adjust the pulse amplitude or pulse width of the DBS pulses with the electrodes implanted in the region of appropriate neural fibers to suppress the seizure in a real-time closed-loop DBS system.

Limitations
In the closed-loop stimulation approaches, different physiological signals are used in controlling the neurostimulator output. These data sources include LFP 79 , ECoG 80 , electromyogram (EMG) 80 and kinematic (acceleration) 81 . Selection of the appropriate biological signals depends on resolution, invasiveness, and relevance to the patient's symptoms. In previous closed-loop DBS methods, the recorded biosignal time series were not directly applied as the input to the controller. In fact, the researchers processed the recorded signals and extracted different features to control the stimulator output. Synchrony of the neural oscillations 48 and sub-band power of the brain signals 82 are examples of the considered features in DBS which are relevant to the disease symptoms. But, in the proposed closed-loop DBS technique, the recorded membrane voltage directly drives the controller without any additional processing stage. In other words, in the proposed method, the control signal varies in a continuous manner such that the output of the plant tracks the desired membrane voltage. Since the stable intracellular recording of the action potentials requires periodic calibration of the experimental setup, the long-term clinical application of this approach is limited. To eliminate the limitation of the tracking-based controller and make it suitable for clinical applications, the time series tracking error should be substituted by an appropriate feature tracking error. A candidate feature can be the synchrony of the LFP signals of the brain nuclei. Another limitation of the proposed control scheme in designing the control input is the delay of the model states. In this paper, to estimate the nonlinear dynamics of the delayed states in the input of the fuzzy estimator it was assumed that the value of the delay is known. To hardware implementation of the proposed control algorithm in a real-time closed-loop DBS system, since the delay in the nonlinear dynamics of the system is unknown it should be approximated online using the time-delay estimation (TDE) methods 83,84 such that the stability of the closed-loop system including TDE algorithm is guaranteed.

Data availability
No datasets were used in the current study and all data needed for the simulation are included in this published article and its supplementary materials file. The MATLAB codes used for the generation of the results are available from the corresponding author on reasonable request.