Domino-like transient dynamics at seizure onset in epilepsy

The International League Against Epilepsy (ILAE) groups seizures into “focal”, “generalized” and “unknown” based on whether the seizure onset is confined to a brain region in one hemisphere, arises in several brain region simultaneously, or is not known, respectively. This separation fails to account for the rich diversity of clinically and experimentally observed spatiotemporal patterns of seizure onset and even less so for the properties of the brain networks generating them. We consider three different patterns of domino-like seizure onset in Idiopathic Generalized Epilepsy (IGE) and present a novel approach to classification of seizures. To understand how these patterns are generated on networks requires understanding of the relationship between intrinsic node dynamics and coupling between nodes in the presence of noise, which currently is unknown. We investigate this interplay here in the framework of domino-like recruitment across a network. In particular, we use a phenomenological model of seizure onset with heterogeneous coupling and node properties, and show that in combination they generate a range of domino-like onset patterns observed in the IGE seizures. We further explore the individual contribution of heterogeneous node dynamics and coupling by interpreting in-vitro experimental data in which the speed of onset can be chemically modulated. This work contributes to a better understanding of possible drivers for the spatiotemporal patterns observed at seizure onset and may ultimately contribute to a more personalized approach to classification of seizure types in clinical practice.


Introduction
Recurrent spontaneous seizures, the hallmark of epilepsy, are characterized by behavioral symptoms alongside abnormal patterns of electrical activity in the brain. These are believed to be caused by an imbalance between excitation and inhibition within neural populations leading to hyper-excitability at the macro scale [1].
The International League Against Epilepsy classification of epilepsy 2017 classifies seizures by type of onset and brain networks affected. A seizure initiated in a local brain region is classified as focal, while a seizure involving activation of brain regions on both sides of the brain is classified as generalized. In the case when the above is not clear or known, a seizure is classified as unknown [2,3]. This division represents a practical separation made on clinical grounds and (where possible) informed by electroencephalogram (EEG) recordings in order to guide treatment. By lumping together a variety of seizure onset patterns into these three broad categories, the diversity of spatiotemporal patterns observed at seizure onset, for example, focal onset of generalized seizures, are typically overlooked. Note that focal onset of a generalized seizure is simply an EEG variation seen in idiopathic generalized epilepsy (IGE) but not an indication of focal epilepsy [4]. However, ambiguity of seizure onset captured in the clinic could result in significant diagnostic delay as well as unnecessary invasive investigations and inappropriate anti-epileptic drug therapy in some cases.
Despite recent advances, there remains insufficient knowledge for a purely scientific classification of seizure types [2,3,5]. Current research identifies seizures as the initiation and recruitment of neuronal populations in distinct but interconnected brain regions. Brain activity can thus be modeled by means of a dynamic network in which each node represents a population of neurons in a specific region of the brain. Understanding the spatiotemporal dynamics driving recruitment of nodes from background to seizure state across a given network is crucial to inform and advance scientific taxonomy of seizure onsets. Local and global changes in functional and structural network properties generate "hyper-excitable" networks, characterized by high frequency switches from background to seizure states [6]. However, the relative contribution to the network level activity from individual or populations of nodes is not well understood.
In this paper, we identify and investigate a variety of domino-like onset patterns observed in EEG recordings in IGE and characterize them by total recruitment time (i.e. the time it takes for seizure activity detected in one electrode to spread to all electrodes) and the lag time between recruitment of sequential electrodes. We show that different patterns of detected activity are clinically observed in the same individual. We develop a framework for modeling seizure onset and use it to demonstrate how different onset patterns arise as a result of interplay between heterogeneous node dynamics and heterogeneous coupling between nodes. We then explore the individual contribution of heterogeneity in node excitability and coupling among nodes.
To this end, we use a phenomenological network model of seizure onset previously studied by us and others [7][8][9][10][11]. Each node in the network has two stable states: background (steady) and seizure (oscillatory). The transition between the two is driven by an external random 'noisy' input. The same model has been used to explore domino-like recruitment of nodes across motif networks [7,9]. Each node is set to the background state (domino standing upright on a table) and transits to the oscillatory state (domino lying down on a table) due to noisy perturbations and inputs from other nodes (the table is wobbled). Depending on how likely it is to fall over (intrinsic property) or how close to other dominoes it is (coupling) you get slow or fast toppling patterns. We show that the domino-like seizure onset patterns observed in the EEG data have fast-domino and slow-domino cascades. Moreover we extend this approach by introducing the novel multiple-domino onset where large lags (gaps) separate groups of sequential electrodes for which the detected activity has similar transition times.
We build on previous work by considering larger networks and introduce heterogeneity to the node excitability and coupling weights. We first apply the model to generalized seizures observed in EEG recordings that exhibit fast-, slow-and multi-domino onset patterns. To illustrate our findings, we present in detail the analysis of a single individual with Juvenile Absence Epilepsy, taken from the data set detailed in [4,12]. Secondly, to unambiguously separate the contribution of node excitability or coupling weight to the observed onset patterns requires a tractable and controllable in-vitro model. To this end, we also apply our modeling framework to the slow onset of seizure-like activity induced along slices of medial entorhinal cortex (mEC) from mice, detailed in [13]. Our results highlight the complex interplay between network properties and provide a compelling case for further investigation into their impact on emergent dynamics in order to better categorize seizures and advance patient-specific diagnosis and treatment.

Clinical EEG data
We use epileptiform events recorded using EEG from patients diagnosed with genetic generalized epilepsy (GGE). Full details of the database can be found in [4,12]; for completeness we summarize briefly. Ambulatory EEG (Compumedics Ltd, Melbourne, Australia) were recorded from 107 patients for a period of 24-hours using 32 gold cup electrodes arranged according to the international 10-20 system. Signals were recorded with a sampling rate of 256Hz. The ProFusion 4 software (Compumedics Ltd, Melbourne, Australia) was used by an experienced EEG reader to review the recordings, identify and extract epileptiform events. Patients were categorized into five groups depending on syndrome: childhood absence epilepsy (CAE); juvenile absence epilepsy (JAE); juvenile myoclonic epilepsy (LME); generalized epilepsy with generalized tonic-clonic seizures only (GTCSO); and genetic generalized epilepsy unspecified (GGEU). Data was manually annotated by the clinician based on the common average reference [4,12,14]. The data set consists of 15 second long epochs each containing an event positioned such that the clinician identified onset time is at 5 seconds. Events from recordings were classified into three types: generalized paroxysms (events lasting longer than 2s), generalized fragments and focal discharges. The onset of generalized paroxysms was further classified into generalized onset or focal onset by a trained EEG reader. Focal onset of a generalized paroxysm is defined as lead-in focal discharges lasting at least 0.2s [12].
In this paper, we consider epochs containing generalized paroxysms, that we will henceforth call seizures. We extract 1267 epochs from the database containing seizures from 59 patients, including 15 seizures classified as focal onset. Note that there are no epochs containing seizures in the GGEU group. We submit each 15s epoch to the seizure onset detection algorithm, detailed below. We illustrate the application of the modeling framework using three consecutive seizures from one person with JAE, which is the most common syndrome within the data set [12].

Mouse mEC recording
Detailed description of the procedures used to collect recordings and analyze the data from the mouse mEC are given in [13]. All procedures were carried out in accordance with the UK Animal (Scientific Procedures) Act 1986 and were approved by the University of Exeter Animal Welfare and Ethical Review Body. Here we provide a brief summary for completeness. Parasagittal brain slices (400 μm thick) containing mEC were collected from male C57/BL6 mice. The slices were transferred to an interface-style recording chamber, and allowed to equilibrate for 30 minutes before a bath application of 4-aminopyridine (4-AP; 100 μM concentration equilibrated). A silicone probe consisting of 16 individual shanks (55 μm wide, 100 μm apart), with a single electrode contact point at the end of each shank (Neuronexus, Ann Arbor, MI; probe catalog number: A16x1-2mm-100-177), was positioned along the dorsal-ventral axis of the mEC; see Fig 1. A 32-channel amplifier (RHD2132; Intan, Los Angeles, CA) coupled to an open-source acquisition board (Open Ephys Inc, Cambridge, MA) was used to take recordings and the data were subsequently band-pass filtered (1-500 Hz) and digitized at 2 kHz.
In this paper, we apply our modeling framework and analysis to six recordings of length between 18 to 50 min, each containing several 'seizure-like' events. From each recording we determine the first onset event in which all 16 channels transit to a seizure-like state; this occurs 9-20 min from the start of the recordings using the same threshold detection algorithm (detailed below) to determine onset in each channel.
In further experiments, the induction of epileptiform activity by bath application of 4-AP was followed by co-application of GABA A receptor modulators as either positive modulation via diazepam (30 μM) or negative with Ro19-4603 (10 nM) [15]. In a final set of experiments considered here, a scalpel blade was used to make a cut in the intermediate mEC, thus anatomically separating dorsal and ventral portions. In this case, data were recorded from two sites, one in the dorsal and one in the ventral end of the slice; the 16 channel multielectrode array was not used in this case [13].

Seizure onset detection
We refer to the moment of transition from background to seizure state as seizure onset time. It is important to note that we detect this change in activity in the same way in both the EEG and mEC data sets, using a threshold detection algorithm as in [13] implemented using custom code in MATLAB (ver. R2018b).
The threshold detection algorithm is as follows: 1. The signal from a single channel or electrode is band pass filtered then z-normalized followed by rectifying, giving the processed signal x(t) for each time point t.
2. An envelope u is determined using spline interpolation over local maxima separated by at least np data points.
3. Seizure onset time τ n for each electrode n is the first time over a threshold defined by s×standard deviation from the mean of the envelope, The choice of band-pass filter, np and s values is different in the two cases (human and animal) and depends on the sampling frequency and seizure duration time in each data case. The values used in this paper are given in Table 1 and the algorithm procedure is illustrated in S1 Fig. For the EEG data, we submit 1267 epochs (described above) to the threshold detection algorithm. We only consider epochs in which onset is detected in all electrodes for further analysis; onset is detected in all electrodes in 1159 seizures out of the total 1267. From the detected onset times for each seizure we calculate the recruitment time defined as t n = τ n − min(τ n ) for each electrode n. We order the electrodes for each seizure by recruitment time from min(t n ) = 0 to max(t n ). The total recruitment time is max(t n ) and we calculate the maximum lag by taking the largest difference between consecutive (ordered) recruitment times. Note here that we are interested in the diverse range of patterns observed at seizure onset and that the seizures are recorded using ambulatory EEG so may occur under different (unknown) conditions e.g. sleep/wake. Therefore, we do not average recruitment times over seizures and instead submit them to a linear discriminant analysis detailed in the next section.
In the case of mouse data we submit the six manually identified mEC epochs containing seizure-like events to the threshold detection algorithm. Onset is detected in all channels of the six events. The data and seizure onset times are shown in S2 Fig. Here, we average the recruitment times t n in each channel over all six events to get the mean recruitment timet n ¼ ht n i and cascade duration (the total mean recruitment time) maxðt n Þ À minðt n Þ. The averaging of recruitment times is representative of the system's dynamics as the in-vitro experimental set up is well controlled compared to the EEG recordings. The cascade duration is consistent for choices of s > 0.65 and np > 16000; see S3 Fig.

Seizure onset pattern grouping
We refer to the patterns of dynamic activity captured by the EEG electrodes and detected by the threshold detection algorithm as onset patterns. We first consider a subset of 27 generalized seizures: 12 seizures clinically classified as generalized onset from one individual and 15 seizures clinically classified as focal onset from multiple individuals. The onset patterns for these 27 events are shown in S4 and S5 Figs; note one focal onset seizure was excluded from further analysis as seizure onset was detected in only 18 out of 19 channels.
Based on initial visual inspection of the detected onset patterns we propose three groups: 1. Fast domino: short total recruitment time and no large lags between detected activity in consecutive electrodes; 2. Slow domino: long total recruitment time and no large lags between detected activity in consecutive electrodes; 3. Multiple domino: long total recruitment time and one or more large lags that separate electrodes into subgroups for which the detected activity has similar recruitment times.
This grouping is inspired by the slow and fast domino effect detailed in [7,9]. Both papers consider the escape of nodes from quiescent to oscillatory states in small motif networks. The fast domino effect is where all nodes escape in quick succession, almost simultaneously: this is our fast domino onset where electrodes are recruited in quick succession and the total recruitment time of all nodes is small. The slow domino effect is seen when all nodes on a network escape but it takes longer between consecutive nodes to escape giving a longer total escape time: this is our slow domino effect in which there are small lags between recruitment of electrodes but the total recruitment time for all electrodes is longer than the fast domino onset. The multiple-domino group is a novel class identified in this paper where subsets of electrodes escape in either a fast or slow domino effect, and there is one or more large lag in time between subsets. This gives the effect of multiple domino cascades that could only be observed in larger networks (> 3 nodes).
To quantitatively delineate between the groups we chose an initial threshold between fast and slow total recruitment time to be 0.5s and the threshold between long and short lags to be 0.4s. We then compute linear discriminants based on this data set using MATLAB. We use these delineations to assign each seizure, including the original 26, to one of the groups listed above.
We test how changing the parameters np and s of the seizure detection algorithm affects the group assignment of the 26 onset patterns from patient 1 and all focal events, described above.  Fig 2) and the median times for varying both np and s. Specifically, for fixed s = 0.6 the median times computed for np 2 [40,80] give the same group assignment as the times for np = 60, s = 0.6 in 81% of the 26 events. For fixed np = 60 the median times computed for s 2 [0.4, 0.8] give the same group assignment as the times for np = 60, s = 0.6 in 96% of the events. The group assignment is robust to the choice of threshold modified by s in the detection algorithm, but is more sensitive to the envelope modified by np.

Mathematical model
Coexistence of stable (attracting) background and seizure states has been implied to underlie mechanisms of transition to epileptic seizures in the context of both generalized and focal seizures [17,18]. Bistability has been used effectively in phenomenological models of seizures driven by a desire to understand the underlying fundamental mechanisms of seizure transitions [8,11].
We consider a network of coupled bistable nodes where each node represents one electrode on either the EEG placed on the scalp or the silicone probe array used in the mEC. The network of nodes is given by a system of stochastic differential equations dz n ðtÞ ≔ ½ f ðz n ; n n Þ þ b X m A n;m ðz m À z n Þ� þ adW n ð1Þ for z 2 C with n = 1, . . ., N [7,9]. For the EEG data N = 19, the number of electrodes placed on the scalp, and for the mouse mEC data N = 16, the number of electrodes on the silicone probe.
The intrinsic bistable dynamics of each node is given by the following function, This is a truncated form of a Hopf normal form (or Bautin bifurcation) analyzed in [8,9] based on [10]. Note here that we substitute ν = 1 − λ in the equation from [8]. We choose the parameter 0 < ν < 1 so that f admits two stable states; the seizure state represented by a periodic orbit with phase ω and the background state represented by a fixed point. The bifurcation diagram for (2) is shown in S10 Fig. Each node starts in the background state (z = 0) and transitions to the seizure state are driven by noise when there is no coupling. Specifically, an independent identically distributed noise process W n representing input from other neural or external stimuli is added to each node. The noise amplitude α is fixed throughout at α = 0.05 as in [7,9]. Moreover, we assume that the time taken to transit back from seizure to background state is large enough to be ignored. When coupling is applied to the nodes, transitions to the seizure state are driven by both the noise and inputs from other nodes. S10 The parameter ν controls how easily each node can be destabilized by noise; hence ν represents the node excitability. When ν is close to 0 small noise perturbations is sufficient to push the node into the seizure state; the node is excitable. Whereas when ν is close to 1 much larger noise perturbations are necessary for the node to transition; the node is more robust, i.e. less excitable. Each node is assigned a value of ν, these values are chosen as described in the Excitability section below. For simulations in which we fix excitability equal for all nodes we use ν = 0.15 in line with [7,9] unless otherwise stated. Note that the transitions do not depend on the phase and so we fix ω = 0 for the simulations throughout unless otherwise stated. In this way we assume synchronization of nodes in the network and the model captures the observed phenomenon of phase locking at seizure onset.
The coupling between nodes is linear diffusive coupling, chosen in line with previous studies [7-9, 19, 20]. Such coupling could be seen as representative of local field potential mediated coupling via e. g. synaptic mechanisms at the level of the recordings obtained in the mouse mEC slice preparation and brain surface electrical potential coupling at the level of the EEG recordings. Nodes are coupled according to a weighted adjacency matrix A where entry A n,m 6 ¼ 0 if there is a connection from node m to node n and A n,n = 0. The overall intensity of the coupling is governed by the scaling factor β; if β = 0 the nodes are disconnected. We fix the scaling factor β = 1 when modeling the mouse mEC data. When modeling the EEG data we vary β according to functional connectivity networks reconstruction based on individual seizure events and specify the values in the text. In the case of EEG modeling A n,m is given by the association matrix computed for each seizure event as explained below.
Network structure. A plausible network structure between electrodes is required to apply the model to the EEG and mEC data. Linking mechanisms for neurons in a cohesive cortical network spanning the brain are not yet known [21].
For the EEG data, we use a standard measure of functional connectivity, namely a nonlinear regression method to compute the association index denoted A n;m 2 ð0; 1Þ that assigns a functional connection strength from node m to node n [22,23]. This association index represents a correlation ratio, i.e. the fraction of the variance of the time series x n from node n that can be "explained" by the input time series x m from node m. This method has been shown to be robust and applied to human intracerebral EEG data for characterizing seizure patterns [24][25][26]. Larger values of A n;m indicate a stronger functional connection between electrode activity and therefore between network nodes. We compute A n;m using the equation from [23]; detailed in S1 Text. The computations are implemented using a custom EEG-analysis toolbox [27] in MATLAB and the association indexes for the three example seizures are shown in S11 Fig. We use each association matrix as the weighted adjacency matrix A in the model.
To test the dependence of A n;m on the time interval used, we first compute A n;m for 5s intervals, 0-5s, 5-10s and 10-15s. Next, we compute the Pearson dissimilarity as defined in [28] between pairs of functional connectivity (or association) matrices. We apply this by comparing the matrices computed for 5s intervals, 0-5s, 5-10s and 10-15s with A n;m computed using the whole interval; see S12 Fig. The functional connectivity matrix characterized by the smallest Pearson dissimilarity is found to be the 5-10s. Note that this is the interval that contains the seizure. Finally, the onset patterns detected using the threshold detection algorithm show that onset ranges from 4-6 s from the start of the epoch. Therefore, to ensure all onset information is captured for each seizure we use A n;m computed based on the whole 15s epoch.
For the mouse mEC data, we choose a network topology based on the geometry of the recording positions of the silicone probe in a comb-like structure (see Fig 1) that we represent mathematically by a one-dimensional chain of nearest neighbor coupled nodes. Due to the fact that the slice is thin, only 400 μm, we ignore long range connections as they would be present (if any) with a low probability. Therefore, we assume a bidirectionally coupled chain where each node receives input from its neighbors A n,n−1 6 ¼ 0, A n,n+1 6 ¼ 0 and A n,m = 0, if m 6 ¼ n − 1, m 6 ¼ n + 1. We also assume that the coupling weight between two consecutive nodes is the same in each direction, so A n,n+1 = A n+1,n .
Excitability estimation from EEG data. We estimate the value for the excitability parameter for each node in the model from EEG data by means of a measure of energy from the signal from each electrode as in [29]. The time-dependent energy is computed using a 1s sliding window with 50% overlap as Here x n (kδ) is the time series of node n at time kδ with δ as sampling time step. The sum is taken over the time interval [t, t + 1]. We then compute the total energy as In particular for the two chosen seizures in the EEG data, we use t 1 = 0s and t 2 = 14s, consistent with whole interval used for the association index calculation. To ensure reasonable computing time for model simulations and to normalize the energy profile range for all events, we scale this energy profile E n to [0.1, 0.2]. The scaled energy profile is denoted by E n and these values are used to compute the node excitability as n n ¼ 0:3 À E n for each n which ranges between 0.1 and 0.2 unless otherwise stated. The energy profiles for the three examples are shown in S13 Fig. Numerical simulations and comparison with data. We simulate the system given by Eqs (1) and (2) using a stochastic Heun method with time step h = 10 −3 implemented in a custom code written in C++. We compute K � 1000 realizations for each set of parameters. The initial condition is z n = 0 for each node n = 1, . . ., N, where N is the total number of nodes in the network. Each realization has a different seed for the added noise processes. We detect the simulated seizure onset time τ n for each node n as the first time a realization crosses the threshold ξ as where the threshold is fixed throughout at ξ = 0.5. We then compute the recruitment time t n = τ n − min(τ n ) for each node. We average recruitment time over realizations T n ≔ ht n i. Note that the node with the lowest recruitment time min(t n ) is not necessarily the same in each simulation and so T s 1 ¼ ht s 1 i > 0, where σ 1 denotes the node with the lowest simulated mean recruitment time.
To compare the simulation results and the data we scale the simulated times by their maximum mean time t n /max(T n ) and denote the scaled mean recruitment time asT n ¼ T n =maxðT n Þ so the largest simulated scaled mean recruitment time maxðT n Þ ¼ 1. We scale the EEG and mEC data in the same way, by taking t n /max(t n ) andt n =maxðt n Þ respectively.
To quantify the model fit to the EEG data and compare between model simulations we compute the distance d between the simulations and the data. To this end, we split the 1000 simulation into 20 groups, order the nodes according to the mean recruitment time and average over all 50 realizations in each group. We then compute the least square distance between each group and the data. We calculate the median distance d and 95% confidence intervals using a standard method from [30]. We use the Mann Whitney U test [31] with the null hypothesis that d values are samples from continuous distributions with equal medians.
For the mEC data, we average over all realizations and calculate the least square distance d between the scaled mean times from the data and model. For the mouse mEC data we also compare the proportion of seizures where seizure-like onset starts in the ventral nodes 1-8 or in the dorsal nodes 9-16. Note, the proportion P V of activity starting in the ventral end has been already estimated in [13] for a larger number of samples (43) as P V � 0.86.

Seizure onset patterns in human EEG
The results of the linear discriminant analysis outlined in Material and Methods identify the following boundaries between the three groups, fast-domino, slow-domino and multi-domino onset. The slow-fast dividing line is given by L 1 = 2.9644 − 1.5236r − 16.5419l and the fastmulti dividing line is given by L 2 = 31.0766 + 2.301r − 97.5312l where r is the recruitment time and l is the maximum lag time. Each of the three groups, fast-domino, slow-domino and multi-domino onset, contain generalized seizures that are clinically classified as focal onset and generalized onset. This grouping goes beyond this clinical classification and so we wish to understand in more detail how these onset patterns arise. To this end, we apply our modeling framework to these three events.
Heterogeneous connectivity and excitability contribute differentially to different onset patterns. We build a network model for each seizure by constructing an all-to-all directed network of 19 nodes (one for each electrode) and explore the effect of changing the weighted adjacency matrix that defines the coupling, in our case the association matrix, and the intrinsic node excitability defined by the energy profile; see Materials and methods section for details.
We use the modeling framework to test the effect of different coupling and excitability in the model. Specifically, we perform three numerical experiments using: (Ex.1) homogeneous coupling strength given by the average of the association matrix elements A n,m with heterogeneous excitability given by the energy profile E n (shown in S13 Fig There is no significant difference between the d-values for Ex. 1 and Ex. 3 for either the slow domino onset (p = 0.78), or for the fast domino onset (p = 0.62). In this case the coupling weights are zero and the nodes escape independently governed by the excitability profile. This indicates that the excitability of nodes rather than the coupling is key to generating the fast and slow onset patterns.
To explore further the interactions observed in the in vivo data, we next apply our modeling framework to recordings of seizure-like neuronal activity from an in vitro experiment on mouse mEC.

Seizure-like onset in mouse mEC
In this section, we demonstrate the applicability of our modeling framework to an experimental system at a different spatial and temporal scale. We consider the recruitment of seizure-like events in six experimental recordings from slices of mouse medial entorhinal cortex (mEC) that has been previously reported in [13] and is summarized in Fig 1. Panels b and c in Fig 1  clearly show a domino effect of the recruitment of nodes to the seizure-like state that takes around 20s to spread across all nodes over a distance around 1500μm. We call this a slowdomino effect as it is much slower than the short interictal-like events propagated across the electrodes in the order of ms; see [13] Fig 3. We note that this domino recruitment is much slower than the total recruitment times for the seizures in the EEG recordings. This in vitro system allows for some simplifying assumptions with regard to the network connectivity and excitability to be made. This in turn enables a more detailed exploration of the relative contributions of excitability and coupling weights in seizure onset patterns. To this

PLOS COMPUTATIONAL BIOLOGY
Domino-like transient dynamics at seizure onset in epilepsy end, we apply our network modeling framework to the mEC data and assess the contribution of excitability and coupling to the domino-like onset patterns, specifically the recruitment times and the proportion of ventral initiation.
To model the onset of seizure-like events in mEC, we build a network consisting of 16 nodes where each node represents one electrode of the silicon probe array; see Materials and methods section. We start with a simple case where coupling weights and excitability are homogeneous and fixed among the nodes; in particular we fix ν n = 0.2, (n = 1, � � �, 16). In this case, we find that the simulated mean recruitment times exhibit symmetry; see S17 Fig. Therefore, in the homogeneous system there is no discernible dorsoventral gradient due to network symmetry and the proportion of simulations with ventral initiation is P V = 0.5. This indicates that homogeneous coupling weights and excitability cannot explain the sequential recruitment of seizure-like initiation along dorsoventral axis observed in the mouse mEC data. Thus, in the following subsections, we systematically explore how heterogeneous coupling or node excitability affect the domino-like recruitment along the chain.
Linear gradient in excitability facilitates slow-domino like sequential recruitment. Here we explore the effect of heterogeneous node excitability and constant coupling weight on recruitment times in the model. Experimental results have revealed the presence of functional dorsoventral gradients in a variety of excitability properties in neurons of the mEC [32]. In our modeling framework such experimental observations could be accounted for by changing the excitability parameter ν for each node linearly (in its simplest form). We note that this is the only parameter in the model that represents intrinsic properties of a node (in this case representing a small group of neurons in the mouse slice located in the vicinity of each of the silicon probes in the array). Accordingly, we define a linear gradient in ν given by with ν 0 , δν > 0 for n = 1, � � �, 16. Here we fix the coupling weights as follows, A n,n+1 = A n,n−1 = A = 0.1.  When there is no gradient (δν = 0), the simulated results are in line with the findings using homogeneous excitability, d � 2.5 and the proportion P V = 0.5. The solid lines delineating d = 0.2 show a band 0.001 < δν < 0.003 across a range of ν 0 where the difference between the model and the data is small. There is a region of overlap between this band and the region of (ν 0 , δν) contained within the dashed lines for P V . Fig 4b also shows that the gradient in excitability δν has a more significant effect compared to the initial level of excitability ν 0 in capturing the observed sequential recruitment. Note that the value of coupling chosen here is important for recruitment speed. When using a lower value of the coupling weight A = 0.03, there is no overlap in the regions of ν 0 and δν where d < 0.2 and within the confidence interval for P V (see S18 Fig).
Modulation of excitability gradient controls recruitment speed. Ridler et al. pharmacologically modulated the excitability levels in the slice [13] once the seizure-like activity was established with 4-aminopyridine (4-AP) by additional bath application of diazepam or Ro19-4603 (both GABA action modulators). It was reported that Diazepam decreased the speed of recruitment whereas Ro19-4603 increased the recruitment speed when compared to baseline (4-AP only). Using our modeling framework we are able to account for these experimental observations as shown in Fig 5a and 5b. Specifically in order to reproduce the chemically induced changes in recruitment speed it is sufficient to vary the excitability gradient δν in the model. For positive modulation we find that increase in δν suffice, resulting in a shallow excitability gradient and large difference in excitability between neighboring nodes; hence slowing down the domino effect. For negative modulation we find that decrease in δν is enough to produce a steep excitability gradient and small difference in excitability between neighboring nodes; hence speeding up the domino effect. Note that the electrode order in panel c is the same as in Fig 1a. As above we fix the coupling weight A = 0.1; we note qualitatively similar results were found for lower values of A. We take δν = 0.002 as the pre-modulation (4-AP only) value, as in Fig 4. When the gradient is much larger, δν = 0.025, the onset pattern has a very slow-domino effect and the total recruitment time is very large; here, due to the large recruitment times the realization is terminated before all nodes escape. This pattern of recruitment captures the very slow-domino onset pattern in panel a2 modulated by 4-AP + diazepam. When the gradient is decreased, δν = 0.001 the total recruitment time is small and the domino-effect is much faster, suggesting more synchronized recruitment. This captures the fast-domino onset shown in panel b2 modulated by 4-AP + Ro 19-4603. For δν = 0.001 the mean recruitment time for node 1 is larger than that for other values of δν and is far from 0 because when simulated node 1 is not always the first node to initiate the recruitment.
Linear gradient in coupling delineates regimes of slow and fast-domino recruitment. Next we proceed to explore the effect on recruitment times using heterogeneous coupling and constant excitability ν = 0.2 for all nodes in the model. It has been experimentally suggested that the dorsal mEC receives a greater number of inhibitory inputs than the ventral mEC [13,33] and is hence more resistant to transitions between different dynamic states. As above we look for the simplest possible way to account for these observations in our modeling framework. One way to do this is by assuming heterogeneous dorsoventral coupling weights. To this end, we incorporate a simple but general linear gradient in coupling weights along the chainlike network structure. Specifically, we use the weighted adjacency matrix where δA, A 0 > 0 for n = 1, . . ., 16. We note that in this case the coupling weights increase from the weakest between nodes 1 and 2 (ventral-fewer inhibitory connections) to the strongest between nodes 15 and 16 (dorsal-more inhibitory connections).  P V = 0.5 is in line with simulated results using homogeneous coupling. In panel b, the lines at d = 0.2 delineate the different coupling regimes identified in [9]. If A 0 small (A 0 < 0.02) then the coupling is weak and the network recruits as if uncoupled. For strong coupling A 0 is large, the recruitment of nodes becomes more synchronized (fast domino effect) and d is large (d > 0.2). The region with lowest difference to the data (d < 0.2) is within the intermediate coupling regime, between the solid lines, in which the slow domino effect can be seen. Panel c shows that for δA 6 ¼ 0 the proportion P V � 1. We find values (e.g. indicated by the red star) of A 0 and δA where d < 0.2 and P V � 0.86 in line with [13].
This shows that the network model with a linear gradient in coupling weights delineates regions of slow and fast domino-like recruitment. Moreover, this model captures the behavior of recruitment times seen in the experimental data which is well explained by the slow-domino effect.
Severed coupling facilitates time lag in recruitment. Finally we test how modulation in coupling affects recruitment patterns. To this end, we consider the anatomical separation of the dorsal and ventral portions of the mEC slice via an incision in the intermediate mEC with a scalpel blade, as reported in [13]. Seizure-like activity was found to initiate in the ventral mEC before the dorsal mEC, but transitions in the dorsal mEC took longer to initiate compared to intact slices. This suggests that the dorsal mEC is less likely to produce seizure-like activity in the absence of the ventral mEC. We model this scenario and predict the mean recruitment times along the 16 nodes in the chain. Specifically, we model the anatomical separation by setting A 8,9 = A 9,8 = 0. Moreover, we fix the coupling weights  When the chain is severed, recruitment is disrupted. The recruitment times form two clusters (the ventral cluster with nodes 1-8 and the dorsal cluster with nodes 9-16) separated by a lag in recruitment time. When the coupling is weak A = 0.03 the lag is larger than the intact but could still be considered as a slow-domino effect. However, the size of the lag is more pronounced when the coupling weight is intermediate A = 0.1 and the onset pattern is clearly multi-domino with two cascades. This result is in line with experimental findings where the dorsal mEC takes a long time to transit to seizure-like state in the absence of influence from Fig 7. Separation of dorsal and ventral nodes produces a multi-domino effect. Panels a and b show simulated mean recruitment times for both intact (green) and separated (blue) chain structures with indicated coupling weights A. In both cases the separated chain structure shows a lag between nodes 8 and 9. Panel a shows that when the coupling strength is weak the separated simulation has a slow-domino or multi-domino onset with one larger lag than the intact onset, whereas in panel (b), A = 0.1, the lag is larger and two domino cascades clearly form.

Discussion
The variety of generalized seizure onset patterns observed in clinical EEG recordings have not been fully explored. EEG is the dominant tool used by clinicians and neurologists to diagnose epilepsy therefore analysis and modeling of the complex patterns observed may provide crucial insight for diagnosis and treatment of this complex disorder [34].
Using an exemplary data set of recordings from individuals with idiopathic generalized epilepsy (IGE) [4,12] we show a variety of spatiotemporal onset patterns for seizure in the data that can be characterized by the total onset (recruitment) time and the lags in onset time between electrodes. Here we use these two features to suggest a novel taxonomy consisting of three types of patterns: fast-domino onset where activity is detected almost simultaneously by all electrodes; slow-domino onset where the onset of activity takes longer to be detected in all electrodes; and multi-domino onset where several domino-like cascades of activity are observed across groups of electrodes. These novel categories transcend the existing classification of generalized or focal onset of a generalized seizure.
We develop and apply a network modeling framework of seizure onset to capture different seizure onset patterns characterized by the total recruitment time and lags between electrodes. We demonstrate, as a proof of concept, that this novel framework can be applied at different scales and used to explore the contribution of node level properties on emergent network behavior.
To this end, we apply this network model to qualitatively capture three seizure onset patterns, one from each onset pattern category. The illustrative examples used here are taken from one individual diagnosed with Juvenile Absence Epilepsy (JAE, the most common syndrome within the data) [4,12]. The precision approach taken here by modeling individual seizures is in line with contemporary approaches to epilepsy modeling based on inference of the network structure and node excitability (propensity to transition) from the data for each seizure [19,20]. Our modeling results show that heterogeneous node properties and heterogeneous coupling weights play non-trivial roles in the formation of the different onset patterns. We show that, of the two heterogeneities, heterogeneous excitability plays a more dominant role in shaping the emergent behavior for all the three onset patterns. Indeed, heterogeneous excitability alone is sufficient to reproduce the fast and slow-domino onset patterns. The multi-domino onset pattern is best captured by incorporating both heterogeneous coupling and excitability, where the heterogeneous coupling provides a refinement to the multi-domino onset pattern.
In order to further investigate the contribution from node excitability heterogeneity or coupling heterogeneity separately, and to demonstrate that the framework can be used at the micro scale, we apply the network model to multi-electrode array recordings from mouse mEC [13]. The mEC data shows an induced domino-like seizure recruitment from the ventral to dorsal mEC. To model this, a chain network structure is assumed based on experimental set up. The effects of gradients in node excitability or coupling weights are explored based on model simulations. We capture experimental findings by showing that steeper gradient in excitability leads to fast-domino onset patterns where nodes are recruited in quick succession. Moreover, we show that shallower gradient in excitability leads to slow-domino onset patterns where for very steep gradients, long recruitment times are found and can cause the propagation to effectively stop. We also show that a discontinuity in the coupling leads to a time lag in recruitment in the chain, creating a multi-domino onset effect; these predicted recruitment patterns could be experimentally tested in the future.
The idea of groups of electrodes which detect similar changes in activity at seizure onset or termination was explored by Proix et al. for focal seizures [35]. The authors considered the spatiotemporal patterns of seizure recruitment and termination detected by stereotactic EEG. They identified that although some seizures begin and terminate synchronously across electrodes, others form distinct clusters where nodes within each cluster are recruited or terminate (almost) simultaneously, with a significant time lag between each cluster. They focused on termination patterns and use a neural field model based on an extension to the so-called "Epileptor" model [36] to show that recruitment correlates to weak (structural) connections. This paper extends the concept to the onset of generalized seizures and our finding that coupling plays an important role in the multi-domino onset pattern is in line with their findings.
Network topology alone is a necessary but not sufficient condition for generalized seizures to emerge. For example, Chowdhury et al. demonstrated that alterations in brain network topology were present in both people with idiopathic generalized epilepsy, as well as the unaffected first-degree relatives in comparison to healthy controls [37]. Moreover, Petkov et al. demonstrate increased mean degree is related to seizure onset in generalized seizures [38]. Consequently, Terry et al. studied a phenomenological model of seizure onset in which each node could transition between background and seizure states [11]. They found that increasing the excitability of a single node in a motif (four-node) network could produce dynamics consistent with either focal, generalized or focal to bilateral ictal activity. In that model, the long term behavior of the system (how many nodes are in the ictal state) was considered. Building on these initial mathematical descriptions, Schmidt et al. extended this framework to consider a modular network of tightly coupled Kuramoto type oscillators representing activity within each node, and a functional network between nodes informed by clinical EEG [39]. Therein, onset patterns either defined by cycles within the functional network or synchronous behavior within specific nodes were described. In the context of generalized epilepsy, this approach was shown by Schmidt et al. to have predictive value for revealing epilepsy from epochs of EEG data that would be considered clinically negative (e.g. free from discharges or other abnormalities) [40]. Very recently, this approach has been further extended by Woldman et al. to account for different network properties characterizing both focal and generalized seizures [41]. Martinet et al. show that network characteristics can influence recruitment dynamics of the neocortex in focal to bilateral tonic-clonic seizures [42]. In particular, recruitment times and the degree of spatial organization during onset have been shown to be patient specific and large variations exist between patients. More recently a cascading failure model was used to simulate the neural networks underlying generalized tonic-clonic seizure [43]. They use functional networks derived from graph theory and initiated seizures via stimulation of the node with the largest number of connections. They do not explore the effect of inherent excitability of the nodes to seizure initiation.
In this framework we use a standard measure of functional connectivity from the EEG to inform our coupling structure and weights. Due to the ambulatory nature of the EEG recordings we chose an all to all connected network structure. For each of the two seizure events studied in detail here we computed and compared the functional connectivity matrices obtained from different sections of the event epoch as well as by using the entire epoch. The differences between matrices computed from within one epoch were found to be smaller than the difference between the two epochs. Therefore we chose to use the matrix computed from the entire epoch in each case. There is evidence that functional connectivity changes with time, particularly before and during epileptic seizure [44][45][46]. A future development of the model could be to include a dynamic coupling structure to explore the role of evolving connectivity in conjunction with excitability on seizure onset patterns including multiple seizure events and on a longer time scale. Furthermore, cortical excitability is known to change before, during and after a seizure [1]. A further development of the model could be to include a dynamic excitability on the nodes. This could be incorporated by making ν a slow variable on the order of minutes or hours compared to the total recruitment time of the seizure onset on the order of milliseconds.
A number of computational and mathematical models have explored the influence of intrinsic node properties on the dynamics of recruitment to a pathological state; see for example [47][48][49]. Notably, Goodfellow et al. showed that spatial heterogeneity is involved in the onset of absence seizures [50]. Building on the computation models of spike-wave discharges observed in EEG during absence seizures by [51] and further extended by [52,53], Goodfellow et al. develop a spatial extended Jansen-Rit model and show that seizure onset may relate to the influence of local heterogeneous excitability. However, they do not explore how the interplay between connectivity and heterogeneous excitability produce different seizure onset patterns.
There are multitude approaches to elucidating brain region excitability measures from neuroimaging data. Most are applied in the problem of epilepsy surgery for Focal epilepsies (rather than generalized epilepsies considered here). The first and second generation Epileptor models have been used to model and predict seizure recruitment of local and distant brain regions in focal epilepsy [20,36,54,55]. The authors characterized regimes of recruitment behavior for different coupling weights [55] but do not consider excitability gradients or recruitment timings. Jirsa et al. take two approaches to identify a brain excitability map for their Virtual Epileptic Patient described in [20]. Having constructed a personalized Epileptor network model using structural connectivity measures from diffusor tensor imaging (DTI), they apply heterogeneous excitability to the nodes first as prescribed by an expert clinician then using excitability elucidated form functional stereotactic EEG data using Bayesian inference methods. Using the first they are able to simulate seizure patterns observed from the patient. With the second they use the excitability measure itself to identify the seizure onset zone to inform surgery for resection. Another approach to approximate excitability was taken by Hutchings et al [56]. The authors elucidate heterogeneous measures of excitability and connectivity using DTI data from individual patients with temporal lobe epilepsy. They then apply a Bautin bifurcation model based on [8] to simulate surgery (resection of nodes) and predict outcome success.
Although our model contains oscillatory (seizure) states it does not take into consideration, for example, the frequency of oscillations as all nodes are synchronized by design. The results presented here do not depend on the choice of node model given by (2). We note that we repeated the model simulations in the mouse mEC section using a simple bistable model with two steady states (rather than one oscillatory) on each node and the results are consistent. Several other mathematical models have been designed to reproduce the activity and transitions between resting and seizure states, including neural mass models [19,48,54]. These could be used to incorporate different levels of biophysical reality. Moreover, a generalized Hopf model that incorporates additional phase terms with phase dependent coupling could be used to explore the effect of synchronization on onset patterns. We consider these beyond the scope of this paper.
A future application of this framework would be to model the non-seizure events that have been observed in resting state EEG recordings and are also generated in mouse mEC experiments [13]. Expanding the framework to investigate the initiation, recruitment and recurrence of inter-ictal spikes may provide further insight into the structural and functional alterations indicative of different types of seizure and inter-seizure events. This expansion would rely on using suitable intrinsic node dynamics in place of Eq (2). The question of which node dynamics to impose relies on making assumptions regarding the mechanism that defines the transition from background activity to a seizure or inter-ictal spikes; see for example the comparison of the choice of node dynamic models on outcome in the context of epilepsy surgery [49].
To assign seizure onset patterns to the fast, slow, and multi-domino onset groups we use linear discriminant analysis. We show that this assignment is robust to choices of the parameters of the seizure onset detection algorithm in S6-S9 Figs and note that the median over a range of parameters could be used as an alternative to the fixed choice of parameters used here. The lines delineating the groups are based on an initial assignment of points by visual inspection (see Materials and methods). The onset patterns could instead be submitted to a machine learning algorithm to formulate a classification. In this study, for consistency with the other measures, we use the entire EEG epoch of the data provided to extract seizure onset pattern. A smaller interval around the seizure onset as identified by the clinician could be used to give a clearer grouping; for example S19 Fig shows the results of submitting the first 8 seconds of each EEG epoch to the seizure detection and linear discriminant analysis. Whilst for the chosen subset of events the groups are somewhat more clearly distinguished, S19 Fig also shows the spread of the generalized onset patterns over all three groups is maintained. The focus of this paper is to highlight the various clinically observed onset patterns and apply a modeling framework to identify key components to how they are generated. An important step towards implementation of this classification in clinical practice would require verification by future studies and identifying connections between the groups and, for example, treatment outcome.
S9 Fig shows that the multi-domino class is actually the most robust, as the median lies in the same classification as the point for np = 60, s = 0.6 for all multi-domino onset events. This may be due to the fact that the multi-domino onset region is the largest of the three in the recruitment/lag plane. We show that the exact order of electrodes (based on detected activity) varies with the parameters of our detection algorithm, but our novel assignment of dominoonset pattern is robust to these perturbations; see S6-S8 Figs.
We are using 19 electrode EEG, which is unsuitable for source localization techniques and so we cannot make conclusions about the excitability of underlying brain regions nor place significance on the precise order of electrodes here. Future work would be to incorporate high density or intracranial EEG recordings to link the electrode order in the domino cascades to the underlying generators of epileptic activity. In particular, we observe for the multi-domino onset example shown in S8 Fig that for the onset patterns classified as multi-domino onset the maximum lag always precedes recruitment of electrode C4. We note that the association matrix shown in S11 Fig, electrode C4 has lower connectivity values than other nodes for the multi-domino event. This is incorporated in our modeling results shown in Fig 3 where row a column 4 shows the best fit for the multi-domino case where the largest simulated lag proceeds electrode C4 (not labelled). This raises the question of how the persistence of this lag relates to the underlying neurodynamics.
We present here, a variety of clinically observed onset patterns and a modeling framework applied at multiple scales to show that the interplay between excitability and coupling is required to generate them. Together these constitute a step toward a more comprehensive understanding of different types of seizures and therefore a more advanced classification of the epilepsies. patient with all the focal onset events for np = 60 and s = 0.6. The fast, slow and multi-domino example are in darker blue, purple and orange, respectively, with L1 and L2 marked; this is the same as Fig 1a but here is plotted on a log scale. The other two panels show these 26 points each with a corresponding median value (+) that indicate the results of varying np for s = 0.6, and varying s for np = 60. For fixed s = 0.6 the median values computed for 40 < np < 80 give the same classification as the points np = 60, s = 0.6 in 81% of the 26 events. For fixed np = 60 the median values computed for 0.4 < s < 0.8 give the same classification in 96% of the events. The algorithm is very robust to the choice of threshold modified by s but is more sensitive to the spline interpolation modified by np. We note that the multi-domino group is robust to all values of both np and s. (EPS) S10 Fig. Panel a shows the bifurcation diagram of Eq (2) for fixed ω = 20. For ν < 0 there is one unstable equilibrium and a stable limit cycle. At ν = 0 there is a subcritical Hopf bifurcation that gives rise to an unstable limit cycle and stabilizes the equilibrium. For 0 < ν < 1 the system is bistable. At ν = 1 there is a saddle node of limit cycles and for ν > 1 there is one stable equilibrium. We restrict our ν values to the bistable regime where the stable equilibrium represents the background state and the stable limit cycle represents the seizure state. Panel b shows an example of node trajectories from the network model given by (1)