Estimating the Parameters of the Epileptor Model for Epileptic Seizure Suppression

Epilepsy is one of the most common brain disorders worldwide, affecting millions of people every year. Given the partially successful existing treatments for epileptiform activity suppression, dynamic mathematical models have been proposed with the purpose of better understanding the factors that might trigger an epileptic seizure and how to mitigate it, among which Epileptor stands out, due to its relative simplicity and consistency with experimental observations. Recent studies using this model have provided evidence that establishing a feedback-based control approach is possible. However, for this strategy to work properly, Epileptor’s parameters, which describe the dynamic characteristics of a seizure, must be known beforehand. Therefore, this work proposes a methodology for estimating such parameters based on a successive optimization technique. The results show that it is feasible to approximate their values as they converge to reference values based on different initial conditions, which are modeled by an uncertainty factor or noise addition. Also, interictal (healthy) and ictal (ongoing seizure) conditions, as well as time resolution, must be taken into account for an appropriate estimation. At last, integrating such a parameter estimation approach with observers and controllers for purposes of seizure suppression is carried out, which might provide an interesting alternative for seizure suppression in practice in the future.


Introduction
Epilepsy is one of the most severe, intriguing brain disorders worldwide, whose causes, despite possibly linked to genetic abnormalities, centrous nervous system infection, traumas Iasemidis (2003) or even external stimuli Zhu et al. (2010), to name a few, are not fully understood yet Vezzani et al. (2011). Some of its consequences are: physical implications (increase/decrease of muscle tone, or contraction) Browne and Holmes (2008), physiological implications (headaches, mental confusion or unconsciousness), psychological or psychiatric disorders (as higher incidence of depression and psychosis), and social issues (social dislocation, embarrassment, and isolation) Fisher and Schachter (2000). Therefore, it is a sensitive subject and a major reason for decreasing the quality of life of patients.
Although there are many possible treatments for epilepsy nowadays, such as pharmacological, surgical, based on electrical stimulation Iasemidis (2003), neurofeedback Walker and Kozlowski (2005) and ketogenic diets D' Andrea et al. (2019), they, in general, achieve only partially successful results. Understanding the possible epileptogenic mechanisms may, thus, help to clarify the triggering factors of seizures and target the best available treatments. In this sense, mathematical models represent interesting options, especially because one of their main advantages consists in providing a powerful tool for studying physical phenomena (in this case, epilepsy) without necessarily conducting experimental research. This might explain why they have been increasingly well accepted by the scientific community Wendling et al. (2015).
Nowadays, there exist dynamic models that describe the neural behavior on different scales 1 , among which Epileptor stands out when it comes to representing such activity, and, specifically, the main features found in time signatures of seizures, given in terms of local field potentials (LFPs) Jirsa et al. (2014). Besides providing a greater description of seizure-like events and their possible triggering factors Jirsa et al. (2014), it is a relatively simple, computationally and mathematically speaking, easily addressable, and generalpurpose model Brogin et al. (2020).
Epileptor was developed based on an experimental model using signals obtained from the hippocami of mice, but under the premise that similar dynamic structures are shared across species, which was evidenced by carrying out a detailed statistical analysis of real electrical activity from mice, fish and humans Jirsa et al. (2014). It has led to many other relevant works, including a bifucartion analysis El Houssaini et al. (2015), seizure propagation Proix et al. (2017); Sip et al. (2021), synchronization Zhang and Xiao (2018), modulation in phase and frequency Pinheiro et al. (2020), and suppression Nagaraj et al. (2017); Brogin et al. (2020). The model can be conveniently written using the state-space representation, thus leading to six states and eight parameters, as discussed by Brogin et al. (2020).
Although Epileptor is designed based on biophysical factors Jirsa et al. (2014), it is composed of some abstract variables and parameters, which may impose some restrictions for practical usage, such as seizure control. However, Brogin and coauthors Brogin et al. (2020) introduced a new approach for designing state feedback-based controllers to suppress epileptic seizures based on this model. The approach considered that all states are known over time, but some of them cannot be measured since their physical meaning is not precisely defined according to current knowledge on this model. Recently, this challenge has been overcome considering the concept of observability Ichalal et al. (2011), which enabled the development of a state observer for estimating all unmeasured states from only the measurable LFPs Brogin et al. (2021).
The approaches developed for designing the controller and the state observer are powerful tools that may lead to promising investigations on the use of active control techniques for epileptic seizure suppression. However, for practical applications, it is necessary to obtain a patient-specific or condition-specific Epileptor model, which represents each individual for whom an active controller is designed. To do so, a specific set of parameters must be found such that Epileptor's dynamics can be defined and represent features of realistic seizures 2 . A recent work in the literature Hashemi et al. (2020) partially addresses this problem by using Baysian inference to estimate the epileptogenicity parameter ( x 0 , which represents the degree of excitability of the system). The authors consider a reduced-order Epileptor, leading to a set of 2 equations ( x 1 , the LFP for the first subsytem, and z, the control variable) and 3 parameters ( I rest1 , 0 and x 0 ). In this work, a more comprehensive approach is proposed for estimating all of the 8 parameters in the full-order Epileptor, as shall be discussed, which can nonetheless be applied to the reduced-order as well to estimate the remaining ones ( 0 and I rest1 ). Besides, this strategy is consistent with all the previous works involving state feedback control/state obversers, which are designed considering all of these parameters.
It is worth mentioning that there is a wide range of seizures with different dynamics. In fact, their categorization has been of particular interest in the fields of neuroscience, and there are works in the literature which propose a taxonomy of seizures based on their types of bifurcations Jirsa et al. (2014) or the so-called dynamotypes Saggio (2020). The latter one even discusses how such dynamotypes evolve during patients' lifetime, display complex variations, as well as how often they may occur (individually or together). Furthermore, as mentioned above, variations of Epileptor exist today, such as the reduced-order model Proix et al. (2017); Hashemi et al. (2020), or a model considering longer interictal periods Proix et al. (2014). The present work focuses on obtaining the set of parameters that was originally defined for Epileptor Jirsa et al. (2014). Nevertheless, as long as any of these systems of equations can be converted into state-space notation, the proposed approach can be applied.
The approach chosen for this article is based on the Nelder-Mead (NM) optimization algorithm, introduced by Nelder and Mead (1965). The application of heuristic algorithms, such as the NM or the Particle Swarm Optimization (PSO) Poli et al. (2007), for example, have been widely used in the fields of engineering for purposes of parameter estimation Neto et al. (2012);Hardt et al. (2021); Altan (2020); Altan and Parlak (2020). Besides, their use in biomedical engineering, especially in the context of epilepsy, is relatively common, involving, for instance, epileptic seizure prediction Dollfuss et al. (2013); Kim et al. (2014);Hamad et al. (2018); Subasi et al. (2019) or classification Abuhasel et al. (2015); Li et al. (2017), based on the use of one of the algorithms alone or hybrid approaches. In fact, the NM algorithm has been even employed in an optimization problem concerning seizure propagation speed in the Epileptor model Proix et al. (2018). However, applications for determining Epileptor's parameters in the time domain, mainly for seizure mitigation purposes, still lack in the literature.
In this work, the NM algorithm is not directly and simply applied; an efficient strategy to obtaining the state-space representation of the Epileptor model for an individual undergoing an epileptic seizure is introduced. The main idea is to measure the LFP during the epileptic seizure and use it as input to the algorithm. The strategy involves a pseudo-random matrix and a normalization criterion to apply the algorithm successively until the best set of parameters is attained to describe the input LFP through the Epileptor model. Furthermore, since there is an inter-dependency between controller, observer, and parameter estimation, the latter one plays a crucial role when conceiving possible therapeutic strategies based on feedback control. In practice, once the parameters of the Epileptor model are properly identified, the gains of the controller and state observer can be accurately computed, as introduced by Brogin et al. (2020) and Brogin et al. (2021), respectively. Therefore, this work presents how an integration process between these three stages can be carried out as well.
At last, although the Nelder-Mead algorithm involves a heuristic approach, which may lead to local minima or maxima 3 , this problem is also addressed in this work by using a successive optimization approach in the domain time, combined with a statistical analysis and median values to assure that the estimated parameters converge towards the reference ones. Therefore, it is regarded herein as a simpler, more straightforward, and only output-dependent approach, among other advantages, as discussed in the next sections. For these specific reasons, it represents an attractive option with a strong practical appeal.
The results show that it is possible to estimate Epileptor's parameters successfully based on an output LFP. Besides, there is a link between the estimation of biophysical factors and the two main conditions described by Epileptor: interictal (healthy) and ictal (ongoing seizure). In general, slow-evolving variables are better estimated before the seizure onset (in the interictal region), whereas fast-evolving variables are more accurately estimated after seizure onset (during the ictal region). Moreover, the number of windows used during the optimizations (i.e., time resolution), the uncertainty factor (which is related to how close from the reference values the estimated parameters are), and the presence of noise, also influence the parameter estimation process.

Methods
The strategy proposed in this article consists of a process with four main stages: i) obtaining the input signals (i.e., the LFP), ii) pre-optimization, defining initial guesses for the parameters, iii) optimization for parameter estimation (modeled as a minimization problem), and iv) validation. The third stage characterizes the optimization itself, during which the unknown parameters of the model are estimated via an objective function based on the LFP. As for the last stage, a validation process is carried out considering several successive optimizations along the LFPs over time, whose median values (as shown herein) tend to converge towards the target ones. Figure 1 depicts an illustration of the aforementioned identification process.
Note that for practical applications of epileptic seizure suppression, the identified state-space model is used for designing the controller and state observer, as introduced by Brogin et al. (2021Brogin et al. ( , 2020. Further information is included herein to clarify the reader.

Obtaining the Input Signals Using Epileptor
Epileptor is a mathematical dynamic model that represents the main features of epileptic seizures, such as fast oscillations and spike or spike and wave events, which are two of the most common features in time signatures of epileptiform activities Jirsa et al. (2014). Under the premise that dynamic structures in such activities are shared across species, it was developed by Jirsa et al. (2014) based on Parameter Estimation Fig. 1 Flowchart for the proposed parameter estimation process, comprising the four main stages: obtaining the input signals as references, pre-optimization used to define the initial guesses of the parameters, the optimization itself, and validation an experimental model considering the electrical brain activity of mice. Tools from nonlinear dynamics, as limit cycles, bursters and phase portrait analysis were largely used to propose such a model, whose study led to a taxonomy of seizure-like events depending on the types of bifurcation found in real signals. In this sense, Epileptor presents, according to the same authors, two specific ones: saddle-node bifurcation (for seizure onset) and homoclinic bifurcation (for seizure offset), modeled in terms of local field potentials, usually given in V. To represent all of these features, a set of six first-order differential equations was proposed, given in state-space notation by: where nl , and are respectively given by where: The matrix is composed of a linear combination between x 1 (t) and x 2 (t) (i.e., −x 1 (t) + x 2 (t) ), which is the observable activity (t) that can be measured over time: A biophysical factor is associated to each of these variables and parameters, which are briefly described as follows Jirsa et al. (2014);Brogin et al. (2020Brogin et al. ( , 2021: First Subsystem ( x 1 and y 1 ) This set of variables is responsible for generating the so-called fast oscillations, that is, the activity with the fastest time scale ( 1 ) among all of the (1) variables. The measurable activity over time is x 1 , whereas the unmeasurable one is given by y 1 . Together, both account for the first important dynamic structure in Epileptor, characterized by an inward-spiralling attractor in the phase portrait. They are associated to neuron membrane electrical activity, modeled as local field potentials, which can be given in V or mV.
Second Subsystem ( x 2 and y 2 ) This set of variables, in turn, is responsible for generating spike and wave events, which is an activity with a relatively slower time scale ( 2 ) compared to the first set. The measurable activity over time is x 2 , whereas the unmeasurable one is represented by y 2 . Similarly, they intend to model membrane potentials, given in V o mV, but their dynamic structure is evidenced by an open limit-cycle constant oscillation.
Permittivity This entity, represented by z(t), is a control parameter; that is, depending on its values, Epileptor may be driven into or out of a seizure, which accounts for its bi-stable behavior. It is possibly related to slow-evolving extracellular processes that occur during epileptiform activity, such as levels of ions, oxygen and energy metabolism. Its time scale is represented by 0 . Therefore, the following relation holds: 0 > 2 > 1 .
Electric Currents I rest1 and I rest2 represent the flow of electric currents that flow inward or outward neural cells, for the first and second subsystems, respectively.
Epileptogenicity The degree of epileptogenicity is a parameter related to excitability, that is, after a critical value Epileptor is considered epileptogenic, thus triggering seizures autonomously Proix et al. (2017); Jirsa et al. (2017); Hashemi et al. (2020). Jirsa et al. (2014), there is evidence that both fast oscillations and spike and wave events are usually coupled. This behavior is thereby captured by the following functions: g ′ 1 and g ′

Pre-Optimization
Based on the Epileptor model, the parameters considered correspond to every entity that likely has a biophysical meaning and is not a variable or a purely mathematical constant (e.g. c from Eq. (7)), or influences significantly its behavior. Therefore, from the point of view of the system's dynamics (Eq. (1)), it is necessary to define 6 initial conditions (for the 6 variables considered in state-space representation) and 8 parameters which, if changed, also change Epileptor's behavior. By doing so, the set of 8 parameters that characterize a local field potential can be expressed as : or, using an equivalent notation, for simplicity, = {p (1) , ..., p (8) } T , where p (1) ≡ 0 and so forth until p (8) ≡ I rest2 . Note that all of them must be known before using the Epileptor model. In practice, this means the set is patient-specific; that is, each individual presents their unique seizure signature, which, in turn, implies that estimating these parameters leads to the full description of their epileptiform activity, for a given time window recorded over time. Therefore, the pre-optimization stage concerns the definition of initial guesses for each of them. In this case, the set of parameters over the n optimizations can be described by where 8 × 8 is an identity matrix, is an uncertainty factor, n is the disturbance matrix, which changes for each n-th iteration. Vector ref contains reference values, which in practice is unknown, but for the proposed approach it is taken as the target or desired parameters, to simulate a real scenario. At last, n is the input matrix for the n-th iteration. Note that, if no iterations are carried out, n = 0 , and 0 ≡ + 0 ref ; therefore, the estimated vector is ≡ 0 , and the disturbance + 0 ≠ is the highest possible. On the other hand, as the iterations proceed, n → N , where N is the defined maximum number of iterations, the estimated vector is n → ref . That is, ≡ n=N is sufficiently close to the reference ref , and + N ≅ .
The reference vector ref can be arbitrarily defined. However, it is recommended to use values typically employed to describe representative epileptic seizures to improve the feasibility of the optimization process. In this sense, without loss of generality, the values used by Jirsa et al. (2014) are considered. The robust identification is achieved by conveniently computing the disturbance matrix n of each n-th iteration by the following equation where diag( r ) is a diagonal matrix constructed with the elements in the vector n = {r 1 r 2 r 3 r 4 r 5 r 6 r 7 r 8 } T , and n is a randomly-generated vector for the n-th iteration. max(| n |) is the highest value in | n | . The highest element r j is obtained in a different position into the diagonal of the matrix n . The disturbance factor is previously defined considering > 0 , such that = 1 is the maximum percentage of variation admitted to the parameter in n in relation to the reference value in ref . For example, if = 0.25 , n includes 25% of variation in p j (i.e., the j-th element into n ) in relation to its reference value defined in ref . Figure 2 shows a schematic illustration of this approach, in which D (j) is the j-th element from the diagonal of n . Note that The process of estimating the actual parameters is also subjected to disturbances related to the presence of noise in experimental signals, which is usually associated to the output variables, i.e., the physical quantities that are measured in the time domain ( (t) ). Thus, along with the uncertainty factor, the influence of noise addition is also considered in this work as: where S noise dB stands for the power of the noise present in the output, S signal dB ≡ LFP ref dB is the power associated to the measurable reference LFP, and SNR dB is a pre-defined

Fig. 2
Visual illustration of the non-null values in n . Initially, the uncertainty factor for the optimization is defined as , along with the pseudo-random matrix n ; then, this matrix is normalized by its maximum value, thus leading to the disturbance n ; at last, the initial guesses n are formed by the disturbance term and the previously defined reference values . Red dashed lines represent the reference values and ⋆ represents each value of p j for the n-th iteration. The superscript j from each n indicates its j-th element, whereas n indicates the number of the iteration process during optimization. The same holds for n signal-to-noise ratio, given in decibels (dB). By applying the appropriate logarithm-linear convertion, the noise-corrupted output can be modeled in the time domain as: is an artificially generated and normally distributed pseudorandom variable, known as additive white Gaussian noise (AWGN) Liu and Lin (2006), given in the same unities as the LFP.

Optimization
The optimization process is employed to get n → from ref .
Each n-th iteration is evaluated by computing the function , the target LFP. Similarly, w n is the LFP obtained by integrating the Epileptor model in the time domain considering the parameters in n , i.e., the n-th set of parameters. The function f (W, w n ) is defined by where W(k) and w n (k) correspond respectively to the k-th value of W and w n , k = 1, ..., K , i.e., their discrete representation for computational purposes is considered.
= |W max − W min | , which is a normalization based on the maximum and minimum values of the reference LFP, and K is the number of time instants considered. Note that the target LFP is the reference local field potential to be described using the identified Epileptor model. In practical applications, it is a recording of epileptiform activity, whose parameters are unknown. The optimization procedure is defined as a minimization problem, as follows: The strategy established in the optimization procedure is as follows: for each n-th iteration, a set of parameters n is evaluated; a computational simulation for the Epileptor model is then run for such a set; in case f (W, w n ) = 0 , then w n ≡ W , and the iteration process is stopped. Typically f (W, w n ) = 0 is not achieved, and because of this, the tolerance value and the maximum number of iterations N are predefined. Then, the following criterion is assessed: f (W, w n ) ≤ , which limits the number of iterations, as is the case in practice 4 . Figure 3 shows an illustration of this iteration process.
The Nelder-Mead algorithm was chosen because it provides a good convergence speed Neto et al. (2012), it is easily addressable, its implementation is relatively simple, and presents a low computational complexity Neto et al. (2012), which makes it an interesting option, where several successive optimizations are performed, as shown in next section. It has been used for different applications, such as the tuning of Kalman filters Powell (2002), or model parameters estimation Neto et al. (2012); Hardt et al. (2021), for example. This algorithm was first established in 1965 Nelder andMead (1965) and it is also known as Downhill Simplex Method Neto et al. (2012). As mentioned in the Introduction, it is based on a heuristic approach, which may lead to local minima or maxima. However, this problem is addressed herein by adopting successive optimizations over time, as presented in the next subsection. Another important advantage is that it can be applied directly to Epileptor's output, in the time domain, which thus requires no previous stages of signal processing or manipulation.
Consider that there are S parameters to be estimated; then, as seen in Gao and Han (2012), the Nelder-Mead algorithm involves a simplex of S + 1 points (or vertices) that attempts to find the best solution in the parameter space based on the objective function. Essentially, this is attained by five main operations: 1) sort, 2) reflection, 3) expansion, 4) contraction and 5) shrink, each of which is briefly explained below 5 . Fig. 3 Visual illustration of the iteration process for a time window of an LFP. Initially, at n = 0 , both target W and optimized w n LFPs present different parameters, whose variation is defined by the uncertainty factor ; as the iterations proceed, this difference is minimized until a maximum number of iterations N or tolerance is achieved 1. Sort: the objective funtion h(.) = f (W, w n ) is evaluated at each of the S + 1 vertices based on the initial guesses given in the pre-optimization stage, and these values are sorted in ascending order such that: for 2 ≤ s ≤ S + 1 , where 0 < < 1 is a shrink factor. Figure 4 shows an illustration of the main stages of the Nelder-Mead algorithm to better understand this calculation flow. For representation purposes, the number of parameters was limited to S = 2 , such that the Simplex method is formed by S + 1 = 3 vertices.

Validation
The validation stage comprises repeating the optimization process considering M different time frames to define W. The idea is to evaluate an entire seizure, from preictal to postictal states, by cutting the signal in M parts. According to the number of points measured in the target LFP, it is possible to consider an overlap as usually done in signal processing techniques Oppenheim et al. (1997). The M repetitions provide (1) , ..., (M) sets of parameters. The set used is the median resulting value of all optimizations, i.e., Since m = 1, ..., M , the estimated parameters are considered as random variables. As such, an intuitive choice to summarize the them would be the mean, but if they are averaged using it, they risk being biased, assuming high values which do not reflect the real parameters Soong (2004), because the algorithm only aims at the best numerical solution. In this sense, medians are more stable and conservative as they are centered samples Soong (2004). Other advantages in doing so are summarized below: • Depending on the requirements of the optimization (in this case, the number of points, for example), it may take a relatively long time for the parameters to converge towards the desired ones. Therefore, shorter time windows enable to perform several shorter, faster evaluations, which leads to a finer resolution in function of the number of windows applied; • Since preictal, ictal and postictal signals have different characteristics, the best estimate for the parameters can be found individually, for each part, or as a whole, based on median values of multiple successive operations; • By estimating the parameters according to pieces of the signals, it may be possible that they are dominated by the local dynamics, which implies that they do not fully reflect Epileptor's true general behavior, a problem commonly associated to overfitting Hawkins (2004). However, even if a set of parameters proves to be the best fit for a part of the LFP but not for another one, by using median values, all of the parameters yielded from the optimization can be analyzed as independent values and converge to the desired ones. For instance: in case f m (W, w n ) > , obtained if the number of iterations exceeds the predefined one or a global minimum that does not tend to zero is achieved for some m, it is possible that one or more of its p (j) diverge considerably from its reference value. Once many time windows are adopted, the overall median result tends to compensate for particular discrepancies like this, because it is supposed to tend towards the target values without being significantly affected by extreme values.
Based on the aforementioned reasoning, "Influence of the Time Window" is designed to specifically address this issue and find the best trade-off between number of time windows and accuracy. Figure 5 depicts an illustration of this approach, considering sequential time windows over time, Fig. 4 Main stages of the Nelder-Mead algorithm. The objective function h(.) = f (W, w n ) is evaluated at each of the simplex vertices based on the the initial guesses for the parameters to be estimated. Then, after consecutive expansion, contraction and shrink operations, the iteration process is stopped considering a tolerance value of maximum number of iterations based on a reference (i.e., the target) LFP. The frame used herein considered that it is the minimum non-repeating cell for reproducing Epileptor's main stages: interictal and ictal. After the specified number of points, the same structure is obtained over and over, which characterizes the system's bistable behavior Jirsa et al. (2014). Moreover, by considering such a frame, the transitions between states are also captured during the optimizations.

Integrating Parameter Estimation and Controllers for Seizure Suppression
The problems of reconstructing Epileptor's unmeasurable variables based only on the measurable ones for purposes of seizure suppression have been addressed in Brogin et al. (2021Brogin et al. ( , 2020. Essentially, both are related to designing observers and controllers, by computing specific gains for each of the state of the model based on Linear Matrix Inequalities (LMIs) Boyd et al. (1994); Tanaka and Wang (2004). The main idea is to obtain weight matrices and such that: where ̂ is the state vector for the copy of the system, ̂ is its associated output, is the matrix related to the control input, h is the desired healthy state (to which the controller drives the system after its activation). Note that, with an appropriate choice for , (t) →̂ (t) , which implies that [ (t) −̂ (t)] → 0 , thus approximating the copy to the original system. Furthermore, as is defined, ̂ → h , driving the system to a seizure-free state.
It must be stressed that matrix nl is nonlinear. For the design of the observer and controller, its nonlinearities must be handled properly. Fuzzy Takagi-Sugeno modeling (FTSM), which aims at obtaining multiple linear subsystems based on an original nonlinear one Tanaka and Wang (2004) has been applied for these purposes Brogin et al. (2020Brogin et al. ( , 2021. To do so, based on the concept of sector nonlinearity, FTSM assures that 6 : where each i is one of the 2 N nl linear sub-models, k ( ) Tanaka and Wang (2004), where (.) indicates a combination between a mininum and maximum MF, and N nl is the number of nonlinear functions present in the dynamic matrix nl .
In the case of Epileptor, N nl = 4 , namely: replaced the nonlinearities f 1 ( (t)) , f 2 ( (t)) , f 3 ( (t)) and f 4 ( (t)) , respectively. To analyze this problem in terms of stability, the Lyapunov theory states that Slotine (1991): if there exists a positive-definite matrix , which composes the quadractive Lyapunov function V( ) = T > 0 , whose derivative is ̇ T + T ̇ < , then the following controlled plant is stable Tanaka and Wang (2004); Brogin et al. (2020): Steps prior to the validation process for the parameter estimation approach using sequential time windows to estimate the unknown parameters of the target LFP. According to Eq. (16), each window m represents a full optimization process, as described in "Optimization", leading to a set of parameters (m) , m = 1, 2, ..., M . The final set of parameters, considered the best representative one, is given by = med{ (1) , ..., (m) , ..., (M) } 6 A detailed formulation on how to obtain such an equivalence can be found in Tanaka and Wang (2004); Brogin et al. (2020Brogin et al. ( , 2021 and Appendix A. 7 In practice, once the model and nonlinear functions are available, they are measured for a specific time interval of interest, so the model is exactly reconstructed for this period of time. Further details can be found in Brogin et al. (2020). where , for simplicity, can be assumed as an identity matrix containing only ones or specific weights to each state in the case of Epileptor Brogin et al. (2020). Alternatively, To compute the gain , the following change of variables is necessary Tanaka and Wang (2004); Brogin et al. (2020): = −1 and x = , which leads to: All of the 2 N nl = 16 LMIs above must be solved simultaneously to assure stability Tanaka and Wang (2004); Brogin et al. (2020). This problem is usually handled by optimization packages such as YALMIP Lofberg (2011)  where Δ is composed of the entries f ( (t)) and f (̂ (t)) associated to the fuzzy model and fuzzy observer, respectively given by: . In the present case, as previously mentioned, i = , i = 1, ..., 16 , which, for convenience, can be an identity or weighted identity matrix, and (t) is the control input. To cope with Δ , the Differential Mean Value Theorem (DMVT) is applied Ichalal et al. (2011): where i is a constant vector in the interval defined by vectors and ̂ . This notation can be expended and converted into the error dynamics equation by considering all entries i, states j and a vector i ∈ ℜ n composed of zeros, except at its ith position Ichalal et al. (2011): Once again, the concept of sector nonlinearity can be used, such that FTSM assures that each partial derivative above can be rewritten in terms of a linear combination between its maximum and minimum values Ichalal et al. (2011): where ã ijk , k = 1, 2 is either a maximum or minimum, and v k ij (z i ) is a membership function (MF) associated to a partial derivative 8 . For convenience, an alternative notation can be introduced Ichalal et al. (2011): where q ≤ 2 n 2 . Therefore: , which is a similar representation of the classical error dynamics for linear systems Ogata (2010). In the case of Epileptor, the partial derivatives lead to 6 nonconstant values and, thus, q = 2 √ 6 = 64 linear models, as previously addressed in Brogin et al. (2021): As an example, ã 132 indicates the maximum value of whereas ã 131 indicates it minimum, for a given time interval 9 . The stability of such an error dynamics is assured by Lyapunov theory Ichalal et al. (2011), now considering obs as a positive-definite matrix, a quadratic Lyapunov function and its derivative, as follows: Converting the above equation into an LMI, it follows that: where = obs . All of the q inequalities above must be solved simultaneously to compute . If there exists , then there also exist obs and , and the error dynamics converges towards zero. This problem is also efficiently handled by YALMIP or PICOS, as previously mentioned. Figure 6 presents the integration of all the stages mentioned in this work. Two major steps are, thus, expected: an offline/preparation one (I), during which LFPs are acquired (or previously simulated), taken as reference or target signals, and used to design the observer and controller − therefore, note that there must be available data sets, which represent well patients' specific seizure signatures, prior to this stage; and an online one (II), where the gains computed in (I) are already defined and used to reconstruct the unmeasurable variables, thus leading to the full definition of Epileptor's states, necessary to activate the controller − which, in turn, proceeds to deliver the inputs for seizure suppression.

Results and Discussion
The proposed approach is evaluated considering a target LFP numerically obtained. An Epileptor model is defined considering the following reference parameters Jirsa et al. (2014) Figure 7 presents the results obtained for the proposed parameter estimation approach using successive optimizations, based on box plots, where the center value is the median ( best ). For I rest1 , I rest2 , 1 , 2 and y 0 , small time windows lead to a smaller region inside the box plot and less scattering around their respective medians, compared to those from longer time windows, which suggests that a finer time resolution is more appropriate in terms of parameter estimation in these cases. A reason for this behavior is that all of these parameters are related to variables that evolve in a faster time scale (as Epileptor's subsystems, which define the epileptiform dynamic features). A more detailed analysis on this topic shall be covered in the next subsection.

Influence of the Time Window
For x 0 , since this parameter is linked to the intensity of a seizure (in terms of amplitude and duration), shorter time windows may benefit the overall performance of the optimizations, but longer ones may comprise the information of the whole ictal state better, which is a possible reason why an intermediate value (25000) seems to be the best result. For 0 , in turn, longer time windows seem to favor its estimation. This parameter is associated to both duration and interval between two consecutive seizures in Epileptor, and thus epochs that capture the global behavior of an LFP better seem to be more appropriate choices. At last, longer time windows also favor the estimation of . This paramter does not present a precise physical meaning, which makes it difficult to relate it to a specific behavior; however, it plays the role of a slow-pass filter in the model, and slower activities Fig. 6 Integration process considering the proposed parameter estimation technique and seizure suppression. In the offline stage (I), previously acquired LFPs are used as reference signals to carry out the optimization process and obtain the best set of parameters , as seen in Fig. 5; these values are then used to design both observer and controller. In the online stage (II), the observer is activated to reconstruct Epileptor's actual states based on the estimated ones ̂ , and the controller uses the reconstructed states to deliver the input stimuli for seizure suppression might be benefited from longer time windows. As a general result, considering that all parameters were fairly well estimated regardless of the time window, the values obtained for the parameters tend to approach those originally defined for Epileptor Jirsa et al. (2014), and a finer time partitioning is beneficial for most cases.
Table 1, in turn, presents the relative mean squared errors (RMSE) between the estimated and actual parameters, as well as their median final objective-function evaluations (after M = 800, 400, 40, 8 and 4 optimizations, respectively, depending on the number of windows). Based on the function evaluation, it endorses the advantage of a relatively smaller time window to improve the estimation of the parameters. However, the RMSE does not necessarily reflect the same results. This may happen because f best (W, w) is obtained as a median final evaluation of all M windows, whereas RMSE is only related to the difference between a simulated reference signal and a new one obtained with the final set of parameters estimated (that is, it is not a median value). Nonetheless, the general behavior is that the finer the resolution (shorter time windows), the better the results are. This is possibly due to the pseudo-random characteristic of the initial guesses: since an unknown combination of them is set at the beginning of the process, slightly higher or lower values are likely to be obtained for the parameters; also, a set of parameters that best fits one time window is not necessarily the best choice for the remaining ones; however, the guesses tend to converge towards the correct values on average, as the number of available signal samples increases. Furthermore, note that, for this particular case, to obtain an RMSE lower than 1%, the threshold that relates the RMS of an LFP f best (W, w) ≤ must be such that = 10 −6 .

Influence of Interictal and Ictal States
In Fig. 8, the same parameters as in Fig. 7 are shown over the number of optimizations (alternatively, the number of time windows), but using a 250-point frame. For the parameters chosen, different characteristics are observed over the number of optimizations (which include information about the LFP over time). For I rest1 , the estimated values remain at an almost unchanged baseline, thus providing an accurate representation of its original magnitude. Such a parameter is only present in vector in Epileptor's equations, which is decoupled from the dynamic matrix, which might be why it is relatively easier to be approximated (the same behavior occurred for I rest2 ).
As for x 0 , two possible reasons may be associated with its behavior: it interacts with other parameters that must be estimated (as 0 , for example), which makes it difficult for the algorithm to yield a precise value throughout the whole signal; besides, it is associated with variable z(t), which represents processes that evolve slowly in time. Note that, for 1 , the opposite behavior is obtained, probably because it is associated (and interacts) with fast-evolving variables and parameters, as x 1 and y 0 , respectively. Therefore, each of the parameters are better estimated according to their biophysical nature in terms of time scale. Furthermore, this endorses the fact that adopting successive optimizations, including interictal and ictal states, is an interesting strategy, even if the resulting values must be averaged.
This result can be generalized, as seen in Fig. 9, where a comparison between the box plots in the interictal and ictal states is presented. For these cases, a 250-point frame is also adopted, such that the resulting box plot represents all optimizations carried out for the parameters (as seen in Fig. 8) in one of the aforementioned conditions. For 0 and x 0 , their best estimation is performed before a seizure occurs. The first one is a time scale parameter that, essentially, influences on the duration and interval between seizures. The second one, in turn, is known as epileptogenicity and influences on the intensity of a seizure, in terms of amplitude and duration. It interacts with 0 according to the state control variable z in the Epileptor model, which is a slow-evolving variable. This Vertical dashed lines represents transitions of states helps explain why their estimation is better in the interictal region, the period between two consecutive seizures.
As for the parameters 1 and 2 , they are time scales that regulate the behavior of both subsystems that define the main features of an ongoing seizure: fast oscillations and spike and wave events, described by x 1 and x 2 , respectively, which are both fast-evolving variables (the following relation always holds: 0 > 2 > 1 Jirsa et al. (2014)). Parameter y 0 , in turn, is a constant value associated to the first subsytem ( x 1 , y 1 ) in the Epileptor model. Since all of them are linked to events that take place in the ictal state, this suggests that they are likely to be more accurately estimated in this region.
Parameters I rest1 and I rest2 can be considered as particular cases: even though each of them is associated to one of the subsystems (which represent fast-evolving activities in the time domain), they were better estimated in the interictal region in the cases analyzed, not the ictal one. However, both are assumed as a constant flow of electric currents that flow inward or outward neural cells in the model Jirsa et al. (2014), and thus they do not belong to the dynamic matrix, which might explain this finding. It should be stressed that the results for the ictal period are fairly good as well and do not significantly compromise their estimation.
Finally, parameter , linked to the dummy variable u(t), presented a slightly better result for the interictal state; however, since it is not directly related to any biophysical entity and works as a low-pass filter, it can be reasonably estimated in both regions.
As a general rule for deciding which periods of time during an LFP activity described by Epileptor should be measured for the best possible estimation of its parameters, that is, which parts of the signal over time are more likely to provide the right estimates, arrows pointing to the best box plot are indicated in both Fig. 9. Essentially, the criterion used is not only how close the median values are to the reference once, but the variation around them: low fluctuations imply in higher chances of estimating the correct values. Figure 10, in turn, depicts a comparison between the estimated parameters obtained for different uncertainty factors . The box plots were obtained using a 250-point frame and combining both interictal and ictal states. Note that, as the uncertainty (that is, the relative difference in magnitude between the estimated and actual values) increases, the parameters tend to move farther from the reference ones. This is particularly evidenced by some visual features: for low , the region inside the box plot is smaller, more concentrated around the medians, whereas the outliers are less spread around them as well; for high , in turn, the opposite behavior is seen, with larger regions inside the box plots and a higher scattering in terms of outliers. It is worth mentioning, however, that even for = 0.75 and = 1.0 , the resulting medians are fairly close to the target values for most parameters analyzed.

Influence of the Uncertainty Factor ı
In practical terms, the prior lack of information about the initial guesses influences directly on the final output; that is, the less information, the worse. Nonetheless, accurate predictions can be made based on the successive optimizations method. At last, Table 2 summarizes these results, showing that not only does the relative mean square error get higher with the uncertainty, but also that the final objectivefunction evaluation is farther from zero. In other words, the resulting value of the objective function, which expresses how close the estimated LFP is to the reference signal, diverges from zero as the uncertainty increases, which leads to a less precise estimation. Figure 11 depicts the influence of noise addition when estimating the set of 8 parameters in both conditions, interictal and ictal. The box plots are obtained considering a 250-point frame and SNR dB = 70 dB. It is worth mentioning that other values of SNR (which shall be presented in this subsection) were tested and similar results were obtained. For 1 , 2 and y 0 , the same behavior from Fig. 9 is observed: since they are  parameters related to the subsystems that define Epileptor's main epileptiform features, they are better estimated in the ictal region. Parameter , in turn, is best estimated in II as before, and the same conclusions hold for this case as well. For the remaining parameters, an interesting behavior occurs: especially for the ones related to slow-evolving variables, they are not better estimated in the interictal region. Note that 0 , x 0 and I rest1 10 present a lower variation around the median values in IC. This may be explained due to the main features of each region in the time domain. Before the seizure onset (first bifurcation), the model essentially behaves as a positive DC trend, without any other strong dynamic structures, as spikes or fast oscillations. When noise is added to this part of the signal, it causes a progressive degradation of the baseline 11 , thus compromising the estimation process. During the ictal period, nonetheless, the main epileptiform structures take place and help the algorithm identificate the correct parameters more easily. Despite this peculiarity, the estimations are fairly good in both regions.

Influence of Noise Addition N(t)
Finally, Fig. 12 shows the results for the parameter estimation process considering different moderate SNR dB . The box-plots were obtained using a 250-point frame and combining both II and IC regions. Essentially, the overall behavior is that, as the power of the noise present in the LFP decreases, the process of parameter estimation is more accurately carried out. This can be particularly seen by the lower region inside the box-plots and fewer outlines as SNR dB increases. A higher SNR implies in less signal degradation Agarwal et al. (2017); Hussein et al. (2018), that is, the intensity of the measured activity, as the LFP in this case, is stronger (or overcomes) than that of the noise that corrupts the output, which thus helps the optimization algorithm find an appropriate solution more easily. Conversely, as the noise progressively kicks in, it corrupts the dynamic structures, thus compromising the estimation process. Nonetheless, for the cases analyzed in this work, even relatively lower SNR yielded fairly good results, which is evidenced by the small distance from the reference values and the median ones.
To obtain a more comprehensive and detailed analysis, especially considering the influence of noise (which is responsible Fig. 11 Influence of noise addition on the process of parameter estimation in both regions, using a 250-point time window and SNR = 70 dB. Blue dashed lines represent the reference values for corrupting the quality of the reference signal and decreasing significantly the accuracy of the optimization process), different sets were tested having one of their parameters changed. In this case, the first set is the same as the one addressed in the previous sections; the second one, in turn, is the same set, except for x 0 ; the third one is again the same set, except for I rest1 . Table 3 presents the results for these selected cases based on a 250-point frame, = 0.25 , and three levels of noise.
For f best (W, w) , it increases as the noise increases, moving farther from the ideal value close to zero. However, the overall result is within 0.2 (or 20%), which is relatively low. A similar behavior occurs to the RMSE: in general, it becomes higher as the levels of noise increase, leading even to percentages that exceed 100%. This can be explained by the fact that it represents a pointwise comparison between the reference LFP and the reconstructed one. That is, the signals are not perfectly synchronized, especially due to

Suppressing Seizures Based On Estimated Parameters
As illustrated in Fig. 6, after carrying out the parameter estimation methodology proposed in this work, the first requirement for applying observers and controllers is met: since both rely on knowing the initially unknown parameters before being designed, once these values are available, the necessary respective gains, for reconstructing the Epileptor's model and driving it out of a seizure Brogin et al. (2021Brogin et al. ( , 2020, can be properly computed. Therefore, the following strategy for seizure suppression is adopted: the estimated parameters are directly used to design both controllers and observers; the gains obtained from the design are then used along with inputs to each of the states of the system, which characterizes the state feedback control approach Brogin et al. (2020). Note that the design is performed with estimated parameters, whereas the activation of both observer and controller is for the system with the actual values of the parameters. It is expected thus that, albeit the set best is different from the real condition, it is sufficiently close to it so that seizure mitigation can be attained effectively. Figure 13 presents the aforementioned strategy, including the reference LFP (W) and the resulting one (w) after the successive optimization process, constructed with the estimated values using a 250-point time window and = 0.25 . Note that the curves present a good fit, and the only major visual difference is related to a horizontal stretch from both sides, leading to a positive delay over time. Nonetheless, such a discrepancy does not significantly compromise the controller's performance: as soon as it is activated, it can efficiently drive Epileptor out of a seizure, back to the healthy condition, expressed by the baseline at a higher LFP level.
It is worth mentioning that, although the relative mean square error for this condition is around zero ( RMSE = 0.63 %, from Table 1) and the difference between the expected output signal and the target one is slight, it can still be visually well identified. This suggests that the threshold must be as low as possible for an accurate estimation of Epileptor's parameters to be made. Due to this fact, close attention should be paid to those whose variation affect significantly slow-evolving processes in the time domain and, thus, are directly related to driving Epileptor into and out of seizures: x 0 and 0 . Mathematically Slotine (1991) and physically Jirsa et al. (2014) speaking, z(t) is a control variable, which explains this behavior, and since it presents both of these parameters linked together (as in Eq. (2)), it might be difficult for the optimization to yield highly precise values. Nonetheless, these results highlight that both interictal and ictal states are relevant when applying parameter estimation approaches to the Epileptor model. Furthermore, since both fuzzy observer's and controller's designs are carried out using estimated parameters, to assure that Epileptor's states are indeed properly reconstructed by the first one, as well as efficiently driven out of a seizure by the latter one, the Euclidean norm of the dynamic error ( ∥ (t) ∥ ) and the one associated to the actual states in relation to the healthy condition ( ∥ (t) − h ∥ ) were calculated and presented in Fig. 14. This approach has been addressed before, and comprehensive results can be found in Brogin et al. (2020). The seizure-free condition h can be defined as the vector composed of the states that are immediately after the seizure offset (desired behavior) Brogin et al. (2020). Both norms are computed based on the results of Fig. 13. Note that, after approximately 1000s, ∥ (t) ∥→ 0 ; also, right after the control input is activated (at 1200s, for 10s), There is a small fluctuation on the error Fig. 13 Comparison between the reference LFP (W) and the one using the estimated parameters obtained from the optimization process proposed in this work (w), as well as the controller's action, over time. The control is activated at 1200 seconds, as indicated by the vertical dashed lines, and kept on for 10 seconds, indicated by the red dotted box, based on the estimated states reconstructed by the observer * * W w Control is on dynamics at the seizure onset (that is, when the behavior is more unstable); however, as seen, this does not significantly compromise the observer's performance, and the error tends to converge towards zero.

Relevance to Real Seizure Suppression
The case analyzed in the previous section evidences the effectiveness of the controller when designed with parameters estimated at an RMSE = 0.63% . In terms of real epileptiform activity, however, a stronger background activity is present along with other artifacts and noise. Although there exist techniques for dealing with artifact removal Chen et al. (2015); Islam et al. (2015), as well as the classic techniques for signal processing and filtering Oppenheim et al. (1997), it is reasonable to assume these activities are inevitably part of real signals, which requires that the FTSM approach works well under more adverse conditions too.
This issue is addressed in Fig. 15, where the action of the controller is simulated for different conditions of the Epileptor model (as in Table 3) and their respective Euclidean norms are calculated. ref is the set of reference values (simulated) and is the set of values used to design the controller. To emphasize only seizure suppression, for simplicity, the set of parameters was manually set according to Eq. (10), and a fixed level of noise was also added to each simulation according to Eq. (13). Based on these cases, it can be seen that Epileptor is driven out of a seizure towards the seizurefree point. This implies that, even if no optimizations were carried out and the only different between estimated and actual parameters was reladed to the uncertainty , it would still be possible to design effective controllers.  The proposed parameter estimation approach aims at the best possible outcome, which is supposed to provide an accurate approximaton of the true parameters. Nonetheless, considering real applications where the number of iterations, time windows and processing are limited, finding the exact actual parameters that represent a seizure seems unrealistic, which makes it necessary to assume a variation around them. From this point of view, designing the controllers based on the uncertainty factor is an interesting strategy and demonstrates its purposes in terms of seizure suppression.
Although the necessary tools for designing controllers and observers for seizure suppression are already available, another relevant aspect on this matter is how the input stimulus can be delivered to the brain region undergoing a seizure. The state-space representation used in this work and previous ones Brogin et al. (2020Brogin et al. ( , 2021 considers the full-order Epileptor with six states, but since the fairly abstract physical nature of the model prevents the real system from receing multi-input stimuli, only one of them must be delivered, whose units are expressed in terms of electrical activity. In that case, underactuated control approaches can be applied Spong (1998); Reyhanoglu et al. (1999), where the number of states used during control design and activation is not the same. This issue remains to be addressed in future work. Nevertheless, the current status of the controller for seizure suppression purposes allows applications in computational neuroscience as of today, along with the study of seizure propagation Proix et al. (2017Proix et al. ( , 2018, for example.

Final Remarks
This work has addressed two relevant challenges involving the Epileptor model: i) estimating its parameters based on a reference (target) LFP: by using a successive optimization technique, which characterizes a parameter estimation approach, it was possible to obtain their values with a good accuracy, thus leading to the reconstruction of the signal in the time domain; ii) integrating this approach with existing methodologies for designing observers and controllers that aim at estimating Epileptor's states and providing input stimuli that suppress its seizures. For the latter one, since such a design relies on knowing the parameters that describe a seizure in the model, its application would be restricted without the contributions achieved in this work.
Considering the different conditions covered herein, some interesting results were observed. In general, a finer resolution in time (assessed using variable size time windows) may entail a better estimation. As it is common in practical applications and adopted in this work, a fixed number of iterations/function evaluations was adopted. Given this scenario, shorter time windows (that is, local behavior) may be easier for the optimization procedure to analyze and reduce the difference between the reference and estimated LFPs more effectively. Although longer epochs might be appropriate to comprise more dynamic features, the successive minimizations compensate for the restricted length of the signal, such that the median values converge toward the desired ones.
Both Epileptor's interictal and ictal periods must be taken into account when applying parameter estimation techniques because, as a general behavior, the parameters linked to slow-evolving variables are better estimated before seizure onset (interictal period), such as 0 and x 0 , which are responsible for the duration/interval and intensity of seizures, respectively, and are not directly related to the ictal period described by the model. Conversely, the ictal period is a more suitable choice when estimating parameters associated with fast-evolving variables, especially those from the two subsystems that describe fast oscillations and spike and wave events, which are the main features of ongoing seizures.
Concerning the uncertainty inherent to the initial guesses of the parameters (denoted as ), when it is set to a considerably high value, it might lead to a higher RMSE; that is, the identification process becomes less precise. As seen in the previous section, this error must be as close to zero as possible for an accurate representation of the parameters to be achieved, otherwise distorted versions of the target LFP may be obtained, which might reduce the effectiveness of observers and controllers designed for seizure suppression. However, for the cases analyzed, a higher degree of uncertainty does not severely compromise the identification process, because even for high values of (as 0.75 and 1.0), fairly good estimations were obtained.
This suggests that if an experimental LFP is taken as an input to the optimization, and assuming that an Epileptor model can accurately describe its dynamics, the unknown parameters (whose variation is also unknown) may be sufficiently well estimated by the proposed approach. Although it is necessary to define boundaries within which each parameter can vary ( ), in clinical applications it is not possible to measure them beforehand. To circumvent this issue, works addressing variations on the parameters exist in the literature El Houssaini et al. (2015), Zhang and Xiao (2018) may be useful for establishing references for the initial guesses, which is a major requirement for the optimization process to work properly and convergence to be achieved.
As for the performance of the proposed methodology in the presence of noise, it was verified that, as the signalto-noise ratio SNR is lower, it corrupts the reference LFP and the estimated parameters present a higher scattering around the median values. In particular, parameters linked to slow-evolving variables, whose suitable region in terms of parameter estimation was II, presented better results in the ictal region (IC), which is explained by the lack of strong dynamic features before the seizure onset, a region that is essentially characterized by a positive baseline. Nevertheless, even for relatively low SNR, the estimations were fairly good. This is particularly evidenced by the consistency of the results considering three differents sets of parameters, all of which presented a similar behavior and final estimated parameters close to their reference ones.
Concerning seizure suppression, although several works demonstrate the possibility of mitigating epileptiform activity by electrical stimulation, the type of stimulus may vary significantly in terms of amplitude, frequency and duration Rizzone et al. (2001); Velasco et al. (2001Velasco et al. ( , 2007; Tellez-Zenteno et al. (2006). As a matter of fact, interestingly enough, even non-periodic stimuli might attain promising results Cota et al. (2009). There exist approaches in the literature that help elucidate these questions by providing how the shape of such stimuli must be for seizure suppression, if the Epileptor model is considered Nagaraj et al. (2017); Brogin et al. (2020).
In this sense, the integration between the parameter estimation technique with observers and controllers has provided evidence that, once the initially unknown parameters that describe a seizure over time are available, that is, if the model is accurately identified as introduced in this work, they can be used to design the gains for both. Since the controller is designed based on a sufficient condition to stabilize the dynamic system, when the input control (electrical stimuli, in this case) is applied, the seizure is suppressed. In practice, real data usually contain noise, background activity and other artifacts, which may make it difficult to represent the experimental dynamics through the Epileptor model. However, the FTSM approach has also shown to be effective under more adverse conditions, where the uncertainty concerning the parameters that describe the system is fixed and no optimizations are carried out.
Nonetheless, as the cases studied in the present work indicate, estimating Epileptor's parameters is possible considering initial uncertainties and the addition of noise, which are situations often found in practical applications. Therefore, if real signals provide a good adherence to Epileptor, this implies that the model represents well the experimental dynamics, and then accurate designs of observers and controllers can be carried out for purposes of seizure suppression, according to the proposed methodology. Testing this strategy with real data may lead to feasible feedback control paradigms in the field of epilepsy, thus contributing to reducing its severe undesired effects and increasing patients' quality of life.

Information Sharing Statement
The data used in this work will be available at https:// github. com/ gmsint upon acceptance of the article. The codes used to implement the Epileptor model, the fuzzy Takagi-Sugeno observer and controller, and the optimization strategy will be available upon reasonable request. For further information, please contact the authors.

Formulation for the Linear Sub-Systems of the Dynamic Matrix A
Consider the following nonlinear dynamic system: For simplicity, if (t) = 0 , each element of ̇ can be expressed as: where f ij ( (t)) represents any function at position (i, j), N is the dimension of the system, f ij ( (t)) = f k ( (t)) , k = 1, ..., N nl < N 2 , represents the N nl nonlinear functions present in the system. According to the Fuzzy Takagi-Sugeno modeling (FTSM), this system can be rewritten as Tanaka and Wang (2004): where M ij are called fuzzy sets, r is the total number of rules, z 1 (t), ..., z p (t) are known as premise variables, and i are linear sub-models that represent the original dynamics locally. To identify them, Taniguchi et al. (2001) proposed writting all of the nonlinear functions as linear combinations between their maximum and minimum values, which leads to their exact representation for a given time interval: where f max k = max[ f k ( (t))] , f min k = min[ f k ( (t))] , and: known as membership functions (MF) Tanaka and Wang (2004), which must comply with the following restrictions: g min k ( ) ≥ 0 , g max k ( ) ≤ 1 , and g min k ( ) + g max k ( ) = 1 . Based (36) (t) = nl (t) + nl (t) where r 1 = 2 (N nl −1) , and N i ( ) = ∏ N nl k=1 g (.) k ( ), which must meet the requirement: ∑ 2 N nl i=1 N i ( ) = 1 . The superscript (.) indicates a combination between the maximum and minimum values for the kth nonlinear function g (.) k ( ) . After substituting each f k according to Eq. (42) into the matrix nl , the system in Eq. (36) is rewritten by Each j is, thus, a linear sub-model. If matrix nl and the input are considered back into the formulation: or, alternatively: (41) For the sake of understanding, Fig. 16 presents an example of the exact reconstruction of a nonlinear function using the membership functions and linear submodels. In this case: f (t) = e −2t sin(2 t) , that is, N nl = 1 . The membership functions g max 1 (x(t)) and g min 1 (x(t)) weigh the maximum and minimum values of f(t) ( f max 1 and f min 1 , respectively) over time such that its exact representation is obtained as f 1 (x(t)).

Formulation for the Linear Sub-Systems of the Error Dynamics e(t)
To obtain the membership functions (MFs) used to model the error dynamics, the concept of sector nonlinearity is applied once again Tanaka and Wang (2004); Ichalal et al. (2011). Each partial derivative can be rewritten as: assuming it is differentiable on the interval [a ij , b ij ] , where a ij and b ij are individual maximum and minimum values for each i and j, respectively. Therefore, the above equation can be conveniently rewritten as: where: which must meet the following requirements of convexity Ichalal et al. (2011):