Model-based design for seizure control by stimulation

Objective. Current brain stimulation paradigms are largely empirical rather than theoretical. An opportunity exists to improve upon their modest effectiveness in closed-loop control strategies with the development of theoretically grounded, model-based designs. Approach. Inspired by this need, here we couple experimental data and mathematical modeling with a control-theoretic strategy for seizure termination. We begin by exercising a dynamical systems approach to model seizures (n = 94) recorded using intracranial EEG (iEEG) from 21 patients with medication-resistant, localization-related epilepsy. Main results. Although each patient’s seizures displayed unique spatial and temporal patterns, their evolution can be parsimoniously characterized by the same model form. Idiosyncracies of the model can inform individualized intervention strategies, specifically in iEEG samples with well-localized seizure onset zones. Temporal fluctuations in the spatial profiles of the oscillatory modes show that seizure onset marks a transition into a regime in which the underlying system supports prolonged rhythmic and focal activity. Based on these observations, we propose a control-theoretic strategy that aims to stabilize ictal activity using static output feedback for linear time-invariant switching systems. Finally, we demonstrate in silico that our proposed strategy allows us to dampen the emerging focal oscillatory sources using only a small set of electrodes. Significance. Our integrative study informs the development of modulation and control algorithms for neurostimulation that could improve the effectiveness of implantable, closed-loop anti-epileptic devices.

iEEG electrodes are commonly implanted in select patients undergoing presurgical assessment, where non-invasive neuroimaging methods failed to provide accurate localization of the epileptic network and the seizure onset zone. Although a wide array of rigorous clinical (e.g., seizure semiology) and neuroimaging (e.g., structural MRI and PET scans) assessments are used to guide the placement of iEEG electrodes, some cases remain in which the electrodes still fail to accurately localize the epileptic networks. Since the successful localization of the seizure onset zone and the epileptic network via iEEG is pertinent to our ability to model and ultimately control the underlying seizure dynamics, here we only focus only on seizure samples from patients labeled as "localizing" through clinical assessment (F.M.). Non-localizing iEEG datasets commonly exhibit patterns of seizure propagation marked by the presence of diffuse and slow seizure onset activity. The absence of the LVFA and sharp rhythmic spiking was the most important feature of the iEEG datasets with poor localization of the SOZ (see SI Figure 1). We provide a list of "localizing" patients, or those with a clear localization of the SOZ, in Table 1.

Changes in the distribution of eigenvectors following seizure onset
Focal seizures in the localizing iEEG datasets are marked by the emergence of epileptiform activity from the electrodes located in the vicinity of the seizure onset zone. As demonstrated in Figure 2 of the main text, eigenvectors of the system estimated in the peri-ictal regime can reflect these emerging focal oscillations. The changes in the maximum value and kurtosis of the eigenvectors following the seizure onset effectively highlight the shift from the disordered preictal regime to the highly ordered and focal ictal onset regime. We used these measures to demonstrate the similarity between focal seizure samples in our dataset.
Here we demonstrate that the seizure onsets are marked by a strong increase in the kurtosis and maximum value of the estimated eigenvectors associated with oscillations across a wide range of frequencies (see SI Figure 2). We tested the statistical significance of these observations by comparing the pre-and post-seizure onset distributions of the normalized kurtosis and maximum value of the eigenvectors estimated from all iEEG samples with well-localized seizure onset zones. As seen in SI- Figure 2, the normalized kurtosis and maximum values are significantly lower (Wilcoxon rank-sum test, p < 0.05) in a 5 sec preictal window compared to values calculated from a window with the same length immediately after both the earliest electrographic and unequivocal electrographic onsets across several frequency bands. These results indicate that the emergence of ictal sources is similarly reflected in the eigenvectors of the system, as focal oscillatory eigenmodes across a wide range frequencies.  Top panels show the average normalized kurtosis calculated from eigenvectors associated with different frequency bands at the earliest electrographic change (blue line) and at the unequivocal seizure onset (red line) across all localizing samples. The x-axis represents the time from seizure onsets. The shaded area represents the standard deviation of the samples. A non-parametric Wilcoxon rank-sum test reveals that seizure onset is marked by a wide-band significant increase (p < 0.05) in the normalized kurtosis values concatenated across all seizure samples over a 10-second peri-ictal time window. The color-coded curves in the bottom panels show the average normalized kurtosis calculated for each patient separately following the unequivocal seizure onset. (B) The top panels show the average normalized maximum values calculated from eigenvectors associated with different frequency bands at the earliest electrographic change (blue line) and at the unequivocal seizure onset (red line) across all localizing samples. The shaded area represents the standard deviation of the samples. Similar to panel A, we used a non-parametric Wilcoxon rank-sum test and observed that the seizure onset is marked by a wide-band significant increase (p < 0.05) in the normalized kurtosis values concatenated across all seizure samples over a 10-second peri-ictal time window. The color-coded curves in the bottom panels show the average normalized kurtosis calculated for each patient separately following the unequivocal seizure onsets.

Generalized pole placement method
The effect of direct electrical stimulation using the linear model can be captured by the following equation where ! ∈ ℝ ! is the input injected at time . The states of the sensors, i.e., the state of the system, can be also identified from the measured output of the system, i.e., cortical activity, and can be described as where ! ∈ ℝ ! correspond to the data collected at time by the intracranial electrodes. We propose using a static output feedback controller of the form and then the closed-loop dynamics is given by As a seizure intervention control strategy, we propose a closed-loop actuation, with the aim of driving the eigenvalues of the system closer to the center, towards the asymptotically stable zone and, thus, in effect damping the focal ictal onset dynamics. This problem can be studied as a particular case of the generalized switching pole placement problem studied in our prior work [1] and formulated as follows.
! Given a time-varying system (1) and (2), static output feedback described by (3), and a collection of closed subsets ! , … , ! ⊂ ℂ where we want the eigenvalues of the closed-loop system to be contained, determine such that is the -th eigenvalue of the closed-loop system described in (4).
Here, we introduce the terminology and basic notions used in the rest of this section for readers. For a given matrix ∈ ℂ !×! , the vectorization operator vec( ) ∈ ℂ !" consists of the columns of stacked below each other. We denote by !×! the zero matrix in ℂ !×! , and, given matrices ∈ ℂ !×! and ∈ ℂ !×! , the Kronecker product between and denoted by the × matrix ⊗ . Finally, we denote by Re( ) ∈ ℝ !×! and Im( ) ∈ ℝ !×! the real and imaginary parts of a matrix ∈ ℂ !×! . Now, let be an element and be a closed (not necessarily convex) set in a Hilbert space ℋ. Any ! ∈ such that ∥ − ! ∥≤∥ − ∥ for all ∈ is referred to as a projection of onto . In our case, we deal with finite-dimensional Hilbert spaces, where there is always at least one such point. It is known that if is a closed convex set then each point in ℋ has only one projection onto . Then, one can introduce the projection operator onto that is a function : ℋ → ℋ such that for each ∈ ℋ it returns the projection of onto , i.e., ( ).
We focus on a particular case of the following problem: given closed sets ! , … , ! in a finite-dimensional Hilbert space ℋ , determine a point in the intersection ! ! !!! (assuming this to be non-empty. Assuming all the sets are closed and convex, then this problem is solvable using the alternating projections method [2]. If the sets are nonconvex, then no global convergence guarantees exist, and also for different initial conditions the alternating projections method may not converge [3]. In our previous work [1], we show that solving problem ! is computationally hard (i.e., NP-hard) and consequently, it is unlikely that an algorithm exists that solves ! polynomially. Therefore in problem ′ ! , we propose a reformulation of ! as the problem of determining the gain that lies in the intersection of several sets.
Let ℒ ! be the set of all possible closed-loop matrices for the different evolution matrices ! , described by ℒ ! = ∈ ℝ !×! : = ! + for some , where = 1, … , , and a finite-time horizon, and let ℳ ! be the set of complex matrices with eigenvalues in the specified regions ! ! , … , ! ! given by Therefore, problem ! can be solved by addressing the following problem.
lies in the intersection of the sets ℒ ! that are convex sets, whose projection can be computed as follows.

Lemma 1 ([4]):
The projection of ∈ ℂ !×! onto ℒ ! is given by Nevertheless, ℳ ! is non-convex, and no easily computable projection scheme is available. Therefore, instead, we propose to use the following approximation ~ℳ ! (!) . Definition 1 ( [4]): Let X = VTV * be the Schur's decomposition, where V ∈ ℂ !×! is a unitary matrix and T ∈ ℂ !×! an upper triangular matrix. The approximate projection mapping ~ℳ ! ( ) of onto ℳ ! is given by otherwise , and * ∈ Σ, with Σ denoting the set of possible permutations, minimizes the following optimization function. ⋄ The mapping in Definition (1) is motivated by the fact that if is a symmetric matrix, then ~ℳ ! ( ) is the best approximation of in the Frobenius norm [4]. Further, the permutation * in Definition (1) can be determined by reducing the problem to a minimum weight maximum matching. This reduction process consists in determining permutation matrices ! and ! such that they minimize the trace( ! Θ ! ) where Θ is the matrix whose entry ( , ) contains the value is the projection of the complex number !! onto the closed convex set ! ! . Therefore, * can be described by the pairs ( , = * ( )) corresponding to the diagonal entries of ! Θ ! .
Next, in order to address ′ ! , in our previous work [1] we propose two variations of the alternating projections method, which involves sequentially computing the projection of a matrix in one of the sets, followed by finding the projection of the former in another set and so forth. Nonetheless, because of the projection on the sets ℳ ! is only approximated by ~ℳ ! , we propose to weight the quality of the projections in the alternating projection method. This quality is captured by a convex combination of the projection onto a convex set and the one onto a non-convex set. More specifically, we bias the approximation towards the projection onto the convex sets, whose projection can be exactly determined in a computationally efficient manner. For more details regarding the proposed algorithms see our previous publication [1].

Closed-loop feedback control using pole placement method
Seizure onset is marked by a transition from the disordered pre-ictal regime to the ordered regime following the ictal onset, marked by prolonged synchronized oscillations between several electrodes [5,6]. Although the ictal-onset oscillations persist for several seconds, system parameters estimated using a sliding-window fluctuate. Here we demonstrate that the calculated feedback gains for stabilizing the estimated system's eigenmodes (see manuscript Figure 3 for details) using the pole placement method [7] only guarantee results for a single window. SI- Figure 3 demonstrates that the eigenmodes of the closed system associated with seizure onset oscillations are unstable (|λ|>1), even after one second, when we use the feedback gains calculated from the seizure onset window (SI- Figure 3.B). Together, these results show the high sensitivity of the outcomes to small fluctuations in the estimated system parameters, when the feedback grains are calculated using pole placement method from a single window. Eigenvalues are sorted based on their frequencies, from highest to lowest. Blue numbered dots represent the empirically estimated values. We simulate the effect of closed-loop static feedback between the few chosen electrodes (marked by red in panel A), by representing the eigenvalues of the closed system (red numbered dots). The same five electrodes selected in Figure 4 of the main manuscript with the highest eigenvector loading at seizure onset (as seen in panel B) were selected as stimulating electrodes to mimic the limited channels of implantable neurostimulation devices. Nevertheless, unlike our proposed method, the pole placement method [7] requires all of the electrodes to also act as sensing electrodes. The static output feedback gains were calculated using the pole placement method with the control-theoretic objective of shifting all the higher frequency (>15 Hz) eigenvalues of the systems that were estimated from ten consecutive sliding widows following the seizure onset to the predefined zones represented by the green circle.      Figure 4 in the main text. The error bars represent the standard deviations around means. The red bars show the same values calculated from the closed system in Figure 4 in the main text. Non-parametric statistical testing (Wilcoxon rank-sum test, p < 0.05) shows that the loading of SOZ channels is significantly reduced at both time points. The green bar shows the same values calculated from the closed system in SI- Figure 3 using the pole placement (PP) method. Note that all high-frequency eigenmodes have low stability values, which explains the absence of the green bar at t = 0. However, the loading grows at the next one-second window and the closed system becomes unstable (as seen in SI- Figure 3). (B) Similar to panel A, bar plots represent the average loading of SOZ electrodes across higher frequency eigenvectors (>15 Hz) with high stability values (>0.6) for the seizure sample in Figure 5. Note that similar to panel A, the eigenvector loadings reduce significantly (Wilcoxon rank-sum test, p < 0.05) at both time points in the closed system.

A.
B. Fig. 2A. (A-D) The frequency (cyan trace) and stability (red trace) of eigenvalues associated with four representative eigenvectors -two associated with high frequencies and two associated with low frequencies -whose evolution is displayed in heatmaps.

Figure 5: Evolution of eigenmodes estimated from a second-order autoregressive model over the ictal period shown in
Note that the second-order model captures the emergence of local ictal sources following the seizure onset similar to the first-order model in Fig. 2   The elements marked by the red '*' show the instances with significant change from the period before to immediately after seizure onset (significance is evaluated by a non-parametric Wilcoxon rank-sum test, p < 0.05, corrected for multiple comparisons using Bonferroni correction). Note that changes in the kurtosis capture the emergence of ictal oscillations across different choices of spatiotemporal downsampling. These results also show that the modeled system is influenced by the sampling rate, the spatial coverage of seizure sources, and the total number of electrodes. B.

C.
All All All All All All All All All All All All

Number of Electrodes
Sampling Frequency Sampling Frequency Sampling Frequency Sampling Frequency show the instances with significant change from the period before to immediately after seizure onset (significance is evaluated by a non-parametric Wilcoxon rank-sum test, p < 0.05, corrected for multiple comparisons using Bonferroni correction). Note that changes in the maximum value capture the emergence of ictal oscillations across different choices of spatiotemporal downsampling. These results also show that the modeled system is influenced by the sampling rate, the spatial coverage of seizure sources, and the total number of electrodes. B.

C.
All All All All All All All All All All All All

Number of Electrodes
Sampling Frequency Sampling Frequency Sampling Frequency Sampling Frequency show the instances where no eigenmode was found in that frequency bands across all samples. Overall, these results show that the spectral profile of the modeled system, especially in the gamma band (25-55 Hz), is influenced by the sampling rate, the spatial coverage of seizure sources, and the total number of electrodes. (C) Panels show the difference between the results in panels B and A. The elements marked by the red '*' show the instances with significant change from the period before to immediately after seizure onset (significance is evaluated by a non-parametric Wilcoxon ranksum test, p < 0.05, corrected for multiple comparisons using Bonferroni correction). Note that although we identify significant changes in the average frequency of high-frequency oscillations (25-55 Hz), overall the eigenmodes' average frequency is not a sensitive marker of the emerging ictal oscillations as it is dependent on the spatiotemporal sampling. B.

C.
All All All All All All All All All All All All

Number of Electrodes
Sampling Frequency Sampling Frequency Sampling Frequency Sampling Frequency show the instances where no eigenmode was found in that frequency bands across all samples. Overall, these results demonstrate that the stability profile of the modeled system, especially in the gamma band (25-55 Hz), is influenced by the sampling rate, the spatial coverage of seizure sources, and the total number of electrodes. (C) Panels show the difference between the results in panels B and A. The elements marked by the red '*' show the instances with significant change from the period before to immediately after seizure onset (significance is evaluated by a non-parametric Wilcoxon ranksum test, p < 0.05, corrected for multiple comparisons using Bonferroni correction). Note that although we identify significant changes in the average frequency of high-frequency oscillations (25-55 Hz), overall the eigenmodes' average stability value is not a sensitive marker of the emerging ictal oscillations as it is dependent on the spatiotemporal sampling.  following the unequivocal seizure onset. Eigenvalues are sorted based on their frequencies, from highest to lowest. Blue numbered dots represent the empirically estimated values. We simulate the effect of closed-loop static feedback between a few electrodes (marked red in panel A), by representing eigenvalues of the closed system (red numbered dots). Only nine electrodes with the highest eigenvector loading at seizure onset were selected as sensing and stimulating electrodes to mimic the limited channels available in implantable neurostimulation devices. The static output feedback gains were calculated using the generalized pole placement method [1]. The control-theoretic objective was to shift all the higher frequency (>12 Hz) eigenvalues of the system, which were previously estimated from ten consecutive sliding windows (4-5s in length) following the seizure onset to the predefined zones (represented by green circles). Non-parametric statistical testing shows that the loadings of stimulating channels are significantly reduced at both time points (for n = 9 stimulating electrodes each of which has n=50 (open) or n=40 (closed) stable higher frequency eigenmodes, a Wilcoxon rank-sum test gives p = 2.5×10 !!" ; for n = 9 stimulating electrodes each of which has n-50 (open) or n=40 (closed) stable higher frequency eigenmodes, a Wilcoxon rank-sum test gives p = 2.9 ×10 !!" , respectively). Broadly, this numerical experiment demonstrates that our proposed algorithm can calculate static feedback that stabilizes the electrographic seizure activity in non-localizing seizure samples.

Choosing the length of the time window in the dynamical stability-based approach
Here, we investigate whether the accuracy of the dynamic stability approach depends on the time window over which we estimate the model. Intuitively we might imagine that if the time window is too large, we would not be able to accurately specify short timescale changes. Therefore, we applied dynamic stability analysis to the surrogate data windowed over a range of temporal lengths (SI - Figure 13). We found that the frequency and stability estimates did indeed depend on the length of the time window over which the model is fit. Specifically, the estimated frequency and damping of the oscillations were reduced progressively for larger window sizes, likely due to the presence of temporally overlapping patterns of oscillations during the time-windows considered for estimation.
Next, we examined the effect of window size on the modeled iEEG time series. We consider 1, 2, 5, and 10 sec time-windows while maintaining 100 msec shifts. The calculations provide similar intuitions to those that we obtained from the synthetic data: we observe a loss of sensitivity to high-frequency dynamics as the length of the timewindow increases (SI- Figure 14). These results provide evidence in support of smaller window sizes. Therefore, in our manuscript, we report the results obtained via 1 sec windows.

Figure 14: Effect of window size on the dynamical stability-based characterization of ictal periods. (A-B)
For the example seizure from Study 020 patient, we show the time evolution of the angles and stability of the eigenvalues associated with the different eigenmodes, for different window sizes (1, 2, 5, and 10 sec; shown across different columns). The red and cyan arrows highlight the onset and offset of the ictal period, respectively. (C-E) The temporal evolution of three representative eigenvectors: two associated with high frequencies and one associated with low frequencies. A.