Epileptic seizure suppression: A computational approach for identification and control using real data

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.

Fig 1. Illustration of the time windows used (whose minimum size is N_min) to obtain the AR models and to which the observers and controllers are applied, ranging from m = 1, 2, …, M.

For each time window, ${\mathbf{A}}_m^{\left\{c\right\}}$ 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.

Fig 2. General flowchart with the 3 major stages of the proposed SI and control approach.

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 (x_h) 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 ($\widehat{\mathbf{y}}$) 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 ($\widehat{\mathbf{y}}$) is compared to the reference behavior (x_h) 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.

Fig 3. Raw and filtered downsampled signals with their respective spectrograms.

(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.

Fig 4. Reference II and PIS signals from SUB and their respective AR models obtained after order selection based on the AIC/BIC/MSE criteria.

(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: $-\square -$ indicates AIC curves; $-\star -$ indicates BIC curves.

Fig 5. Model extensions calculated for the whole time series based on the AR models found using the first N_min samples.

$-$ indicates the filtered downsampled signals; $--$ indicates the fitted models; $\dots$ indicates the mean values of the MSE. Each cell contains N_min samples.

Fig 6. Discrete and continuous dynamic matrices obtained from the AR coefficients in state-space notation.

(A) and (B) represent ${\mathbf{A}}_m^{\left\{d\right\}}$ and ${\mathbf{A}}_m^{\left\{c\right\}}$, from Eqs 6 and 9, respectively. The set of matrices ${\mathbf{A}}_m^{\left\{c\right\}}$, 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.

Fig 7. Action of the observer and controller over time and behavior of the error dynamics and Euclidean norms.

In (A), only the observer is activated while the controller is not. In (B), both are activated while setting x_h = 0. In (C), both are activated while setting x_h ≠ 0. For the observer, x(0) (states of the real system) and $\widehat{\mathbf{x}}\left(0\right)$ (states of the copy of the system) are set with different values and, as it progressively weighs the outputs $\widehat{\mathbf{y}}$ and y, the trajectories over time converge within a limited period of time. When x_h = 0, the controller attenuates the epileptic spikes, but attempts to mitigate them towards zero. When x_h ≠ 0, in turn, it delivers an additional input to simultaneously attenuate these spikes while driving the system to the II condition. In this case L_PIS and G_PIS are used.

Fig 8. Action of the observer and controller over time analyzed in terms of the normalized Euclidean norm ∥x(t) − x_ref∥ and frequency content over time (spectrograms).

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 L_PIS and G_PIS are used. Before the controller is activated, x_ref is the PIS signal being observed; after it is activated, x_ref is the non-seizure II signal.

Fig 9. Action of the observer and controller over time analyzed in terms of the normalized Euclidean norm ∥x(t) − x_ref∥ and frequency content over time (spectrograms).

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 L_II and G_II are used. Before the controller is activated, x_ref is the PIS signal being observed; after it is activated, x_ref is the non-seizure II signal.

Fig 10. Effect of tuning α and B₁₁ on the input response of the signal in the time domain.

(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 L_PIS and G_PIS and L_II and G_II are used, respectively. The parameters before tuning are: α = 1680, 1/γ₂ = 1 and 1/γ₃ = 1.55 (non-hybrid), and α = 810, 1/γ₂ = 1 and 1/γ₃ = 2 (hybrid). The parameters after tuning are: α = 840 and 1/γ₂ = 1 and 1/γ₃ = 1.55 (non-hybrid), and α = 675 and 1/γ₂ = 1 and 1/γ₃ = 1.5 (hybrid). For all cases, 1/γ₁ = 1/10.

Fig 11. Comparison between the reference observed II and PIS signals and the uncontrolled/controlled ones using PSDs and PCA after activating the non-hybrid controller.

(A) contains the PSD results: $-$ represent the PSDs obtained for each of the N_min-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 N_min-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.

Fig 12. Comparison between the reference observed II and PIS signals and the uncontrolled/controlled ones using PSDs and PCA after activating the hybrid controller.

(A) contains the PSD results: $-$ represent the PSDs obtained for each of the N_min-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 N_min-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.