https://orcid.org/0000-0002-1266-5902
Figures
Abstract
Epilepsy affects millions of people worldwide every year and remains an open subject for research. Current development on this field has focused on obtaining computational models to better understand its triggering mechanisms, attain realistic descriptions and study seizure suppression. Controllers have been successfully applied to mitigate epileptiform activity in dynamic models written in state-space notation, whose applicability is, however, restricted to signatures that are accurately described by them. Alternatively, autoregressive modeling (AR), a typical data-driven tool related to system identification (SI), can be directly applied to signals to generate more realistic models, and since it is inherently convertible into state-space representation, it can thus be used for the artificial reconstruction and attenuation of seizures as well. Considering this, the first objective of this work is to propose an SI approach using AR models to describe real epileptiform activity. The second objective is to provide a strategy for reconstructing and mitigating such activity artificially, considering non-hybrid and hybrid controllers − designed from ictal and interictal events, respectively. The results show that AR models of relatively low order represent epileptiform activities fairly well and both controllers are effective in attenuating the undesired activity while simultaneously driving the signal to an interictal condition. These findings may lead to customized models based on each signal, brain region or patient, from which it is possible to better define shape, frequency and duration of external stimuli that are necessary to attenuate seizures.
Citation: Brogin JAF, Faber J, Reyes-Garcia SZ, Cavalheiro EA, Bueno DD (2024) Epileptic seizure suppression: A computational approach for identification and control using real data. PLoS ONE 19(2): e0298762. https://doi.org/10.1371/journal.pone.0298762
Editor: Patricia Wollstadt, Honda Research Institute Europe GmbH, GERMANY
Received: December 22, 2022; Accepted: January 31, 2024; Published: February 28, 2024
Copyright: © 2024 Brogin et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: The data used in the present work can be found on Zenodo at: https://doi.org/10.5281/zenodo.5039479. The codes developed for all analyzes are available on GitHub at: https://github.com/jafbrogin/identification-control-epileptiform-activity.
Funding: This work was supported by: Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - CAPES (code 001, granted to E.A.C., and 88887.481049/2020-00, granted to J.A.F.B.), Conselho Nacional de Desenvolvimento Científico e Tecnológico - CNPq: 442563-2016/7, granted to J.F., CNPq/MCT-Instituto Nacional de Neurociência Translacional (INNT): 573604/2008-8, granted to E.A.C., Universidad Nacional Autónoma de Honduras (CU-O-041-05-2014, granted to S.Z.R G.). The sponsors/funders did not play any role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript. There was no additional external funding received for this study.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Epilepsy affects millions of people worldwide [1, 2], leading to several implications that reduce patients’ quality of life [3–5], which defines this disorder as a very sensitive and relevant topic of research. From the neurophysiological point of view, it is associated to the progressive hypersynchrony of neural firing over time, generated by endogenous or exogenous agents [1, 6]. However, the precise triggering mechanisms that disrupt a healthy brain leading to a seizure are not sufficiently understood yet [7]. Currently, different papers have been focusing on the development of computational models able to identify the beginning of seizures and their main time signatures [8, 9]. These models help to explain experimental observations involving seizure onset and propagation [8, 10], mechanisms of seizure synchronization and control [11–13].
From the physical point of view, epilepsy is a multivariable and multiparameter dynamic system, described as differential equations or state-space notation [8, 13]. The parameters in these equations intend to represent biophysical factors that, if modulated, may change its overall behavior [8, 14, 15]. Computationally speaking, by setting them, it is possible to pre-define the response of the system and analyze several scenarios (i.e., how the response is influenced by different parameters). In practice, if one or more of these simulations match a specific signature of real brain activity (usually given in the time domain), this implies that the set of parameters used for this case is consistent with the real activity. In system identification (SI) approaches [16, 17], this strategy is known as parameter estimation [18].
SI is a widely studied topic in some areas of knowledge, such as structure and control engineering [16, 17], where each case depends on the knowledge about the system. For example: there exist white-box models, which rely on having prior knowledge and deterministic equations that represent it [19]; gray-box models, which rely on prior knowledge and some available data (in this case, real signals from brain activity) [16, 19]; and black-box models, which only take the data into account, without any previous knowledge [16, 19].
In the case of epileptiform activity, especially involving dynamic models, SI applications in the literature are not commonly found. There exist some parameter estimation ones applied to the Epileptor model [8], whose development was an important breakthrough for a better understanding of the possible origins of seizures: Ref [20] proposed a Baysean framework to estimate the epileptogenicity, a parameter responsible for the excitability of the model; Brogin, J. A. et al. [Unpublished], in turn, introduced a successive optimization approach to estimate all of its parameters considering uncertainties and noise. A limitation of the latter approach is that it considers only a reference output of the local field potential (LFP), not real signals. Furthermore, the application is restricted to seizures that can be accurately represented by the Epileptor.
Works involving linear or nonlinear autoregressive models (AR/NAR), typical black-box modeling tools, have been mostly proposed in the literature for seizure detection, prediction and synchronization [21–25], but approaches considering system identification have also been used to describe the so-called NAR fingerprints of patients [26, 27]. However, the latter ones do not provide a systematic procedure to chacterize the AR equations in terms of dynamic models. If this is done, seizures can be described in state-space notation and used for purposes of control [13], for example.
Considering this, the aim of the present work is twofold: (1) to propose an SI approach, which comprises the use of real brain activity, recorded from human brain slices; (2) to provide a strategy for reconstructing real signals and applying control inputs to attenuate the epileptiform activity based on state-space plants. For the first one, AR models [28–30] are used, making the identification process simpler, since they only take into account the output of the system, which is often the case in practice.
For the second one, state-space (SS) models are obtained based on rearrangements of the AR equations [29] and converted into continuous-time models [31, 32]. Once the SS plants are available, they make it possible to design both observers and controllers for purposes of reconstructing the signal over time and then applying inputs to attenuate the undesired activity, respectively. These tools are commonly found in control engineering. Observers are useful for estimating abstract/unmeasurable variables from mathematical models based on the measurable ones [33, 34]; controllers, in turn, play the role of external inputs whose objective is to drive the system to a desired condition or state, thus mitigating undesired oscillations, for example [33, 34]. In this sense, a seizure can be mathematically represented according to each specific signal and activity (named in this work as customized models), and computationally attenuated.
Furthermore, the second objective of this work is analyzed considering two other different cases: the SS plant used to design both observer and controller is obtained from a signal undergoing epileptiform activity (named here as periodic ictal spiking − PIS) and the control inputs are applied to the same signal; the plant is obtained from a signal that represents interictal activity (named here as II) and applied to the PIS signal, which is taken as a hybrid controller in this work. The main purpose is to assess the possibility of designing observers and controllers without necessarily having access to seizure events, only the interictal ones. This is regarded as an interesting and important study in the sense that one does not need to rely on awaiting seizures to occur to be prepared for it. Although these events are cyclic [8, 35], they are often unpredictable [36, 37], leading to interictal periods that are typically much longer than the ictal ones [38], so hybrid controllers may compensate for the possible lack of PIS activity recorded.
The electrophysiological data used in this work were recorded from slices of hippocampal subfields that were surgically resected from patients with pharmacoresistent temporal-lobe epilepsy (TLE) [39]: Dentate Gyrus (DG), Cornu Ammonis 1, 2, 3 and 4 (CA1, CA2, CA3, CA4) and Subiculum (SUB). The results show that AR models can, indeed, represent real epileptic seizures fairly well, which is particularly helpful to provide customized models for each patient or time signature. The observers can accurately recreate the real activity over time based on the identified SS plant converted from the AR model, regardless of the approach adopted (hybrid or non-hybrid). The controllers are capable of attenuating the amplitudes of the undesired spikes (representative of periodic ictal spikes) while simultaneously driving the system back to the interictal condition (II), being the non-hybrid one more effective. At last, further implications about having these patient-specific models as well as the advantages and disadvantages of the proposed approach are also discussed in detail to evaluate the contributions achieved in this work in future practical clinical applications.
This article is organized as follows: the Methods section presents all the concepts and techniques employed for the proposed SI approach, divided into three stages: design, simulation and comparison; the next sections, Results and Discussion, comprise all the results yielded during these stages, as well as the main contributions and limitations intended in this work; at last, the Final Remarks section summarizes the main findings obtained.
Materials and methods
This work comprises 3 major stages: 1) design, 2) simulation and 3) comparison. For the first one (1), the SI strategy is carried out to obtain the AR models, convert them from discrete state-space (DSS) to continous state-space (CSS) models and design both observers and controllers. Design is composed of 10 main steps: i) obtaining raw signals, which are the brain activity directly measured during clinical procedures (either II or PIS); ii) signal processing, which in this case consists of filtering and downsampling; iii) frequency analysis; iv) AR model definition and existence; v) AR model order selection and fitting; vi) model extension; vii) discrete state-space definition (DSS); viii) convertion to continuous state-space model (CSS); ix) design of the observer; x) design of the controller.
Simulation leads to 2 other steps: i) applying the observer to reconstruct the real signal, which is a copy of the PIS/II activity; and ii) applying the controller to attenuate the epileptiform activity while driving the system (i.e., brain region) to the II condition, where no seizure takes place. This stage is performed such that the estimated states are used as input variables for the controller. At last, the comparison stage comprises comparisons between the controlled signal (i.e., mitigated PIS driven back to II) and the original II one. For this purpose, the Euclidean norm, spectrograms/power spectral density (PSD), cross-correlations, principal component analysis (PCA) and appropriate statistical tests are considered. In this way, a comprehensive time-frequency analysis is carried out to assess the efficiency of the controller.
Ethics statement
The present study was approved by the Ethics Committe of Universidade Estadual Paulista “Júlio de Mesquita Filho” (CAAE 60470222.9.0000.5402), thus allowing for the use, analysis and sharing of the dataset.
Design
Electrophysiological in vitro recordings.
The signals analyzed in this work were previously recorded from the granule cell layer of the dentate gyrus (DG) and in the pyramidal cell layer of the CA1, CA2, CA3, CA4 and subiculum (SUB), by using glass electrodes placed on slices of 12 human hippocampal specimens that were surgically resected from patients with pharmacoresistent temporal lobe epilepsy. This process was approved by the Ethics Committe of Universidade Federal de São Paulo (CAAE 47551015.1.0000.5505) with written informed constent from all the patients, and then published in [39]. For further details about the procedures, please refer to the same study.
Signal processing.
The raw signals were properly filtered to remove any uncorrelated activity (i.e., shifting baseline, spurious harmonics), as carried out in [39]. By filtering them, the resulting activity consists mainly of only epileptic spiking or bursting behavior over time. Then, a downsampling of factor 10 is applied (from 10kHz to 1kHz), which helps significantly to reduce the computational cost and time, without compromising its overall behavior.
Frequency analysis.
By analyzing each signal in the frequency domain through its spectrograms, it was possible to find out which frequency components are more pronounced in the brain activity, and an appropriate time window can be selected based on them, covering a range from fmin to fmax. Therefore, the minimum time window that must be used for the identification process is T = tmin = 1/fmin, leading to Nmin samples. It must be stressed that this value is only the first reference. As shall be shown in the next steps, sufficiently long time windows are applied as a conservative approach.
AR model definition.
AR models are equations whose current value depends on previous observations of a time series. They try to predict present states based on its time history. The classical AR model can be expressed as [16]:
(1)
where n = 1, …, N indicates a discrete-time series, ϕi, i = 1, …, k, are coefficients that weigh each past observation y(n − i), i = 1, …, k, of variable
, which in this case represents the output of the system (i.e., the elecrophysiological recordings), and ϕ0 is the mean value of the time series (or the intercept).
In this work, one of the main assumptions is that the time series is stationary [29], despite local variations represented by the spikes. As shall be shown herein, with an appropriate number of AR coefficients, the epileptiform activity can be adequately reconstructed.
In practice, no AR model is capable of fully describing the dynamics of a time series, because the number of coefficients is usually limited, and for certain time intervals there might be local non-stationary behaviors. Therefore, if the error considered during the modeling is taken into account, Eq 1 becomes:
(2)
where εn is the error associated to the n discrete time. One possible way to evaluate model consistency is by calculating the autocorrelation function (ACF) of its residuals: [30]:
(3)
where, for simplicity of notation, εn = ε(n) and εn−i = ε(n − i), cov(.) is the covariance operator, and var(.) is the variance operator. If ρi is relatively steady and close to zero for any i > 0, then εn can be assumed as a Gaussian error of the type N ∼ (μ, σ2), where μ is a zero mean and σ2 is a unity variance. If this proves to be true, then the AR model has captured the significant dynamics of the time series, and the residuals are merely random fluctuations.
AR model order selection and fitting.
Some criteria that are often used in the literature to define the order of an AR model are: the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC). Essentially, both attempt to obtain a good trade-off between the order k and the sum of squared errors (SSE) between the resulting model and the original time series. Their equations can be defined, respectively, as [29, 30, 40]:
(4)
(5)
Based on these equations, when a minimum is found on the k vs AIC/BIC plane, it represents the best order for the model. Another possible way to assess order selection is by the mean squared error (MSE), which is closely related to the previous approach. In this work, MSE is also tested to balance the order of the model and fix a limited error, and N = Nmin is used as the minimum time window.
Model extension.
This step comprises the extension of the AR(k) model to more samples from the measured signal. That is, if Nmin samples were used to fit it, N > Nmin samples are now included. In practice, since it is not feasible to measure an infinite time of epileptiform activity, this strategy is useful to extend the fitted model to future time instants. Naturally, this extension involves errors; however, as demonstrated in this work, they do not significantly compromise the system identification approach.
Discrete state-space definition (DSS).
Once the AR models are properly obtained, they can be converted into state-space notation according to the difference equation xk′ in function of the states xk (which contains the dynamic matrix composed of the AR coefficients), and the measurable activity in the time domain, represented by yk [29]:
(6)
(7)
where A{d} is the discrete SS model and C is the matrix corresponding to the measurable states in time (in this case, only the first one, yn, can actually be measured, whereas the remaining ones are delayed states, thus not having a specific physical meaning). Since low-frequency components are filtered out, b = 0. Besides, this AR modeling identifies the discrete plant and the error is evaluated a posteriori. In this sense, in practice εk = 0 is adopted for design purposes.
Convertion to continuous state-space model (CSS).
The previous step presents the AR model in terms of a difference equation, which is a common practice for a one-step ahead prediction, for example [28]. Although digitized signals are finite and taken as an ensemble of samples, a more realistic approach is to consider a continuous-time model that accounts for all time instants, as usually done in engineering applications [33, 34]. Besides, because the application proposed in this work requires observers and controllers, a numerical integration procedure is carried out in SS notation, which also justifies the need for a continuous-time plant. For this purpose, the formulation becomes [31, 32]:
(8)
where
is the continuous-time model, R = (A{d}−I)(A{d}+ I)−1 and I is an identity matrix of the same size as A{d}. By adopting the Direct Truncation Method [31], the continuous-time model can be approximated, for simplicity, by:
(9)
Note that this model avoids terms of higher order, thus reducing the magnitude of the elements in the continuous-time matrix. It can be further adapted by considering a scaling factor such that , where γ1 ≠ 0. This is particularly interesting for some reasons: (1) local or quasi-non-stationary behavior of the signals leads to higher AR coefficients in the discrete model [30], which in turn is reflected on the continuous one; thus, scaling down the models is proposed in this work to mitigate this problem; (2) this problem is endorsed by the sampling frequency fs of the signals, which is considerably high even after downsampling (1kHz) due to the time scale of neuron firing being relatively fast; (3) since the proposed approach focus on AR models of low order, local fluctuations may not be captured by them, which means that the difference between the dynamics of the recording being read and the one reconstructed by the identified AR model may not be compensated for by the gains of the observer.
Design of the observer.
Suppose that, for each of the m = 1, 2, …, M time windows of size Nmin presented in Fig 1, there exists an identified plant: for m = 1, during the Nmin first points; for m = 2,
in the interval] Nmin, 2Nmin], and so forth. Then, the local models in state-space notation become:
(10)
where t = n/fs represents the continuous time variable, such that x(n/fs)≅x(t) [31]. This system can be understood as the mathematical representation of the real PIS/II signatures in time. However, controllers cannot be directly applied to it since the signals are only samples stored computationally, not dynamic systems. In other words, there is no access to the true system set of deterministic equations that describes the electrophysiological brain activity in hippocampal subfields or as a whole, only the signals that have been previously measured. To tackle this problem, the concept of observability needs to be applied. A similar approach has already been carried for the Epileptor model [41], where the system’s equations are known [8]. Essentially, based only on the measurable local field potential, the unmeasurable variables from the model were reconstructed [41].
For each time window, is defined and, as the observer/controller procedes in time, models are switched to adequately represent the window. The gains, however, are single for the whole signal: G and L.
A state observer is commonly used when some of the states in vector x(t) do not have a precise physical meaning or cannot be measured for any reason in real applications [33]. As previously discussed, although C = [100…0]T is consistent with the fact that only the first state is measurable and the other ones are merely artificial states, for purposes of recreating the signal in time and applying a control input afterwards, it is necessary to implement an observer such that the numerical integrator recreates its behavior artificially. For this purpose, a “copy” of the identified plant is required, represented by the following equation [33, 34]:
(11)
During the numerical integration in the time domain, the matrix L weighs the relative difference between the output of the stored PIS signal y(t) and the estimated output (which is in fact the recreated signal) such that, as
,
, and the copy becomes an exact approximation of the actual output. In practice, y(t) is not a stored signal, but the real activity measured, for instance, by electrodes positioned on specific brain regions. To assess how accurate this approximation is, the error can be calculated by
[33, 34]. By taking the derivative of this equation, the error dynamics becomes [34]:
(12)
It is thus expected that as t → mNmin/fs, according to each m = 1, 2, .., M. Although the time windows switch over time, it is an interesting strategy to obtain a single L that is capable of observing the whole signal, i.e., one that is valid for all time windows. To do so, it is necessary that this gain be stable across all the identified plants
. This issue is addressed in this work by taking advantage of the concept of Fuzzy Takagi-Sugeno Modeling (FTSM), which states that a nonlinear plant can be exactly reconstructed using linear submodels represented by the so-called Linear Matrix Inequalities (LMIs) [33]. LMIs can be formulated considering the stability criterion by Luapunov theory [42] applied to the error dynamics:
(13)
By substituting Eq 12 into the derivative, it follows that:
(14)
Since the norm of the error remains positive still convergence towards zero, the equation above can be rewritten as
(15)
Therefore, in the case of the observers, the LMIs are stated as the following optimization problem [33]:
(16)
where
and Pobs is a positive-definite matrix [42]. In this optimization problem, all of these inequalities are solved simultaneously to find L before the numerical integration in the domain. If there exists a solution, then there exists M, thus also assuring the existence of L. Once L is found, it is applied in the state-space equation from. Eq 11.
These inequalities are formulated based on the Lyapunov theory, so solving them implies that there exist sufficient conditions that assure stability [42]. An important difference is that, while FTSM assumes that nonlinear functions can be rewritten in terms of linear combinations between their maximum and minimum values [33], the present approach only considers the several linear models obtained from the CSS AR indentified models. It is expected that if L exists, it can be applied to the whole signal observing it every m time window as the models switch over time.
Design of the controller.
Considering again the m = 1, 2,.., M time windows observed and recreated by the responses of the observer, a control input can be applied to each of them according to the following equation (a similar approach has already been carried for the Epileptor model too [13]):
(17)
where u(t) is the input entry for control (which in this case represents an active electrical stimulus, in terms of volts [13]) and B is the input matrix, whose elements can be set according to the design requirements or applications. The input u(t) has been proposed to be of the following form [13]:
(18)
Note that it is actually composed of two superposed stimuli: the negative term attempts to attenuate the undesired dynamics while the positive term induces the desired one, based on a predefined non-seizure behavior xh and gains coming from the matrix G [13]. Thus, the system can be rewritten as:
(19)
If xh = 0, then the controller supposes the stability point is a zero baseline; otherwise, if xh ≠ 0, it attempts to drive the current state to the desired condition while simultaneously mitigating the undesired epileptiform activity. In this work, xh is taken as a reference signal in the II condition, with the first state being one of its samples at each time instant and the remaining ones equal to zero. To be consistent with the output matrix C, it is considered as B = diag{B11, 0, 0, …, 0}, which means that the only target state (the one that actually receives an input) is the first one, whose magnitude can be conveniently adjusted [13].
Similarly, the design of the controller can also be carried out considering a single G for all of the m = 1, 2, …, M time windows over time, thus assuring stability based on the Lyapunov criterion, but now considering the states, not the error dynamics: V(e(t)) = x(t)TP x(t)>0. By applying the same procedure described for the observer, the LMIs for the controller are [13, 33]:
(20)
where Gx = G P−1. Considering that the system of interest in this case where physical inputs are intended to be applied is a brain region, close attention must be paid to the intensity of the gain G, so that it is properly delivered to suppress only the epileptiform activity with an acceptable input intensity [13]. This can be carried out considering a concept from FTSM, namely decay rate, defined as follows [13, 33]:
(21)
The decay rate α influences the positiveness of P, thus forcing the optimization to look for stronger/weaker gains [13, 33]. In this work, a good trade-off between this parameter and B is analyzed. For purposes of design, matrix B can be any convenient matrix defined by the designer. The simplest one is the identity B = I, for example. However, to be consistent with the fact that only one state is measurable in this work, the application of the controller during numerical integration requires B = diag{B11, 0, 0, …, 0}. This approach is known as underactuation in the control engineering [43]. This set of inequalities, as well as the previous ones for the observer, can be simultaneously solved using specific optimization tools, such as YALMIP [44], for MATLAB, or PICOS [45], for Python. Without loss of generality, the computational package adopted in this work is PICOS.
Hybrid and non-hybrid observers and controllers.
This work investigates the possibility of not only designing observers and controllers based on the epileptiform/PIS activity to suppress PIS activity, but also designing them using II activity and then applying them to attenuate PIS activities. This is regarded as an interesting study in the sense that one does not need to rely on awaiting seizures to occur to take action. The two types of designs carried out are: non-hybrid (design performed using PIS to suppress PIS) and hybrid (design performed using II to suppress PIS), respectively. Essentially, this implies that the internal dynamic structure of changes according to the signal from which the AR models were obtained.
Remarks on the controller.
The input can be adapted as [46]:
(22)
where γ2, γ3 ≠ 0 are constants that are applied to tune its intensity. Therefore, during the time integration the controlling effect is given by:
(23)
In sum, the continuous-time models and control inputs are weighted over time normalized by the combined constants γi, i = 1, 2, 3. The simulations showed that this approach provides more flexibility for a finer tuning of the controller, especially considering the limited order of the AR model.
Simulation
Applying the observer.
Once the gain of the observer L has been computed and stored, it can be applied during the numerical integration process. When the signal is used for observation it represents the output y(t) of the output of the system being measured (i.e., epileptiform/PIS activities that are stored or recorded online from electrodes in practice). Therefore, the copied signal (PIS) is reconstructed over time
as
.
Applying the controller.
In parallel of the G processing and storage, an II activity is also stored to be used as a reference xh. Then, as the PIS signal is observed over time, a control input Fc is applied to it in order to attenuate the undesired epileptiform spikes. At the same time, the controller introduces a second type of input based on xh. At the end of the simulation, it is expected that the reconstructed signal with the control input tends to the reference II signal:
.
Comparison
To assess how close the reconstructed controlled signal is to the desired reference (), some criteria are adopted in this work. The first one performs a pointwise comparison between both signals according to the Euclidean norm. Then, spectrograms and PSDs are applied, which are traditional tools in engineering applications to analyze the frequency content of signals. Considering this, it is expected that the Euclidean norm
and the frequency contents from both signals are similar to each after the controller is activated. Furthermore, the normalized cross-correlation is applied in the time domain; this procedure aims to analyze the maximum correlation value between two signals at different lag positions. At last, PCA is carried out using the previously obtained PSDs to provide a more robust evidence that controlled signals are closer to II ones whereas uncontrolled signals are closer to PIS ones. The characterization of similarity is performed using the Kruskal-Wallis test at a 95% confidence level, and post-hoc Dunn’s test: if statistical significance can be inferred, the signals are considered to be different; if no statistical significance can be inferred, they are taken as similar.
Fig 2 presents a flowchart with all the 3 major stages proposed in this work. Note the interdependency between observer and controller: as the first one recreates the real signal computationally, a control input is applied to the recreated artificial signal to verify whether its epileptiform activity can be mitigated. This integration approach has been proposed elsewhere, but for the Epileptor model, and consistent results were obtained by Brogin, J. A. et al. [Unpublished]. It must be emphasized that the present approach is only carried out because there is no access to the true system. In real applications, a control would not be applied to an artificial signal, but to a real brain region. It is thus expected that the simulations reflect real conditions sufficiently well to reproduce the results in clinical environments.
Design: PIS or II signals are obtained beforehand (thus defining the non-hybrid or hybrid approach) so that the AR modeling can be applied; several AR models are obtained from evenly-spaced time windows, whose size is defined during steps iii and iv; the models are rearranged into DSS and then converted to CSS; the same number of plants for designing the observers and controllers is used to generate single individual gains: G and L, respectively (it is important thus to emphasize that: several models from one single lead to single gains). Simulation: II signals (xh) are taken as the reference behavior to which drive the system after the control is activated (i.e., no seizure condition); the real PIS activity is reconstructed based on the output matrix C (identifying the measurable state) and the relative difference between measured (y) and estimated () states, weighed by the gain L; with the gain of the observer G, control inputs are applied to mitigate the PIS activity. Comparison: the controlled activity (
) is compared to the reference behavior (xh) using specific techniques: Euclidean distance, spectrograms and PSDs, normalized cross-correlation and PCA. Matrices A and B correspond to the identified CSS AR (which varies over time) plant and control input one (identifying the controlled states), properly defined in the next subsections.
Results
This section comprises the results from all stages involving the system identification approach proposed in this work, organized in subsections. For simplicity, some of them have been merged. Besides, for purposes of better explaining each of the steps, the analysis is carried out using only 2 signals, one from each condition: PIS and II. Nonetheless, other results considering more signals can be found in Fig C in S1 Appendix.
Design
Raw signals and filtered, downsampled signals.
Fig 3 presents examples of different signals obtained from the same hippocampal subfield (in this case, the subiculum − SUB) in two conditions: interictal (Fig 3A) and ictal states (Fig 3B). For the raw signals, there is a moving baseline over time, which is likely due to an external source of voltage that is not inherent to the spiking activity. This issue is addressed by implementing a high-pass filter, which in this case was at 0.5Hz.
(A) and (B) contain examples of raw and filtered downsampled signals from SUB in the II and PIS conditions, respectively. For the filtering, a high-pass first-order filter is applied at 0.5Hz. For the downsampling, a factor of 10 is applied in this case, which means that the initially 10000Hz-sampled series becomes a 1000Hz-sampled series. (C) presents the spectrograms for the same signals, both obtained after filtering and downsampling. A higher range of frequency is reglected because it does not contain further strong frequency components.
This process is not only important to unveil the spikes and bursts, but also necessary so that the mean properties of the signal are better kept over time, which is a condition for the AR models to be properly estimated [29]. Although the downsampling effect cannot be explicitely seen, a factor of 10 was applied to all signals, such that the sampling frequency becomes 1000Hz instead of the previously defined one, equal to 10000Hz.
The main difference between ictal (PIS) and interictal (II) spikings concerns the interspike intervals: while the PIS activity presents more regular evenly-spaced spikes, the spikes of the II activity are more irregular, switching from local bursts to isolated spikes over the whole signal. A previous statistical analysis and classification of the signals used in this work has been carried out [39]. Essentially, according to the authors, the patterns exhibited in PIS signals indicate a more deterministic activity and underlying mechanism, whereas II signals suggest a more complex/chaotic one. Further considerations on this matter are out of the scope of the present work, but some of these findings are important and shall be recovered in the next sections.
Frequency analysis.
Fig 3 also shows the spectrogram of both signals analyzed above, after filtering and downsampling (Fig 3C). There is an increase in terms of power for the ictal activity over the whole period of time in question, which is evidenced by the darker pixels. However, in terms of frequency components, their ranges remain relatively the same for both conditions, comprising approximately from ∼0.5Hz up to ∼40Hz. A higher range of frequency is not shown in the same figure because it does not contain strong components.
As a conservative approach, it is assumed in this work that: fmin = 0.5Hz. Therefore, tmin = 2 seconds, at least, to assure that all spectral activity is properly represented during the AR model definition, which leads to Nmin = 2000 points considering fs = 1000Hz, the sampling frequency for the downsampled series.
AR model definition, order selection and fitting.
Since the definition of an AR model is closely linked to its fitting, which, in turn, is based on order selection, these steps are carried out together. To define the AR model, Nmin comes in handy, because it has already been calculated as an adequate time resolution. As a conservative approach, this number of points is increased by a factor of 2.5 in this work, leading to Nmin = 5000. By doing so, it is also possible to comprise different features over a longer period of time, thus providing a source of heterogeneity to the model.
Fig 4 presents the results for the k vs AIC, BIC or MSE criteria, considering the two previously studied signals, in the interictal (Fig 4A) and ictal conditions (Fig 4B). Note that, as the number of k increases, all of the measures for the criteria decrease. Additionally, the difference between previous and past values tends to become less and less abrupt, indicating that the AR(k) coefficients of higher order have a weaker effect in improving significantly the quality of the model. At last, both criteria provide similar results, which is evidenced by the superposition of both curves.
(A) contains the results for the II case whereas (B) constains the results for the PIS case: represents the reference II/PIS signals;
represents the respective fitted AR(k = 6) models. Order selection:
indicates AIC curves;
indicates BIC curves.
From the same curves, k = 6 has been selected for specific reasons: i) as mentioned, the decrease after k = 6 is smoother, less steep; ii) for higher values of k, some coefficients for the interictal model may assume a magnitude such that |ϕ| > 1, which theoretically is not adequate due to the stationarity assumption [30]; iii) the main goal is to keep both models as simple as possible. It is worth mentioning that, for both conditions, the overall MSE is less than 0.05%, indicating a good representation of both spiking activities.
Fig 4 also shows the models obtained for both conditions considering k = 6. Overall, the AR models capture well the dynamics involved in the signals, from steadier behaviors, as the fluctuations around the baselines close to zero, to the abrupt changes that arise as spikes occur. Further results concerning the autocorrelation functions of the residuals are presented in Fig A in S1 Appendix, where similar conclusions can be drawn. Furthermore, the six coefficients found are different between the conditions. This result is important because it provides evidence that relatively simple models, as autoregressive ones, might be an interesting option to describe and distinguish types of seizures, thus providing patient-specific models with their own features.
Model extension
Once the AR(k) has been obtained considering Nmin, the same equations can be applied to the remaining samples of the whole time series (N > Nmin). Fig 5 presents both AR models fitted to their respective signals (interictal and ictal), as well as their MSE over time. To obtain this error, equally spaced moving windows of size Nmin were taken as references so that MSE could be calculated for each of them.
indicates the filtered downsampled signals;
indicates the fitted models;
indicates the mean values of the MSE. Each cell contains Nmin samples.
The errors exhibit a quite regular behavior, whose means are below 0.01%, the value that was previously found during order selection. This result is important especially envisioning more realistic practical situations, where storing an infinite number of samples is not feasible. Therefore, epochs of interictal or ictal signatures can be acquired to be used as templates for the AR models, and then used for the upcoming brain activities.
Discrete/continuous state-space and design of the observer/controller
Fig 6 presents the matrices for the discrete-time models (Eq 6), obtained directly from the AR coefficients, and the continuous-time models (Eq 9), after the appropriate conversion. For illustration purposes, only the m = 1, 2, M models are shown, i.e., the ones obtained from the two first time windows and the last one. Light colors indicate that the magnitude of the element inside a matrix increases positively, whereas cold colors indicate this magnitude increases negatively.
(A) and (B) represent and
, from Eqs 6 and 9, respectively. The set of matrices
, m = 1, 2, …, M, is used to compute single gains for the observer (L) and controller (G) in one of the two conditions (II or PIS). For illustration purposes, only the m = 1, 2, M models are shown, i.e., the ones obtained from the two first time windows and the last one.
All m = 1, 2, …, M together in the continuous form are used to design the gains of both observer and controller for one of the proposed schemes: hybrid (LII and GII, Fig 6A) and non-hybrid (LPIS and GPIS, Fig 6B). In this way, single gains assure dynamic stability when applied to the whole signal. Note that their structure is quite similar to each other, regardless of the condition, mostly having magnitudes of the type |ϕ|<1 (for discrete models), which is desired in the case of AR modeling [30].
For clarity of understanding, Table 1 shows the AR coefficients for some of the aforementioned matrices. There are differences from the coefficients calculated in both conditions, which indicates that AR models can capture the dynamic transitions inherent to the real signals, which is reflected on their coefficients. This distinction is particularly evidenced by the fluctuation in the values between each condition, especially for II, where the first term (ϕ1) is less stable across the windows presented (characterizing the irregular pattern of II compared to a more regular PIS activity [39]). This characteristic is specific to this type of modeling and has been used in pattern recognition approaches involving, for instance, structural health monitoring [47–50] as well as analysis and validation of computational models in describing II and PIS activities [39], since the AR coefficients are sensitive to changes in the signal dynamics [48].
Means are calculated considering all of the time windows (m = 1, 2, 3, …, M), not only the ones presented. Each line corresponds to the first line of the matrices (discrete case) from Fig 6. The table is restricted for illustration purposes. For further results, please see Fig B in S1 Appendix.
In fact, Ref [39] used the same signals present in this work, and by applying nonlinear/nonparametric classification techniques, such as logistic regression and confidence intervals, fairly good classifications were obtained. Moreover, within each condition, the coefficients are relatively bounded over the time windows (for further results, please see Fig B in S1 Appendix), which suggests that not only the patterns of II/PIS are sustained differently, but also that the coefficients do reflect this consistency.
Simulation
This section comprises the results from the simulation of the action of the observer and controller proposed in this work. A fourth-order Runge-Kutta was applied as the numerical integrator, considering the time step equal to the one from the real signal (that is: dt = 0.001s). This procedure is implemented on the Spyder platform, using Python 3.2, which is free and open-source.
Applying the observer/controller
Fig 7 presents the results obtained after the application of both observer and controller. For the sake of understanding, it was divided into three parts showing the first 5 seconds of simulation (and a zoomed in pannel during the first 0.1s): for the first one (Fig 7A), only the observer is activated while the controller is kept off. To verify the efficiency of the observer, the copy of the system is set with different initial conditions from the ones measured in the real system. It is expected that the gain L weighs the relative diffence between states such that both time signatures converge. This is evidenced in the same by the fact that the error e(t)→0 in a relative short time (first 0.1 seconds).
In (A), only the observer is activated while the controller is not. In (B), both are activated while setting xh = 0. In (C), both are activated while setting xh ≠ 0. For the observer, x(0) (states of the real system) and (states of the copy of the system) are set with different values and, as it progressively weighs the outputs
and y, the trajectories over time converge within a limited period of time. When xh = 0, the controller attenuates the epileptic spikes, but attempts to mitigate them towards zero. When xh ≠ 0, in turn, it delivers an additional input to simultaneously attenuate these spikes while driving the system to the II condition. In this case LPIS and GPIS are used.
In this case, the numerical integrator reconstructs the stored signal computationally over time. In practice, however, the process of signal acquisition would be carried out from a real electrode, positioned on the brain region of interest. In other words, by designing a state observer, the real PIS/II activities can be reconstructed (or “observed”) artificially. It must be stressed that this process is different from simply recording the signals: if there were more states in vector , all of them would be recreated, regardless of their physical interpration. Besides, this approach enables having an artificial dynamic model where control inputs can be applied.
Fig 7B presents the case in which both observer and controller are activated from the beginning of the simulation, while setting xh = 0. This implies that the control attempts only to drive the system to a constant baseline around zero. As shown, its action is efficient in attenuating the amplitudes of the spikes over time. Nonetheless, this is not the most realistic scenario, since it is desired not only to mitigate the PIS activity, but also drive the system back to the condition it was before a seizure started [13]. Thus, the reference state must necessarily be of the type xh ≠ 0, as presented in Fig 7C. Note that not only are the amplitudes of the spikes mitigated (from the PIS activity), but the controller progressively delivers the reference input to bring the signal back to the II condition.
At last, considering the normalized Euclidean norms between the reference state xh and actual states x(t), values fluctuate close to zero but do not reach this level in both applications of the controller. This is also expected, firstly because for xh = 0, the difference is calculated from the signal and zero, and since a complex activity is going on, a zero baseline is never achieved. Similarly, for xh ≠ 0 the spikes are attenuated but not canceled, and since it is desired that x → xh, the progressive addition of the reference xh to the controller makes it interact with the PIS activity reconstructed from the observer. The final result is not an exact version of the reference, but a close one. Therefore, it is expected that the controlled signal presents a similar behavior to the reference globally, not locally.
Fig 8 presents the action of both observer and controller, but now considering the whole signal from SUB. The observer remains activated all the time while the controller is only turned on at the middle of the signal (half of the entire time of simulation). In the time domain, note that there is again a clear difference in the observed signal as soon as the controller is activated. The amplitude of the undesired activity (i.e., PIS) is significantly reduced and the controller keeps delivering the input signal to drive it back to the II condition (xh), which is indicated by the irregular spiking pattern of higher amplitude that arises.
The controller is activated at half of the entire time of the signal, whereas the observer remains on during all the time. In this case LPIS and GPIS are used. Before the controller is activated, xref is the PIS signal being observed; after it is activated, xref is the non-seizure II signal.
However, such a difference is not reflected on the Euclidean norm, as previously discussed, resulting in no significant improvement. By plotting the spectrogram it becomes clear that the uncontrolled signal has a similar frequency content pattern as the PIS signal, whereas the controlled one presents higher similarity with the II signal (Fig 3). This is the global similar behavior expected when the controller is activated.
As part of the hybrid strategy also covered in this work, Fig 9 presents similar results, but using LII and GII, that is, both observer and controller are designed using a SUB signal in the II condition. A different signal from xh is used to match the amplitude of the epileptiform activity without significantly increasing the gains of the controller. Futher details can be found in Fig D in S1 Appendix. When the controller is activated, the amplitude of the undesired activity is mitigated and the peaks of the reference signal can be seen standing out among the others, as previously. This indicates that the design can be carried out without the need of having epileptiform activity recorded, only the pre-ictal states. This is an interesting result aiming at practical applications, where seizures often occur spontaneously and unpredictably − which makes it difficult to measure them several times.
The controller is activated at half of the entire time of the signal, whereas the observer remains on during all the time. In this case LII and GII are used. Before the controller is activated, xref is the PIS signal being observed; after it is activated, xref is the non-seizure II signal.
In this work, for illustration purposes, the controller remained activated for half of the duration of the signal. In approaches where the system’s equations are deterministic, such as the one applied to Epileptor [13], the control can be activated for a few seconds to suppress a seizure. In this work, only the signals are known, that is, the partial representation in time of a broader system, so control inputs that last longer must be applied to better visualize their effects. This is also discussed in detail in the next sections. In practice, it is not expected to act for a long time since this might not be the ideal case considering that patients would have to wear devices from which the inputs are delivered, and also for economy reasons, because if the controller is activated for short periods time, not continuously, its lifespan is increased. Besides, its intensity must also be taken into account for the same reasons. In this sense, the amplitude of such inputs must be adequately tuned so as not to compromise the health of desired brain activities. Fig 10 presents the responses of the signals in time when the parameters of the controller (α and γ) are tuned to decrease the intensity of the input (for both hybrid and non-hybrid cases). For illustration purposes, only the first 20 seconds of activation (and zoomed in windows) are shown. For the non-hybrid controller (Fig 10A), when the amplitude is limited it changes from its previous highest peak (based on the time window analyzed) at ∼1480mV to 740mV (in modulus), a reduction of 50%. Naturally, this implies that amplitudes of the undesired activity are not as efficiently mitigated as before, but a considerable attenuation can be seen, specially in the zoomed in epochs of the signal.
(A) contains results for the non-hybrid case, whereas (B) contains the results for the hybrid case: represents the reference PIS signal (from SUB);
represents the controlled signal before tuning the parameters;
represents the controlled signal after tuning them. The input forces are indicated as:
represents the control input without tuning and
represents the same input after tuning. In this case LPIS and GPIS and LII and GII are used, respectively. The parameters before tuning are: α = 1680, 1/γ2 = 1 and 1/γ3 = 1.55 (non-hybrid), and α = 810, 1/γ2 = 1 and 1/γ3 = 2 (hybrid). The parameters after tuning are: α = 840 and 1/γ2 = 1 and 1/γ3 = 1.55 (non-hybrid), and α = 675 and 1/γ2 = 1 and 1/γ3 = 1.5 (hybrid). For all cases, 1/γ1 = 1/10.
Concerning the shape of the control stimulus, two components can be distinguished: the first one acts in the opposite sense of the PIS activity, presenting spikes of inverse magnitude; the second one, in turn, is responsible for delivering spikes of irregular time intervals (similar to the II activity). This finding is consistent with previous observations made in the literature, where continuous periodic stimulation/trains of pulses [51–53] or randomly-spaced ones [54] have been relatively efficient in suppressing seizures. Thus, an interesting investigation in practice would be the use of combined types of stimulation to verify their controlling effects in real epileptiform activity. As for the inputs of the hybrid controller (Fig 10B), the main difference compared to the first one is that its amplitude decreases from ∼850mV to 595mV (a reduction of 30% in modulus). For both of them, a significant attenuation is observed nonetheless.
Comparison
To better assess the efficiency of the controller and endorse the previously presented results of Euclidean norm and spectrograms, this section presents results for the comparisons using PSDs, cross-correlations and PCA. For the sake of understanding: reference II represents the stored II signal (xh); reference PIS represents the stored PIS signal to be observed/controlled; the uncontrolled signal represents the observed PIS signal after activating the observer; the controlled one, in turn, represents the signal after activating the controller.
Power spectral density (PSD) and principal component analysis (PCA)
Figs 11 and 12 summarize all of these metrics for the non-hybrid and hybrid controllers, respectively. To generate the PSDs, several Nmin-long windows were considered such that they are taken as random variables, which makes it possible to calculate mean PSDs with standard deviations and confidence intervals (Figs 11A and 12A). By analyzing the means (Figs 11B and 12B), two interesting patterns can be observed concerning the controlled signal compared to the uncontrolled one: there is a slight decrease in power up to ∼10Hz and an increase after this frequency, consistent with the fact that II represents more complex activities, thus covering a large frequency band, while PIS is a more deterministic periodic activity [39]. This suggests that upon activation of the controller, the signal tends to move towards the reference II and leave the PIS condition, regardless of the control strategy.
(A) contains the PSD results: represent the PSDs obtained for each of the Nmin-long windows;
represents the mean PSD of the reference/observed II signal;
represents the mean PSD of the controlled signal;
represents the mean PSD of the PIS reference/observed signal;
represents the mean PSD of the controlled signal; … represents the confidence intervals at a 95% level. (B) represents all the mean PSDs together for visual inspection; (C) contains the maximum cross-correlation values between the combinations reference II/PIS and uncontrolled/controlled signals using all of the Nmin-long windows and their respective confidence intervals at a 95% level. (D) and (E) present the principal components of PCA highlighting the clusters individually and after considering II/controlled and PIS/uncontrolled groups.
(A) contains the PSD results: represent the PSDs obtained for each of the Nmin-long windows;
represents the mean PSD of the reference/observed II signal;
represents the mean PSD of the controlled signal;
represents the mean PSD of the PIS reference/observed signal;
represents the mean PSD of the controlled signal; … represents the confidence intervals at a 95% level. (B) represents all the mean PSDs together for visual inspection; (C) contains the maximum cross-correlation values between the combinations reference II/PIS and uncontrolled/controlled signals using all of the Nmin-long windows and their respective confidence intervals at a 95% level. (D) and (E) present the principal components of PCA highlighting the clusters individually and after considering II/controlled and PIS/uncontrolled groups.
Figs 11C and 12C show the results for the analysis of maximum cross-correlation coefficients. Each of the Nmin-long window was taken as a sample; then, uncontrolled and controlled signals were both compared to references II and PIS. As expected, PIS activities correlate better with uncontrolled ones, whereas II activities correlate better with controlled ones, which again evidences the effectiveness of the controller, but analyzed in the time domain − for this metric, the non-hybrid controller has a better performance, with ∼0.74 vs ∼0.48. In the 3-dimensional principal component space, when each of the Nmin-long PSD windows is projected four different clusters are revealed, corresponding to references II and PIS and uncontrolled/controlled signals. When plotted together, regardless of the control strategy used, the samples from uncontrolled signals are spatially closer to PIS references, reference II is farther from the groups, and samples from controlled signals fall in between. When plotting only two categories (relating uncontrolled to PIS samples and controlled to II ones), the difference between these clusters becomes more pronounced.
To assess statistical similarity between these four clusters, a bilateral Kruskal-Wallis test was carried out considering only the second principal component (PC2), which proved to be the best one to distinguish them. This test yielded H = 137.64 and for the non-hybrid controller and H = 105.39 and
for the hybrid one, thus indicating the samples are statistically different. Next, the post-hoc Dunn’s test was also performed to evaluate individual differences, whose results are summarized in Table 2. For the non-hybrid controller, this test can efficiently distinguish reference PIS from reference II and uncontrolled from controlled signals. Besides, it does not distinguish between reference PIS and the uncontrolled signal, and between reference II and the controlled signal at 95%, which endorses the effect of the controller in restoring the signal to the II condition. For the hybrid controller, as the non-hybrid one, the test distinguishes well between references PIS and II as well as controlled from PIS references, endorsing the effectiveness of the controller in attenuating PIS. However, it does attest the similarity between controlled and reference II, demonstrating that its performance compared to the non-hybrid one is moderate.
The following notation is adopted: if p < 10−6.
Discussion
This work addressed a relevant challenge in the context of epileptic seizures, which is obtaining a data-based model that describes its signatures over time. System identification approaches involving epileptiform activity (or from interictal to ictal states) are not commonly found in the literature. An important implication of doing so is to provide models that are patient-specific or seizure-specific; that is, models that adept to each type of time signature or patient. From these models, observers and controllers can be designed for purposes of recreating the epileptiform dynamics and controlling it artificially.
From the time-frequency analysis, the total number of samples was reduced, which considerably relieves the computational burden without compromising the signals’ inherent dynamics. In addition, the main frequency components were identified and used as a reference from which to define the templates (several time windows over time) for the AR models with an appropriate fixed time length. This strategy enabled obtaining AR models of reduced order: in this work, AR(k = 6) provided a good trade-off between accuracy and complexity, resulting in low values of MSE and other relevant criteria, such as AIC and BIC. Naturally, as the order gets higher, models become increasingly closer to their reference signals. However, one of the main goals of the proposed approach is to make the models as simple as possible (specially envisioning the later stages of state-space representation), besides avoiding potential problems of overfitting. In this sense, the models can sufficiently approximate the inherent dynamics of periodic ictal activity (PIS) and interictal activity (II), distinguishing them according to the AR coefficients obtained, thus endorsing previous results using the same signals [39].
For the model extension step, the AR(6) model calculated from the first Nmin samples was reapplied to the whole signal (N > Nmin), in both conditions (PIS and II). Since the evolution of the MSE over time was stable, this suggests that the identification from restricted time intervals can be used in other parts of the signal too, a fact that is also endorsed by the consistency of AR coefficients over time windows presented in Table 1 and Fig B in S1 Appendix. This result is particularly interesting because, in practice, obtaining arbitrarily long time signals of epileptiform activity, especially from humans, usually involves complex procedures, or recordings may not be readily available. Therefore, even with reduced information about the time signatures, it is possible to infer about their remaining dynamics over time, which is also important for the next steps of designing observers and controllers.
When applying the state observer, it was verified that the dynamic error e(t)→0 during the first instants of the simulation. This makes it possible to recreate the time signature over time accurately from stored signals, as the case of this work. This result is consistent with a previous study that reconstructed the unmeasurable activity of the Epileptor model based on its measurable activity, according to Brogin, J. A. et al. [Unpublished]. The main difference is that this time the reference signal is not dynamically simulated through the numerical integrator, but approximated based on observations of stored signals. In practice, it is expected that the observed signal is not a stored one, but brain activity measured in real time. As long as such activity can be properly recorded, the proposed approach can be applied.
As for the controllers, its activation could attenuate significantly the amplitude of the undesired activity, i.e., periodic ictal spiking (PIS), while simultaneously driving the system back to the interictal condition (II). To do so, it is required to record II events previously so as to provide a dataset from which to choose xh. This is an interesting result because, since seizures present a wide variety, the reference for each case and design is seizure/patient-specific as well. The design of hybrid controllers (state-space plant built using II signals applied to PIS signals) proved to be relatively effective, which means that it is still possible to design such a hybrid control in the absence of ictal events, a good advantage considering that seizures cannot be easily measured and often occur abruptly.
The metrics used to evaluate how close controlled signals are to their reference desired ones provided complementary results. The Euclidean norm could be reasonably applied during the first seconds of simulation to visualize the efficiency of the controller along with the dynamic error; however, for longer time windows, the norm did not prove to be the best measure for this purpose. In other words, the norm is an interesting tool to assess small-scale events, not global ones. On the other hand, both spectrograms and PSDs provided relevant information regarding the frequency content of signals: they revealed that PIS activities have higher power in lower frequency bands (∼10Hz) whereas II activities spread over the spectrum. Upon the activation of the controller, the controlled signal tends to recover a significant part of power lost in higher bands, which is also reflected on the similarity between both II and controlled signal’s spectrograms.
Concerning the normalized cross-correlation, both controllers provided significant results in the sense that PIS activities correlate better with uncontrolled signals and II activities correlate better with controlled ones. This effect is, however, stronger for the non-hybrid controller. At last, PCA confirmed that both controllers yield easily separable clusters in the principal component space, and they can be distinguished statistically.
Control inputs could be tuned and constrained, which is essential in practical applications not to deliver unnecessarily strong stimuli to the brain region of interest. The shape of the input stimulus is also an important result in this work. As seen, previous approaches considering trains of pulses or continuous stimulation for purposes of seizure suppression are commonly found in the literature [51–53]. Alternative approaches considering irregular time intervals for the simulation can also be found [54]. This raises an important question on which type of simulation is the most effective. On the one hand, applying pulses of inverse magnitude to “cancel out” the spikes of PIS activities seems an acceptable premise, which may explain the similarity of the results obtained here with the first studies. On the other hand, the reasoning behind aperiodic stimulation lies on attempting to disengage the hypersynchronism and reverberating effects of neural networks [54], which is also promising considering epilepsy as a progressive hypersynchronization process [1] and consistent with the definition of xh.
From a physical point of view, epilepsy can be interpreted as the progressive loss of complexity due to a hypersynchronized activity [1, 39], so considering this context and the results obtained here for the control inputs, it may be inferred that periodic stimulation acts on reducing the on-going undesired PIS activity, whereas aperiodic stimulation (x(t)→xh) acts on progressively re-establishing the complexity of the brain region that had been oversynchronized. In this sense, combined types of stimuli remains an open topic for future investigation.
Main contributions and advantages
Advances in system identification and control of seizures: the fact of providing a state-space representation that accurately describes epileptiform activity is a major step towards representing real brain activities. Previous dynamic models designed considering tools from nonlinear systems theory, such as Epileptor and Epileptor-2 [8, 9], found important patterns/classifications in epileptic seizures, and experimental observations have been consistent with some predictions. However high their fidelity and representativeness may be, these models remain fairly abstract to a certain extent. In the present work, the proposed approach is different because it obtains a mathematical set of equations directly from real brain activity. Furthermore, rearranging them into state-space notation enables the application of observers and controllers for seizure suppression purposes. It must also be considered that the signals used in this work were obtained from slices of the hippocampal tissues resected from human patients [39]. Therefore, controlling the activity of these signals is highly relevant in the context of epilepsy. A possible extension of this study is to analyze and build experimental set-ups where such brain tissues can be adequately stimulated and verify the outcomes. AR modeling and variations usually work for any kind of time series, so this step can be carried out using animal tissue, not necessarily human, which considerably facilitates experimental verification.
Seizure suppression on different subfields: by using observers and controllers, the signal studied in this work proved to be controllable in the sense that its spikings were attenuated while it was driven back to the reference II activity. Nonetheless, as seen in S1 Appendix, different signals (from subfields of the hippocampus, such as DG, CA1-CA4 and SUB) were also tested in this work and similar results were obtained. This demonstrates the flexibility and robustness of the proposed control approach regardless of the region of interest. Considering this, future strategies focused on controlling seizures locally can be carried out (on a smaller scale, specific to each subfield). Moreover, it may also be possible to design controllers taking into account time signatures from several subfields at the same time, thus providing a more global suppression approach. These are interesting topics for future investigation.
Possible applications in seizure classification and prediction: once the autoregressive models are available in practice, they can be used, for instance, to identity each type of seizure, or even for distinguishing ictal from interictal states. This is possible because their coefficients are sensitive to the inherent dynamics existing in real signals, and then will likely indicate a difference in behavior if it exists. Besides, given that the main characteristic of AR models is trying to foresee next states based on the time history, it is a predictor by design, a powerful and helpful feature to be taken into account. As a matter of fact, seizure prediction is a widely studied topic in neural engineering, but the methods employed for this purpose usually rely on machine learning [55–61]. Although regarded as a black-box method too, AR modeling for system identification in the context of seizure prediction still needs to be better explored, and may be an alternative to such algorithms.
Main limitations and disadvantages
The issue of stationarity and AR order: considering that AR models are most suitable to stationary time series, the baseline shifts represented by the spikings from the real brain activities analyzed may influence their quality, which helps explain why the PIS activity is not completely attenuated: since the modeling consists of a regression of the signal onto itself, features of its inherent dynamics are expressed by the AR coefficients; if such a dynamics is not fully captured during the fitting process, the resulting state-space representation may not be accurate enough. Nevertheless, consistent representations of the real signals have been obtained at a relatively low order, which indicates that, although simple, AR models may be an alternative for system identification in the context of epilepsy. In this sense, alternative autoregressive techniques can be used to enhance the present results, such as ARMA, ARIMA or NARMA [62–66]. In particular, the latter one is a nonlinear version of AR and might thus be appropriate to capture PIS/II dynamics efficiently.
The issue of noise and background activity: another interesting issue that must be taken into account is the level of noise and background activity in the signals. The proposed approach was partially successful due to a certain level of noise: it grants the signals a more stochastic/stationary nature, thus compensating for local abrupt variations. However, it may also mask important dynamic structures that would otherwise be more evident and meaningful for the identification and design procedures. In this work, filters were mostly applied to remove low-frequency components, so a great part of the remaining dynamics contains seizure-unrelated activity. On the one hand, this is necessary envisioning realistic applications; on the other hand, this might render the observers/controllers less efficient. Besides, this work assumes that the existence of such low-frequency components is negligible, as previously carried out [39], but they may actually lead to modifications in the future designs too.
The issue of controlling signals, not systems: since this work considers a reconstruction of signals that were previously measured and stored, the present strategy is not based on a full system control approach. In other words, there is no access to a more comprehensive mathematical/dynamical model that represents epilepsy in detail, such as the Epileptor [8], for example. Rather, control inputs were applied to data-driven models that were observed in time, which makes the present approach signal-specific (or partially patient-specific), but not generalizable. In other words, the controller is effective in attenuating the epileptiform in signals (or, in the present context, regions of slices from hippocampal subfields), which are partial representations of a broader and intricate dynamics of the whole brain. For this reason, while the controller is activated, it attenuates the undesired activity; as soon as it is deactivated, this effect cannot be visualized anymore. Therefore, the exact duration of a control input is yet to be determined in the case of real activity. This issue remains a challenge in the current literature, but considering the promising results obtained here, data-driven modeling may be a relevant resource in this direction.
Final remarks
This work proposed an efficient approach for attenuating the epileptiform activity of electrophysiological recordings considering the artificial reconstruction of the stored signals based on AR modeling and an integration between state observers and controllers. Several stages/steps of design, application and comparison are involved, thus characterizing a systematic and straightforward procedure that has not been reported elsewhere.
The results provided evidence that attenuating the undesired effects of real seizures can be computationally emulated, being efficient for different hippocampal subfields (DG, CA1, CA2, CA3, CA4 and SUB). Both non-hybrid and hybrid approaches can be used for this purpose: the first one provides stronger input stimuli to mitigate periodic ictal spikings while simultaneaously driving the signal to the interictal-like behavior; the second one, despite causing slightly milder attenuation effects, is also able to generate effective outcomes. This endorses the possibility of designing controllers in the absence of sufficient PIS recordings, which is considered an advantage in this work given the inherent difficulty in acquiring signals with seizure events.
The shape of the control input seems consistent with previous approaches that adopt regular and irregular stimulation individually for seizure suppression. A major contribution in this sense found in this work is that they may be more effective if combined, applied as a superposed one consisting of periodic pulses of inverse magnitude (similar to the regular stimuli) and chaotically-behaved ones (closer to the irregular stimuli). These findings may be of great interest for real applications comprising the modeling and description of seizures (characterized by the AR models), and seizure mitigation and suppression.
Supporting information
S1 Appendix. Complementary results including the autocorrelation function of the errors, AR coefficients and attenuation of epileptiform activity in DG, CA1 and CA2.
https://doi.org/10.1371/journal.pone.0298762.s001
(PDF)
Acknowledgments
To prof. Carla A. Scorza, Prof. Ricardo S. Centeno and Prof. Elza M. T. Yacubian, who contributed to the execution of previous phases of the project and the possibility of obtaining the electrophysiological recordings analyzed.
References
- 1. Iasemidis LD. Epileptic seizure prediction and control. IEEE Transactions on Biomedical Engineering. 2003;50(5):549–558. pmid:12769431
- 2. Dua T, De Boer HM, Prilipko LL, Saxena S. Epilepsy care in the world: results of an ILAE/IBE/WHO global campaign against epilepsy. Epilepsia. 2006;47(7):1225–1231. pmid:16886987
- 3. Surges R, Thijs RD, Tan HL, Sander JW. Sudden unexpected death in epilepsy: risk factors and potential pathomechanisms. Nature Reviews Neurology. 2009;5(9):492. pmid:19668244
- 4. Ryvlin P, Nashef L, Tomson T. Prevention of sudden unexpected death in epilepsy: a realistic goal? Epilepsia. 2012;54:23–28.
- 5. Fisher RS, Schachter SC. The postictal state: a neglected entity in the management of epilepsy. Epilepsy & Behavior. 2000;1(1):52–59.
- 6. Zhu D, Bieger J, Molina GG, Aarts RM. A survey of stimulation methods used in SSVEP-based BCIs. Computational Intelligence and Neuroscience. 2010;1–13. pmid:20224799
- 7. Vezzani A, French J, Bartfai T, Baram TZ. The role of inflammation in epilepsy. Nature Reviews Neurology. 2011;7(1):31. pmid:21135885
- 8. Jirsa VK, Stacey WC, Quilichini PP, Ivanov AI, Bernard C. On the nature of seizure dynamics. Brain. 2014;137(8):2210–2230. pmid:24919973
- 9. Chizhov AV, Zefirov AV, Amakhin DV, Smirnova EY, Zaitsev AV. Minimal model of interictal and ictal discharges epileptor-2, PLoS Computational Biology. 2018;14(5):1–25.
- 10. Proix T, Bartolomei F, Guye M, Jirsa VK. Individual brain structure and modeling predict seizure propagation. Brain. 2017;140(3): 641–654. pmid:28364550
- 11. Zhang H, Xiao P. Seizure dynamics of coupled oscillators with epileptor field model. International Journal of Bifurcation and Chaos. 2018;28(03):1850041.
- 12. Nagaraj V, Lamperski A, Netoff T. I. Seizure control in a computational model using a reinforcement learning stimulation paradigm. International Journal of Neural Systems. 2017;27(07):1750012. pmid:28030999
- 13. Brogin JAF, Faber J, Bueno DD. An Efficient Approach to Define the Input Stimuli to Suppress Epileptic Seizures Described by the Epileptor Model, International Journal of Neural Systems. 2020;30(11):2050062–2050062. pmid:32938259
- 14. Hindmarsh JL, Rose RM. A model of neuronal bursting using three coupled first order differential equations. Proceedings of the Royal society of London. 1984;221(1222):87–102. pmid:6144106
- 15. Izhikevich EM. Simple model of spiking neurons. IEEE Transactions on Neural Networks. 2003;14(6):1569–1572. pmid:18244602
- 16.
Ljung L. (1999). System identification. Wiley encyclopedia of electrical and electronics engineering. 1–19.
- 17. Noël JP, Kerschen G. Nonlinear system identification in structural dynamics: 10 more years of progress. Mechanical Systems and Signal Processing. 2017;83:2–35.
- 18.
Van den Bos A. (2007). Parameter estimation for scientists and engineers. New Jersey: John Wiley & Sons.
- 19. Duun-Henriksen AK, Schmidt S, Røge RM, Møller JB, Nørgaard K, Jørgensen JB, et al. Model identification using stochastic differential equation grey-box models in diabetes. Journal of Diabetes Science and Technology. 2013;7(2):431–440. pmid:23567002
- 20. Hashemi M, Vattikonda AN, Sip V, Guye M, Bartolomei F, Woodman MM, et al. The Bayesian Virtual Epileptic Patient: A probabilistic framework designed to infer the spatial map of epileptogenicity in a personalized large-scale brain model of epilepsy spread. NeuroImage. 2020;217:116839. pmid:32387625
- 21. Franaszczuk PJ, Bergey GK. An autoregressive method for the measurement of synchronization of interictal and ictal EEG signals. Biological cybernetics. 1999;81(1):3–9. pmid:10434388
- 22. Chisci L, Mavino A, Perferi G, Sciandrone M, Anile C, Colicchio G, et al. Real-time epileptic seizure prediction using AR models and support vector machines. IEEE Transactions on Biomedical Engineering. 2010;57(5):1124–1132. pmid:20172805
- 23.
Yu PN, Naiini SA, Heck CN, Liu CY, Song D, Berger TW. (2016, August). A sparse Laguerre-Volterra autoregressive model for seizure prediction in temporal lobe epilepsy. In 2016 38th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC) (pp. 1664–1667). IEEE.
- 24. Ouyang CS, Yang RC, Chiang CT, Wu RC, Lin LC. EEG autoregressive modeling analysis: A diagnostic tool for patients with epilepsy without epileptiform discharges. Clinical Neurophysiology. 2020;131(8):1902–1908. pmid:32599273
- 25. Attia A, Moussaoui A, Chahir Y. Epileptic seizures identification with autoregressive model and firefly optimization based classification. Evolving Systems. 2021;12(3):827–836.
- 26. Schiff ND, Victor JD, Canel A, Labar DR. Characteristic nonlinearities of the 3/s ictal electroencephalogram identified by nonlinear autoregressive analysis. Biological Cybernetics. 1995;72(6):519–526. pmid:7612723
- 27. Schiff ND, Victor JD, Canel A. Nonlinear autoregressive analysis of the 3/s ictal electroencephalogram: Implications for underlying dynamics. Biological cybernetics. 1995;72(6):527–532. pmid:7612724
- 28.
Ljung, L. (1998). System identification: theory for the user. Prentice hall.
- 29.
Durbin J, Koopman SJ. (2012). Time series analysis by state space methods. Oxford university press.
- 30.
Fox J. (2015). Applied regression analysis and generalized linear models. Sage Publications.
- 31. Shieh LS, Wang H, Yates RE. Discrete-continuous model conversion. Applied Mathematical Modelling. 1980;4(6):449–455.
- 32.
Åström KJ, Wittenmark B. (2013). Computer-controlled systems: theory and design. Courier Corporation.
- 33.
Tanaka K, Wang HO. (2001). Fuzzy control systems design and analysis. John Wiley & Sons Ltd.
- 34.
Ogata K. (2010). Modern control engineering (Vol. 5). Upper Saddle River, NJ: Prentice hall.
- 35. Fisher RS, Acevedo C, Arzimanoglou A, Bogacz A, Cross JH, Elger CE, et al. ILAE official report: a practical clinical definition of epilepsy. Epilepsia. 2014;55(4):475–482. pmid:24730690
- 36. Kalitzin SN, Velis DN, da Silva FHL. Stimulation-based anticipation and control of state transitions in the epileptic brain. Epilepsy & Behavior. 2010;17(3):310–323. pmid:20163993
- 37. Kalitzin S, Koppert M, Petkov G, Velis D, da Silva FL. Computational model prospective on the observation of proictal states in epileptic neuronal systems. Epilepsy & Behavior. 2011;22:S102–S109. pmid:22078510
- 38. Proix T, Bartolomei F, Chauvel P, Bernard C, Jirsa VK. Permittivity coupling across brain regions determines seizure recruitment in partial epilepsy. Journal of Neuroscience. 2014;34(45):15009–15021. pmid:25378166
- 39. Araújo NS, Reyes-Garcia SZ, Brogin JAF, Bueno DD, Cavalheiro EA, Scorza CA, et al. Chaotic and Stochastic Dynamics of Epileptiform-Like Activities in Sclerotic Hippocampus Resected from Patients with Pharmacoresistant Epilepsy. PlOS Computational Biology. 2022;18(4):e1010027 pmid:35417449
- 40.
Gordon RA. (2015). Regression analysis for the social sciences. Routledge.
- 41. Brogin JAF, Faber J, Bueno DD. Burster reconstruction considering unmeasurable variables in the Epileptor model. Neural Computation. 2021;33(12):3288–3333. pmid:34710900
- 42.
Slotine JJE, Li W. (1991). Applied nonlinear control (Vol. 199, No. 1). Englewood Cliffs, NJ: Prentice hall.
- 43. Reyhanoglu M. Dynamics and control of a class of underactuated mechanical systems. IEEE Transactions on Automatic Control. 1999;44(09):1663–1671.
- 44.
Lofberg J. YALMIP: A toolbox for modeling and optimization in MATLAB. in: 2004 IEEE international conference on robotics and automation, eds. IEEE, (2011), pp. 284–289.
- 45.
Sagnol G. Picos documentation. A Python interface to conic optimization solvers. Release 0.1, 1. 2012.
- 46. Bueno DD, Sandoval Goes LC, Gonçalves PJP. Control of limit cycle oscillation in a three degrees of freedom airfoil section using fuzzy Takagi-Sugeno modeling. Shock and Vibration. 2014.
- 47. Yao R, Pakzad SN. Autoregressive statistical pattern recognition algorithms for damage detection in civil structures. Mechanical Systems and Signal Processing. 2012;31:355–368.
- 48. Entezami A, Shariatmadar H. An unsupervised learning approach by novel damage indices in structural health monitoring for damage localization and quantification. Structural Health Monitoring. 2018;17(2):325–345.
- 49. Yang JH, Lam HF. An innovative Bayesian system identification method using autoregressive model. Mechanical Systems and Signal Processing. 2019;133:106289.
- 50. Yano MO, Villani LG, da Silva S, Figueiredo E. Autoregressive model extrapolation using cubic splines for damage progression analysis. Journal of the Brazilian Society of Mechanical Sciences and Engineering. 2021;43(1):1–15.
- 51. Velasco F, Velasco M, Velasco AL, Menez D, Rocha L. Electrical stimulation for epilepsy: stimulation of hippocampal foci. Stereotactic and functional neurosurgery. 2011;77(1-4):223–227.
- 52. Tellez-Zenteno JF, McLachlan RS, Parrent A, Kubu CS, Wiebe S. Hippocampal electrical stimulation in mesial temporal lobe epilepsy. Neurology. 2006;66(10):1490–1494. pmid:16554495
- 53. Velasco AL, Velasco F, Velasco M, Trejo D, Castro G, Carrillo-Ruiz JD. Electrical stimulation of the hippocampal epileptic foci for seizure control: a double-blind, long-term follow-up study. Epilepsia. 2007;48(10):1895–1903. pmid:17634064
- 54. Cota VR, de Castro Medeiros D, da Páscoa Vilela MRS, Doretto MC, Moraes MFD. Distinct patterns of electrical stimulation of the basolateral amygdala influence pentylenetetrazole seizure outcome. Epilepsy & Behavior. 2009;14(1):26–31 pmid:18824246
- 55. Faust O, Acharya UR, Adeli H, Adeli A. Wavelet-based EEG processing for computer-aided seizure detection and epilepsy diagnosis. Seizure. 2015;26:56–64. pmid:25799903
- 56. Acharya UR, Oh SL, Hagiwara Y, Tan JH, Adeli H. Deep convolutional neural network for the automated detection and diagnosis of seizure using EEG signals. Computers in Biology and Medicine. 2018;100:270–278. pmid:28974302
- 57. Yuan Q, Zhou W, Xu F, Leng Y, Wei D. Epileptic EEG Identification via LBP Operators on Wavelet Coefficients. International Journal of Neural Systems. 2018;28(8):1850010. pmid:29665725
- 58. Li Y, Cui W, Luo M, Li K, Wang L. Epileptic seizure detection based on time-frequency images of EEG signals using Gaussian mixture model and gray level co-occurrence matrix features. International Journal of Neural Systems. 2018;28(07):1850003. pmid:29607682
- 59. Sharma P, Khan YU, Farooq O, Tripathi M, Adeli H. A wavelet-statistical features approach for nonconvulsive seizure detection. Clinical EEG and Neuroscience. 2018;45(04):274–284.
- 60. Yuan S, Zhou W, Chen L. Epileptic seizure prediction using diffusion distance and bayesian linear discriminate analysis on intracranial EEG. International Journal of Neural Systems. 2018;28(01):1750043. pmid:28950741
- 61. Ansari AH, Cherian PJ, Caicedo A, Naulaers G, De Vos M, Van Huffel S. Neonatal seizure detection using deep convolutional neural networks. International Journal of Neural Systems. 2019;29(04):1850011. pmid:29747532
- 62. Landau ID, Karimi A. A recursive algorithm for ARMAX model identification in closed loop. IEEE Transactions on Automatic Control. 1999;44(4):840–843.
- 63. Jansson M. Subspace identification and ARX modeling. IFAC Proceedings Volumes. 2003;36(16):1585–1590.
- 64. Baldacchino T, Anderson SR, Kadirkamanathan V. Computational system identification for Bayesian NARMAX modelling. Automatica. 2013;49(9):2641–2651.
- 65. Liu H, Zhu L, Pan Z, Bai F, Liu Y, Liu Y, et al. ARMAX-based transfer function model identification using wide-area measurement for adaptive and coordinated damping control. IEEE Transactions on Smart Grid. 2015;8(3):1105–1115.
- 66. Retes PFL, Aguirre LA. NARMAX model identification using a randomised approach. International Journal of Modelling, Identification and Control. 2019;31(3):205–216.