Seizure pathways change on circadian and slower timescales in individual patients with focal epilepsy.

Personalized medicine requires that treatments adapt to not only the patient but also changing factors within each individual. Although epilepsy is a dynamic disorder characterized by pathological fluctuations in brain state, surprisingly little is known about whether and how seizures vary in the same patient. We quantitatively compared within-patient seizure network evolutions using intracranial electroencephalographic (iEEG) recordings of over 500 seizures from 31 patients with focal epilepsy (mean 16.5 seizures per patient). In all patients, we found variability in seizure paths through the space of possible network dynamics. Seizures with similar pathways tended to occur closer together in time, and a simple model suggested that seizure pathways change on circadian and/or slower timescales in the majority of patients. These temporal relationships occurred independent of whether the patient underwent antiepileptic medication reduction. Our results suggest that various modulatory processes, operating at different timescales, shape within-patient seizure evolutions, leading to variable seizure pathways that may require tailored treatment approaches.

Personalized medicine requires that treatments adapt to not only the patient but also changing factors within each individual. Although epilepsy is a dynamic disorder characterized by pathological fluctuations in brain state, surprisingly little is known about whether and how seizures vary in the same patient. We quantitatively compared within-patient seizure network evolutions using intracranial electroencephalographic (iEEG) recordings of over 500 seizures from 31 patients with focal epilepsy (mean 16.5 seizures per patient). In all patients, we found variability in seizure paths through the space of possible network dynamics. Seizures with similar pathways tended to occur closer together in time, and a simple model suggested that seizure pathways change on circadian and/or slower timescales in the majority of patients. These temporal relationships occurred independent of whether the patient underwent antiepileptic medication reduction. Our results suggest that various modulatory processes, operating at different timescales, shape within-patient seizure evolutions, leading to variable seizure pathways that may require tailored treatment approaches.
focal epilepsy | seizure dynamics | functional connectivity | within-patient variability | intracranial EEG F ocal epilepsy is characterized by spontaneous, recurrent seizures that arise from localized cortical sites (1). An unresolved question is how much seizures themselves can vary in individual patients. Past studies suggest that seizures within a single patient share common features (2)(3)(4)(5)(6) and evolve through a similar sequence (7), or characteristic pathway (8), of spatiotemporal neural dynamics. However, there is also evidence that various aspects of seizures can differ within the same patient. Clinically, some patients have multiple seizure onset sites that each produce their own characteristic seizure dynamics (9), and long-term electroencephalographic (EEG) recordings suggest that a subset of patients have multiple types of seizure evolutions (8,(10)(11)(12). Ictal onset patterns (13,14), the extent of seizure spread (15)(16)(17), and seizure recruitment patterns (18) can also differ in the same patient. This variability may arise from fluctuations in the underlying brain state (17,(19)(20)(21)(22)(23), suggesting that background neural activity affects not only seizure likelihood (20,24) but also seizure evolution. Crucially, a given treatment may only address a subset of a patient's seizures: for example, a single neurostimulation protocol may not control the complete repertoire of seizures (19), and a single prediction algorithm may fail to forecast all seizures (10,12,25). Consequently, seizure variability has important implications for clinical management in these patients.
To design optimal and comprehensive treatments, we therefore need to understand the prevalence and characteristics of within-patient seizure variability. Do seizure pathways vary in all patients? How are different seizure pathways distributed in time? To answer these questions, we must objectively quantify the similarity of seizure pathways. This task is challenging due to the complexity of seizure dynamics: a variety of spatiotemporal features change independently during seizure evolution. Although some studies have quantitatively compared withinpatient seizures (26)(27)(28)(29)(30)(31), the current gold standard remains visual inspection of ictal EEG by trained clinicians. This latter approach is time-consuming and subjective and can miss important features, including functional network interactions, that are difficult to detect visually.
Such functional network dynamics, also known as functional connectivity patterns, describe relationships between the activity recorded by different EEG channels. Temporal changes in network dynamics play important roles in seizure initiation, propagation, and termination (2,23,(32)(33)(34)(35)(36)(37)(38)(39)(40)(41), in part due to dynamic changes in the connectivity of the seizure onset zone (7,(42)(43)(44). Past work suggests that in some patients, the brain consistently progresses through a specific sequence of finite network states during seizures; however, other patients had unexplained variability in their seizure network evolutions (7). To fully understand how functional interactions support ictal processes, we must also understand if and how multiple seizure pathways, representing different ictal network evolutions, can coexist in an individual patient. Such diversity would reveal that the same Significance Epilepsy is an episodic disease characterized by brief periods of abnormal brain activity, known as seizures, that often have clinical correlates. In many patients, seizures preferentially happen during certain stages of daily and multiday cycles. However, it is unclear whether and how seizures themselves change over time, even though such variability may have clinical implications. To address this knowledge gap, we quantitatively analyze the nature of within-patient variability in seizure network evolutions. Contrary to common expectations, we find seizure variability throughout our cohort. Moreover, we demonstrate that seizures do not change randomly; instead, they also appear to fluctuate over daily and slower timescales. Ultimately, we may improve treatments by tailoring interventions to the full repertoire of seizures in each patient.
neural regions can variably interact to produce a variety of pathological dynamics.
In this paper, we therefore focus on quantifying and characterizing within-patient variability in seizure network evolutions; however, our approach can be adapted to compare the temporal evolutions of any ictal feature of interest. We first visualized and compared the within-patient seizure network evolutions of human patients with focal epilepsy (recorded for 43 to 382 h). In each patient, we analyzed all seizures with clear electrographic correlates (i.e., both subclinical and clinical seizures). In total, we quantitatively analyzed 511 seizures (average 16.5 seizures per patient), allowing us to characterize the nature of within-patient variability in these dynamics. In each patient, we found variability in seizure network evolution, revealing that within-patient seizures are not well represented by a single characteristic pathway through network space. However, seizures can share parts or all of the same pathway, with recurring dynamical elements across seizures. Furthermore, we explored how seizure pathways change over different timescales, providing insights into the temporal changes of within-patient seizures. Our analysis revealed that seizures change on circadian and/or slower timescales in each patient, suggesting that different modulatory processes shape seizure pathways.

Results
We analyzed seizure evolution in 31 human patients (511 seizures total, mean 16.5 seizures per patient) with focal epilepsy who underwent continuous intracranial electroencephalographic (iEEG) recordings as part of presurgical evaluation. Patient details are provided in SI Appendix, Text S1. In this section, we first demonstrate how we visualized seizure evolution through network space and quantified the dissimilarity of these seizure pathways. Importantly, our analysis captured differences in network interactions during seizures, which did not necessarily correspond to anatomical differences in the location and spread of seizure activity. We then describe the amount, form, and temporal patterns of within-patient seizure variability. Finally, we hypothesize how underlying processes occurring on different timescales could drive the observed changes in seizure pathways. We provide a visual guide to our approaches, along with the goals of each analysis, in SI Appendix, Text S2.

Visualizing and Quantifying Variability in Within-Patient Seizure
Pathways. Our first goal was to objectively compare within-patient seizure network evolution. For each patient, we extracted the seizure iEEGs (Fig. 1A) and computed the sliding-window functional connectivity, defined as band-averaged coherence in six frequency bands (Fig. 1B). Thus, each seizure time window was described by a set of six connectivity matrices that captured interactions between iEEG channels in each frequency band. We additionally normalized the magnitude of each connectivity matrix to focus on the evolving patterns of network interactions, rather than gross changes in the global level of coherence. The set of all possible connectivity patterns created a high-dimensional space, in which each location corresponded to a specific network configuration. As such, each time window could be represented by a high-dimensional data point, and the evolution of a seizure's network dynamics formed a pathway in this high-dimensional connectivity space. By transforming seizures in this manner, we framed our comparison of seizures as a comparison of seizure pathways (or trajectories) through the high-dimensional network space.
Due to the high dimensionality of this network space, it was infeasible to directly visualize seizure pathways. However, seizure pathways could be approximated in a two-dimensional (2D) projection using multidimensional scaling (MDS), a dimensionality reduction technique that attempts to maintain the distances between high-dimensional data points in the lower-dimensional space (Fig. 1C). This technique has been previously used to visualize ictal and interictal network dynamics (43). In our case, MDS placed seizure time windows in a 2D projection based on the similarity of their network configurations; each time window was represented by a single point, and points corresponding to time windows with more similar network dynamics were placed closer together. While imperfect, this approximation of the network space nonetheless provided an intuitive visualization for comparing seizure pathways in the same patient. For example, in patient 931, the projection demonstrated that two seizures may follow approximately the same pathway (seizures 6 and 8), part of the same pathway (seizures 8 and 9), or completely distinct pathways (seizures 2 and 10) through the network space, in agreement with visual impressions of the EEG.
To quantify these visual observations, we developed a seizure dissimilarity measure that provided a distance between two seizures based on their pathways through network space. Importantly, our approach recognized similarities in seizure pathways, even if the seizures evolved at different rates, by first applying dynamic time warping (45) to each pair of seizure functional connectivity time courses (SI Appendix, Text S3). Dynamic time warping nonlinearly stretches each time series such that similar points are aligned, thus minimizing the total distance between the two time series. We then defined the dissimilarity between two seizures as the average difference between the seizure pathways across all warped time points. The seizure dissimilarity matrix then summarized the dissimilarity between all pairs of seizure pathways in the same patient (Fig. 1D). In patient 931, seizures with similar pathways therefore had a low dissimilarity (e.g., seizures 6 and 8, dissimilarity 0.49); seizures with distinct, distant pathways had high dissimilarity (e.g., seizures 2 and 10, dissimilarity 3.21); and seizures with partially overlapping pathways had an intermediate level of dissimilarity (e.g., seizures 8 and 9, dissimilarity 1.75). Again, our measure of seizure dissimilarity agreed with intuitive comparisons of seizures based on visually assessing the iEEG (Fig. 1A) and MDS projections of the seizure pathways (Fig. 1C).
It is important to note that both seizure dissimilarity matrices and MDS projections were patient-specific: due to different electrode implantations, we could not compare seizures across patients using these network features. However, because we normalized the magnitude of the functional connectivity, we could compare seizure dissimilarity values across patients, even if the patients had different numbers of recording electrodes. In the remainder of the paper, we will focus on the across patients results, while using patient 931's seizures as examples. The seizure variability analysis of all patients is available on Zenodo (46) and summarized in SI Appendix, Text S4.
Seizure Variability Is a Common Feature in All Patients. Using our measure of seizure dissimilarity, we compared seizure pathways through network space in each patient. We first determined if seizure variability was present in all patients by visualizing the seizure dissimilarity matrix of each patient as a distribution of seizure dissimilarities (see Fig. 2A for an example). Note that in these distributions, each point corresponds to the difference in network evolutions of a pair of seizures, rather than a feature of a single seizure. Fig. 2B shows the distribution of seizure dissimilarities in each patient, with patients sorted from lowest (patient 934) to highest (patient I002 P006 D01) median dissimilarity. Although the average level of variability differed between patients (Fig. 2C), it is apparent that all patients had variability in seizure network evolutions. Even in patients with more consistent seizures, such as patient 934, there were pairs of seizures with high dissimilarity, indicating dissimilar seizure pathways. Many patients, including patient 931, had varying levels of differences between pathways, with only a few pairs of similar seizures. In all patients, network differences across all frequency bands contributed to the observed seizure dissimilarities, revealing that variability in seizure network evolutions was not limited to a narrow frequency range within a given patient (SI Appendix, Text S5). Additionally, we found that in the majority of patients, the observed variability was best described as a spectrum of seizure pathways, rather than distinct groupings of different seizure pathways (SI Appendix, Text S6). Thus, in most patients, the full diversity of seizure pathways could not be captured by a few archetypal seizures.
We also explored if the observed seizure variability was related to the available clinical information for each patient. First, we found that within the same patient, seizures of the same clinical type (subclinical, focal, or secondarily generalized) tended to be more similar than seizures of different clinical types; however, there was still a large amount of variability within a given seizure type (SI Appendix, Text S7). Thus, seizure variability in our patients was not solely explained by the presence of different clinical seizure types. This finding was expected given that seizures of the same clinical type may have different features in the same patient (16,47,48). Additionally, we found no association between postsurgical seizure freedom and measures of seizure variability (SI Appendix, Text S8). Likewise, higher levels of seizure variability were not associated with a particular seizure onset site (SI Appendix, Text S8). These findings suggest that the level of seizure variability is not associated with certain patient pathologies or treatment outcomes; instead, other factors may be more crucial for determining the extent and form of the variability.
Seizures with More Similar Pathways Tend to Occur Closer Together in Time. Many time-varying factors, such as sleep (22,24,47,49,50) and hormones (51)(52)(53)(54), are thought to influence seizure likelihood and dynamics. Additionally, during presurgical monitoring, antiepileptic medication is reduced in many patients, impacting brain dynamics (55). We therefore explored whether there is a temporal structure to how seizure pathways change over time in each patient. Fig. 3A shows the pathways of patient 931's seizures, as well as the time that each seizure occurred relative to the patient's first seizure. From this visualization, we see that the pathways gradually migrated through network space as the recording progressed, creating the observed spectrum of network evolutions. Moreover, looking at the seizure timings, we also see that seizures with similar pathways, such as seizures 6 to 8, tended to occur close together in time.
To quantify this temporal relationship, we defined a temporal distance matrix as the amount of time elapsed between each pair of the patient's seizures (Fig. 3B). Patient 931's seizure dissimilarity and temporal distance matrices have strikingly similar structures: groups of seizures with low dissimilarity tended to occur together in a relatively short time interval. In this patient, there was a strong and statistically significant positive correlation between these features (Spearman's ρ = 0.69, P = 0.001, one-tailed Mantel test), indicating that seizures with more similar pathways tended to occur closer together in time. Fig. 3C summarizes the relationship between seizure dissimilarities and temporal distances across all patients. In almost all patients, there was a positive Spearman's correlation between seizure dissimilarities and temporal distances (range, −0.10 to 0.83; mean, 0.45). This association was significant in 21 patients (67.7%) after false discovery rate correction. In these patients, we also observed that the average level of dissimilarity tended to increase with the time between the two seizures ( Fig. 3D). Interestingly, there was no association between whether antiepileptic medication was reduced and whether the correlation between seizure dissimilarities and temporal distances was significant (χ 2 test, P = 0.96) (SI Appendix, Text S9). Therefore, although medication levels may affect seizure occurrence and dynamics (9, 16, 56, 57), medication changes alone could not explain the observed shifts in seizure pathways, suggesting that other factors also play a role in shaping seizure features.
Seizure Pathways Change on Different Timescales. The observed temporal associations of seizure dissimilarities reflected gradual changes in seizure network evolutions across the length of each recording. In other words, we observed relatively slow shifts in seizure pathways over the course of multiple days. However, we also hypothesized that seizure pathways may change on shorter timescales due to, for example, circadian rhythms. Such rhythms would create timescale-dependent relationships between seizures; in particular, there would be a positive correlation between seizure dissimilarities and temporal distances on shorter timescales, but this association would be destroyed over longer timescales.
Therefore, to explore the possibility of different timescales of changes in seizure pathways, we scanned the correlation between seizure dissimilarities and temporal distances on different timescales T ranging from 6 h to the longest amount of time between a seizure pair (Fig. 4A). For example, for T = 3 d, we computed the correlation between seizure dissimilarities and temporal distances for all pairs of seizures that occurred within 3 d of each other. We refer to this set of correlations as a temporal correlation pattern of seizure pathways. Fig. 4A shows the temporal correlation pattern of patient 931's seizures. As we determined earlier, there was a positive correlation between seizure dissimilarities and temporal distances when all seizures were included in the computation (T = 5 d) as a result of the observed gradual changes in seizure pathways. At shorter timescales, however, the temporal relationship fluctuated; for example, the correlation was relatively low at T = 1 and 2.5 d and higher at T = 0.75 and 2.5 d. These fluctuations were signs of additional, timescaledependent changes in seizure pathways.
To investigate how these temporal correlation patterns arose, we modeled different patterns of seizure variability and the corresponding temporal correlation patterns (see Materials and Methods and SI Appendix, Text S10, for modeling details). For each patient, we then determined which pattern(s) of changes were most likely to reproduce the observed dynamics. In particular, we classified patients as having 1) linear changes in seizure pathways (Fig. 4B, Left), which corresponded to the slower, gradual shifts in seizure evolutions; 2) circadian changes (Fig. 4B, Middle), which represented dynamics modulated by daily rhythms; or 3) some combination of both the linear and circadian changes (example combination shown in Fig. 4B, Right). Fig. 4B demonstrates how seizure dissimilarities were simulated using patient 931's seizure times and example models from each of the above categories. In each model (Fig. 4B, Top), the x axis value of each seizure gives the seizure's time, relative to the first seizure. These values are the same across all three models because they are the empirically observed seizure times of patient 931. Thus, the x axis distance between a pair of seizures measures the amount of time, or temporal distance, between them. Based on the seizure times, each model then predicted how seizure pathways would change from seizure to seizure; specifically, the distance between two seizures along the other dimension(s) corresponds to the simulated dissimilarity of each pair of seizures. Each model additionally included noisy dynamics that allowed for further, random fluctuations in seizure pathways and thus seizure dissimilarities (SI Appendix, Fig. S10.1).
From these temporal distances and simulated seizure dissimilarities (Fig. 4B, Middle), we then computed the corresponding temporal correlation patterns (Fig. 4B, Bottom) using the same process shown in Fig. 4A. A linear change in seizure pathways produced a positive temporal relationship that was stronger at longer timescales. Meanwhile, a circadian model only produced strong, positive temporal correlations at timescales shorter than   1 d. Finally, a combination of the linear and circadian factors created both the short-term temporal relationships and a positive temporal correlation at the longer timescales. Note that there were also some additional fluctuations in the temporal correlation patterns due to noisy changes in dynamics; these effects differed depending on the outcome of the noisy simulation. To fully explore these noisy effects, we therefore additionally varied the level of noise added to the models. The tested combinations of noisy, linear, and circadian contributions are provided in SI Appendix, Table S10. For each combination of these factors, we simulated temporal correlation patterns 1,000 times using different noise realizations to produce a series of possible temporal correlation patterns for each model. Fig. 4C shows the type of model (linear, circadian, or linear + circadian) most likely to underlie the observed temporal correlation pattern of each patient. As a measure of model fit, we defined model likelihood as the percentage of model simulations that reproduced the patient's observed temporal correlation pattern. Model likelihood ranges from 0 to 100%, and higher values reveal that the modeled changes in dynamics consistently produced temporal correlation patterns that were similar to the patient's observed temporal correlation pattern. We additionally required the selected model to 1) outperform noisy simulations alone; 2) clearly distinguish between the linear and circadian models; and 3) in the case of the linear + circadian model, clearly outperform one of the simpler models. Using these criteria, 17 patients' temporal correlation patterns were best explained by the linear model, 3 by the circadian model, and 7 by the linear + circadian model. Thus, most patients (77.4%) required a linear component to explain the observed changes in seizure pathways, while 32.3% of the patients were well matched by a model incorporating circadian changes in pathways. Notably, model likelihood tended to be higher for patients with higher number of seizures, reflecting greater model certainty in cases with larger sample sizes (SI Appendix, Fig. S10.4). These different classifications of seizure variability were not associated with surgical outcomes (SI Appendix, Text S8) or whether the patient's medication was reduced during presurgical monitoring (SI Appendix, Text S9). Four patients' temporal correlation patterns could not be assigned to a model, either because the linear model and circadian model performed similarly (patient Study 038) or the best model did not outperform noise alone (patients Study 017, Study 033, and 1163). Additionally, in some patients (e.g., patients Study 020, 756, and 1196), only a small percentage of the simulations matched the observed temporal correlation patterns, indicating that reproducing the observed dynamics required specific patterns of noise. In these cases, other models may therefore provide a better explanation for the patient's changes in seizure pathways. In particular, many of these patients had strong positive correlations at timescales longer than 1 d but less than the length of the recording, suggesting multiday fluctuations in seizure pathways.

Discussion
We have quantitatively compared seizure network evolutions within individual human patients with focal epilepsy, revealing that seizure variability is a common feature across patients. We often observed pairs of seizures with relatively low dissimilarity due to their largely conserved pathways through the space of possible network dynamics, suggesting that seizure evolution is not purely random. However, we likewise found that a single dynamical pathway cannot comprehensively represent all of a patient's seizure evolutions. Interestingly, seizure pathways changed over time in most patients, with more similar seizures tending to occur closer together in time. Our modeling results indicate that in most patients, a combination of circadian and/or slower changes in seizure pathways may underlie the observed Linear seizure dissimilarity temporal distance  variability, suggesting that factors operating on different timescales modulate within-patient seizure evolutions. We investigated variability in seizure functional network evolution due to the importance of network interactions in ictal processes (2,7,23,(32)(33)(34)(35)(36)(37)(38)(39)(40)(42)(43)(44) and build on previous work by demonstrating within-patient variability in these pathological network dynamics. However, in future work, the framework we present could easily be adapted to compare other features that highlight different aspects of seizure dynamics. For example, a univariate feature that captures the amplitude and frequency of ictal discharges may be better suited for comparing the involvement of different channels, similar to how clinicians visually compare EEG traces. Data from other recording modalities, such as microelectrode arrays, could be analyzed to evaluate consistency in neuronal firing patterns between seizures (4, 5). Meanwhile, although we do not perform biophysical modeling of seizure dynamics in this work, other studies have used model inversion to hypothesize how the activities of different neuronal populations change during seizures (8,58,59). Comparing the parameter time courses of such models could reveal different patterns of changes in neural activity during a patient's seizures. Finally, due to patient-specific recording layouts, we focused on comparing seizure pathways within individual patients. However, comparing seizures across patients, either using spatially independent features or common recording layouts, in future studies could uncover common classes of pathological dynamics (8,60).
To quantify within-patient variability in seizure pathways, we developed a seizure dissimilarity measure that addresses the challenges of comparing diverse spatiotemporal patterns across seizures. A few previous studies have attempted to quantitatively compare seizure dynamics using either univariate (27,28,30,31) or network (26,29) features computed from scalp or intracranial EEG. These earlier dissimilarity measures were based on edit distance, which captures how many replacements, insertions, and deletions are required to transform one sequence into another. Importantly, unlike this previous method, our dynamic time warping approach recognizes that two seizures are equivalent if they follow the same pathway, even if they do so at different rates. Despite this difference, those past studies also reported both common and disparate dynamics across within-patient seizures; however, their analysis was limited to a small number of patients and/or seizures per patient. Our work provides insight into the prevalence and characteristics of seizure variability by analyzing over 500 seizures across 31 patients. Finally, we expand on previous work by using seizure dissimilarity to characterize temporal changes in seizure evolutions.
Previous work has found that within-patient seizures have similar dynamics (2)(3)(4)(5)(6)(7)(8), although variability may be introduced through different rates of evolution (4,61) or early termination in the seizure pathway (6,8). In our cohort, we observed that subsets of within-patient seizures follow approximately the same dynamical pathway through network space, and such similar groups of seizures likely underlie these past findings. However, we also found that the complete repertoire of within-patient seizure network evolutions was poorly characterized by a single, characteristic pathway. Notably, we also found that a patient with different seizure pathways did not necessarily have subsets of distinct pathways; instead, small variations between seizures often produced a spectrum of pathways. An intriguing possibility is that various decision points, existing on the framework of potential seizure pathways, produce a repertoire of seizure evolutions. Future studies are needed to map these potential seizure pathways and uncover the factors that determine how individual seizures evolve.
The crucial question is then how these different seizure pathways arise from the same neural substrate. In theory, a range of changes before or during the seizure can affect its network evolution. We hypothesize that spatiotemporal changes in the interictal neural state produce seizures with different characteristics. Past studies suggest that neural excitability (20,57,62), inhibition (61), and network interactions (23,63) influence certain spatiotemporal seizure features, such as the rate and extent of seizure propagation. These changes in brain state may be driven by various factors, including sleep (22,47,49), hormones (51)(52)(53)(54), and medication (55). If interictal dynamics indeed shape how seizures manifest, future research will need to determine how specific interictal features relate to seizure characteristics. One possibility is that elements of seizure networks are activated during interictal states (43); thus, seizures with different network features could be preceded by preictal periods with corresponding network structures. Researchers could also relate preictal/interictal networks to other seizure features, such as seizure propagation patterns, perhaps by investigating how the underlying structural connectome relates to functional networks (64) and mediates seizure spread (63). Importantly, the relationship between preictal network dynamics and seizure features could be limited to a specific frequency band (23), which could in turn suggest possible physiological mechanisms for the observed changes in seizure dynamics (65,66). Additional aspects of interictal dynamics, such as the pattern of highfrequency oscillations (21) and band power changes (17), may also be linked to changes in seizure features. Overall, a better understanding of both functional and pathological fluctuations in interictal dynamics could suggest ways in which the background brain state alters seizure evolutions.
Toward this goal, prolonged recordings of patients with focal epilepsy may provide insight into how pathological brain dynamics change over time and influence seizure features. In particular, recent studies using such data have shown that the rates of epileptiform discharges and seizures fluctuate according to both circadian and patient-specific multidien (approximately weekly to monthly) cycles (50,67). An intriguing possibility is that the same factors that rhythmically modulate seizure likelihood may also influence seizure evolution. Consistent with this hypothesis, we found that the majority of observed temporal patterns of seizure variability were well explained by models incorporating circadian and/or linear changes in seizure pathways. In particular, the linear component of the model may reflect gradual changes in pathways on slower timescales, ranging from weeks to months. These simple models provided an initial hypothesis for the observed patterns of changes in seizure evolutions. Some patients' seizure patterns may be better explained by more complex models that capture different dynamics, such as multistability or multidien cycles. Ultimately, it is likely that various factors, with differential effects on seizure evolution, interact to produce the observed repertoire of seizure pathways. Analyzing within-patient seizure variability in long-term recordings could provide additional insight into such patterns of changes in seizure pathways.
Many of the patients in our study underwent antiepileptic medication reduction as part of presurgical monitoring, making it difficult to disentangle the effects of changing drug levels from other potential slow-varying modulators of seizure pathways. Changes in antiepileptic medication can impact neural excitability (68)(69)(70), and medication tapering increases seizure likelihood in most patients (16,71); however, it is controversial whether it also affects seizure patterns (9,16,56,71). In some cases, it appears that medication tapering reveals latent seizure pathways that are suppressed by medication (9) or allows existing pathways to further progress (e.g., the secondary generalization of typically focal seizures) (16). It is possible that the impact of medication reduction on seizure dynamics is drug-, patient-, and dose-dependent and may ultimately depend on how well the medication controls neuronal excitability (57). However, medication changes alone cannot account for the observed seizure variability in our cohort, as we observed temporal associations of seizure pathways in patients that did not undergo medication reduction. In future work, associating medication levels with differences in seizure pathways could help untangle the different factors shaping seizure dynamics.
Another confounding factor in our data is that the surgical implantation itself could artificially alter seizure dynamics. Using chronic recordings of epileptic canines, Ung et al. (72) found variability in seizure onset and interictal burst dynamics, with the most stable dynamics emerging approximately a few weeks after electrode implantation. In agreement with their work, we found that earlier seizure evolutions often recurred later in the recording, making it unlikely that gradual changes in the recording quality or an acute reaction to the surgery underlay the observed variability. Instead, Ung et al. (72) hypothesized that seizure variability results from transient, atypical dynamics as the brain recovers from surgery, with later dynamics representing a truer epileptic network. Other stressors, such as medication withdrawal, could similarly elicit abnormal dynamics. Nevertheless, a large number of our patients had good surgical outcomes, suggesting that their recorded seizures accurately represented their epileptic networks. Additionally, clinicians often note that patients have typical seizures during iEEG recordings, as compared to preimplantation reports, despite the effects of surgery and medication withdrawal (16). As such, the observed seizure dynamics in our cohort may be part of their usual repertoires of seizure evolutions, even if some dynamics are only elicited by strong stressors. Further analysis in chronic human recordings, such as the NeuroVista dataset (8,12), is needed to determine whether and how seizure pathways vary in a more naturalistic setting.
Contrary to the expectation that high levels of seizure variability may worsen surgical outcomes, we found no association between these patient features. It may be that only some types of variability, such as multifocal (9) or secondarily generalized (73) seizures, impact the likelihood of seizure freedom following surgery. Importantly, variability in the seizure onset network state does not indicate that a patient has multifocal seizures, as different network configurations can be associated with the same apparent ictal onset zone. Additionally, variability in seizure pathways may not be inherently deleterious, as long as it is observed and accounted for when planning the surgical resection. Indeed, due to the short presurgical monitoring time and limited spatial coverage of the recording electrodes, some potential seizure pathways may not have been captured (11,72), leading us to underestimate the level of variability in some patients.
Although the amount of seizure variability was not associated with postsurgical seizure freedom, it may have implications for clinical treatments. First, regardless of the source of the observed seizure variability, the different seizure dynamics observed during presurgical monitoring provide crucial information for guiding surgical resection. For example, recent studies suggest that seizure network properties can help identify epileptogenic tissue (7,74,75); however, we must determine if seizures with different network evolutions provide equivalent localization information. Seizure variability may also have implications for seizure prediction. In particular, in that same patient, seizures with different pathways may have distinct preictal signatures, making seizure prediction more difficult (10,12). A successful seizure prediction algorithm would either need to recognize multiple signatures or find common features among the disparate preictal dynamics. Finally, neurostimulation offers a promising new approach for controlling seizures; however, in rodent models, the effectiveness of a given stimulation protocol depends on the preictal brain state (19). Thus, such interventions may need to recognize and adapt to the specific characteristics of each corresponding seizure evolution in order to control all seizures. Importantly, our cohort was limited to patients with medication refractory focal epilepsy who were candidates for surgical resection. The characteristics and clinical implications of seizure variability may be different in other patient cohorts.
More generally, our work adds to the growing literature on within-subject variability in brain dynamics and other physiological states (76) in both health and disease. In particular, there is increasing interest in developing improved, time-sensitive treatments that adapt to the patient's changing state. Rather than delivering a continuous or regular therapy, such treatments would be modified and/or timed to improve their efficacy while also reducing treatment side effects. Treatment parameters may be tuned to biological rhythms (50,67,(77)(78)(79) or respond directly to fluctuating conditions within each patient (80)(81)(82). To investigate temporal fluctuations within each patient and determine how treatments interact with these changes, researchers may draw inspiration from spatiotemporal analyses in other fields, such as ecology (83), genetics (84), and engineering (85,86), as well as develop new techniques that address specific dataanalytical challenges.
In summary, we have shown that there is within-patient variation in seizure network evolution in patients with focal epilepsy. Temporal changes in seizure evolution suggest that a combination of circadian and slow-varying factors shape these seizure pathways, perhaps by modulating the background brain state. Further research is needed to determine whether and how preictal dynamics influence seizure pathways. Uncovering these mechanisms could provide novel approaches for predicting and controlling seizures that are tailored to the complete repertoire of pathological neural dynamics in each patient.

Materials and Methods
Patient Selection and Data Acquisition. This work was a retrospective study that analyzed seizures from 13 patients from the Mayo Clinic and the Hospital of the University of Pennsylvania [available on the IEEG Portal, https://www. ieeg.org/ (87,88)] and 18 patients from the University College London Hospital (UCLH) who were diagnosed with refractory focal epilepsy and underwent presurgical monitoring. Patients were selected without reference to the cause or other characteristics of their pathology. All IEEG Portal patients gave consent to have their anonymized iEEG data publicly available on the International Epilepsy Electrophysiology Portal (https://www.ieeg. org/) (87,88). For the UCLH patients, their iEEG was anonymized and exported, and the anonymized data were subsequently analyzed in this study under the approval of the Newcastle University Ethics Committee (reference number 6887/2018).
For each patient, the placement of the intracranial electrodes was determined by the clinical team, independent of this study. Ictal segments were identified and extracted for the analysis based on clinical seizure markings. To be included in the study, each patient was required to have had at least six seizures suitable for the analysis. This threshold was chosen to allow examination of seizure variability in a broad cohort of patients, while still ensuring that enough seizures were observed to draw conclusions about the characteristics of seizure variability in each patient. Seizures were excluded from the analysis if they did not have clear electrographic correlates (with clear onset and termination), if they were triggered by/occurred during cortical stimulation, if they had noisy segments, or if they had large missing segments. Periods of status epilepticus and continuous epileptiform discharges were also excluded. However, electrographic seizures without clinical correlates (i.e., subclinical seizures) were included in the analysis as they may have either similar or disparate dynamics (relative to clinical seizures) that convey clinically relevant information (89). Additional information about each patient and the analyzed seizures is shown in SI Appendix, Text S1.
iEEG Preprocessing. For each patient, if different seizures were recorded at multiple sampling frequencies, all of the recordings were first downsampled to the lowest sampling frequency. Noisy channels were then removed based on visual inspection. In the remaining channels, short sections of missing values were linearly interpolated. These sections of missing values were <0.05 s with the exception of one segment in seizure 2 of patient Study 020, which was 0.514 s. All channels were rereferenced to a common average reference. Each channel's time series was then bandpass filtered from 1 to 150 Hz (fourth-order, zero-phase Butterworth filter). To remove line noise, the time series were additionally notch filtered (fourth-order, 2-Hz width, zero-phase Butterworth filter) at 60 and 120 Hz (IEEG Portal patients) or 50, 100, and 150 Hz (UCLH patients).
Computing Functional Connectivity. To compute the time-varying functional connectivity of each seizure, a 10-s sliding window, with 9-s overlap between consecutive windows, was applied to each preprocessed ictal time series. The same sliding window parameters have previously been used to estimate time-varying coherence in ictal iEEG data (90). For each window, the coherence between each pair of iEEG channels was computed in six different frequency bands (delta 1 to 4 Hz, theta 4 to 8 Hz, alpha 8 to 13 Hz, beta 13 to 30 Hz, gamma 30 to 80 Hz, and high gamma 80 to 150 Hz). The coherence in each frequency band was computed using band-averaged coherence, defined as where f 1 and f 2 are the lower and upper bounds of the frequency band; P i,j (f) is the cross-spectrum density of channels i and j; and P i,i (f) and P j,j (f) are the autospectrum densities of channels i and j, respectively. In each window, channel autospectra and cross-spectra were calculated using Welch's method (2-s sliding window with 1-s overlap). Note that band-averaged coherence is equivalent to coherence but filtered in the frequency domain to the frequency band of interest. As such, band-averaged coherence ranges from 0 to 1 and will be higher when the two signals have a consistent phase and amplitude relationship in the specified frequency band.
Thus, in a patient with n iEEG channels, the functional connectivity of each time window was described by six symmetric, nonnegative, n × n matrices, in which each entry (i,j) gives the coherence between channels i and j in the given frequency band. Each matrix was then written in vector form by rearranging the upper-triangular, off-diagonal elements into a single column vector of length (n 2n)/2. Each vector was normalized so that the L1 norm (i.e., sum of all elements) was 1, thus ensuring that differences between connectivity vectors captured a change in connectivity pattern rather than gross changes in global levels of coherence. This normalization step also allowed the magnitude of seizure dissimilarities to be compared across patients with different numbers of electrodes. For each time window, the six connectivity vectors were then vertically concatenated together, forming a single column vector of length 6*(n 2n)/2. Each patient's ictal connectivity vectors were subsequently horizontally concatenated together to form a matrix V containing 6*(n 2n)/2 features and m observations, where m is the total number of ictal windows across all seizures.
Dimensionality Reduction and Visualization. Small fluctuations in the functional connectivity due to noise would create a high baseline dissimilarity between seizures. Therefore, to reduce noise in the connectivity matrices, nonnegative matrix factorization (NMF) (91) was used to approximately factor each patient's ictal time-varying connectivity matrix V into two nonnegative matrices, W and H, such that V∼W × H (details provided in SI Appendix, Text S11). The matrix W contained patient-specific basis vectors, each of which had 6*(n 2n)/2 features that captured a pattern of connectivity across all channels and frequency bands. Each original ictal time window was summarized as an additive combination of these basis vectors, with the coefficients matrix H giving the contribution of each basis vector to each time window. These factorizations were patient-specific since the basis vector features depended on the iEEG electrode layout in each patient. The optimal number of basis vectors, r, was determined using stability NMF (84).
For each patient the selected factorization was then used to create V* = W × H, a lower-rank approximation of the original time-varying seizure functional connectivity (SI Appendix, Text S11). This return to the original feature space is necessary since NMF basis vectors are not orthogonal, and distances in NMF basis vector space are therefore not equivalent to distances in feature space. Each reconstructed connectivity vector was then renormalized to have an L1 norm of 1, ensuring that differences in reconstruction accuracy did not affect the distances between different ictal time points. To visualize the connectivity vectors of patient 931's seizures in Fig. 1C, all time seizures windows were projected into a 2D embedding using MDS (specifically, Sammon mapping) based on their L1 (cityblock) distances in the highdimensional reconstructed feature space.
Computing Seizure Dissimilarities. Following the NMF-based reconstruction of the seizure connectivity, the network evolution of each seizure was described by a multivariate time series with 6*(n 2n)/2 features. To compare network evolutions across within-patient seizures, a seizure dissimilarity matrix was created for each patient. Each pair of seizure functional connectivity time series was first warped using dynamic time warping, which stretches each time series such that the total distance between the two time series is minimized (SI Appendix, Text S3). This step ensures that 1) similar network dynamics of the two seizures are aligned and 2) the warped seizures are the same length. We chose to minimize the L1 distance between each pair of seizures as this metric provides a better measure of distances in high-dimensional spaces (92).
Following dynamic time warping, the L1 distance between the pair of warped time series was computed, resulting in a vector of distances capturing the dissimilarity in the seizures' network structures at each time point. The seizure dissimilarity between the two seizures was defined as the average distance across all warped time points. The seizure dissimilarity matrix contains the dissimilarities between all pairs of the patient's seizures. Note that seizure dissimilarity is not a metric distance, because the triangle equality does not necessarily hold; however, it performs similarly to alternative metric distances of seizure dissimilarity (SI Appendix, Text S12).
Comparison to Temporal Distances. For each patient, we computed a temporal distance matrix containing the amount of time elapsed (measured in days) between the onset times of each pair of seizures. Spearman's correlation was computed between the upper triangular elements of the seizure dissimilarity matrix and the temporal distance matrix of each patient. Since the distances in each matrix were not independent observations, the Mantel test (93) was used to determine the significance of each correlation. Briefly, the rows and columns of one matrix were randomly permuted 10,000 times. The correlation between the two sets of upper triangular elements was recomputed after each permutation, resulting in a distribution of correlation values that described the expected correlation if there were no relationship between seizure dissimilarities and temporal distances. The P value of the association was then defined as the proportion of permuted correlation that was greater than or equal to the observed correlation. To correct for multiple comparisons, the Benjamini-Hochberg false discovery rate (FDR) correction (94) was applied to the set of P values computed across all patients (31 total tests). The correlation was considered significant if the associated adjusted P value was less than 0.05.
Computing Temporal Correlation Patterns. To quantify how seizure dynamics change over different timescales in each patient, Spearman's correlation between seizure dissimilarities and temporal distances was computed only for seizure pairs with temporal distances less than or equal to timescale T. T was scanned from 0.25 d up to the patient's largest temporal distance in steps of 0.25 d. A timescale was excluded from the analysis if fewer than seven pairs of seizures occurred within the given timescale or if no new seizure pairs were added when the timescale was increased. The resulting set of correlations across various timescales was referred to as a temporal correlation pattern.
Modeling Seizure Dissimilarities and Temporal Correlation Patterns. To determine the underlying processes that could produce the observed temporal correlation patterns, changes in seizure dynamics were modeled using the functions f l ðtÞ ¼ 1 7 t ða line with a slope of one per weekÞ, f c ðtÞ ¼ sin 2 πtða sine wave with a period of 1 dÞ, f n ðtÞ Nð0,1Þ ðGaussian noise with a mean of zero and SD of 1Þ, where t is time in days. For each function, a simulated distance matrix D was then defined for the patients' seizures, with where t i is the time of seizure i, t j is the time of seizure j, and f(t) is the corresponding function. The dissimilarity of the two seizures was then defined as where l, c, and n are scalar parameters controlling the relative contributions of the linear, circadian, and noise functions, respectively. The relative contributions of the linear, circadian, and noise functions were scanned by varying the levels of l, c, and n. For each set of parameters, seizure dissimilarities were simulated 1,000 times using different noise realizations (and correspondingly changing the noise distance matrix, D n ), and the resulting temporal correlation patterns were computed for each set of simulated dissimilarities. Note that because temporal correlation patterns only depend on the order of the dissimilarities, only the relative magnitudes of the parameters l, c, and n affected the modeling results. A model was termed a linear model if c = 0, a circadian model if l = 0, and a linear + circadian model if l > 0 and c > 0.
To determine if a patient's seizure dynamics could be categorized as linear, circadian, or linear + circadian, the simulated temporal correlation patterns were compared to the patient's observed temporal correlation pattern by computing the mean squared error (MSE) of each simulated pattern. Simulated temporal correlation patterns with MSE ≤ 0.02185 were defined as good matches to the observed dynamics. This threshold was chosen because it was the fifth percentile of the set of all MSEs, across all patients, and based on visual inspection of simulated temporal correlation patterns with different MSEs. The likelihood L of a given parameter set was then defined as the percentage of good matches produced by the 1,000 noisy simulations of seizure dissimilarities at those parameter values. For each class of model (linear, circadian, or linear + circadian), the model's likelihood (L l , L c , or L l+c , respectively) was the highest likelihood among the model type's parameter sets, and the best model was the model with the highest likelihood. L n was also defined as the highest likelihood of the parameter sets without any linear or circadian contributions (l = 0, c = 0, n > 0). This best model with likelihood L max was then used to categorize the patient's dynamics if it outperformed all competing models. Specifically, we required that 1) the best model clearly outperform noise alone (L max ≥ 2L n ); otherwise, the patient's dynamics were classified as other/indeterminate; 2) the performance of the linear model and circadian model were clearly distinguishable (L l ≥ 2L c if the linear model was best; L c ≥ 2L l if the circadian model was best); otherwise, the patient's dynamics were classified as other/ indeterminate; and 3) if the best model was linear + circadian, it clearly outperform the two simpler models (L l+c ≥ 2L l and L l+c ≥ 2L c ); otherwise, the patient's dynamics were classified as the simpler model (if one simpler model performed comparably by this criterion) or as other/indeterminate (if both simpler models performed comparably).
See SI Appendix, Text S10, for additional modeling details and the selected models for each patient.
Code and Data Availability. All data were analyzed using MATLAB version R2018b. To perform NMF, we used the Nonnegative Matrix Factorization Algorithms Toolbox, available at https://github.com/kimjingu/nonnegfacmatlab/, which implements the alternating nonnegative least squares with block principal pivoting algorithm (95,96). For the remainder of the analysis, we used MATLAB implementations of standard algorithms (multidimensional scaling [Sammon mapping]: mdscale, criterion "Sammon"; dynamic time warping: dtw; hierarchical clustering: linkage; Torgerson's multidimensional scaling: cmdscale; gap statistic: evalclusters; FDR correction: mafdr) or custom code. The iEEG time series of all IEEG Portal patients is available at https://www.ieeg.org/. The NMF factorizations of all analyzed seizure network evolutions, along with the code for producing the primary downstream results (seizure dissimilarity matrices, clustering, and temporal analysis) and figures, are published on Zenodo (46).