Epileptor network model

The epileptor network model has been described in several previous studies (e.g. [6, 11]). For completeness we include here the equations following closely the notation in [10]. For an N-node patient-specific epileptor network model, the dynamics are given by (1) (2) (3) (4) (5) where (6) and for i = 1, 2, …, N. The terms ξ_i(t) and η_i(t) correspond to stochastic inputs. (We effectively interpret the above as stochastic differential equations in the Itô calculus sense; see Model simulations). We set the parameters I₁ = 3.1, I₂ = 0.45, γ = 0.01, τ₀ = 6667, τ₁ = 1, τ₂ = 10.

The coupling weights W_ij are obtained from patient-specific connectivity matrices derived from white-matter tractography with corresponding axonal transmission delays τ_ij (see Structural network connectivity).

The parameter x_0,i denotes the neural excitability in the i^th-node. An epileptogenic zone (EZ) can be instantiated by setting x_0,i = −1.6 in the corresponding node. This excitability level ensures that an isolated EZ node undergoes a seizure. We set x_0,i to non-epileptic values in the non-EZ nodes, specifically x₀ ∈ [−2.3, −2.09] as the critical excitability for an isolated node is , i.e. a node i is considered epileptogenic (an EZ node) if . The parameter w corresponds to the global connectivity strength.

Briefly, given the time scale separation τ₀ ≫ τ₂ ≫ τ₁, slow and fast oscillations in typical epileptic focal seizures are captured by the systems (x_1,i, y_1,i) and (x_2,i, y_2,i), respectively, while the slower time-scale permittivity variable z_i might reflect the generic dynamics of changes in extracellular ionic concentrations (e.g. potassium), metabolic factors, etc. The diffusive coupling [x_1,j(t − τ_ij) − x_1,i] incorporates the hypothesis that a seizure spreads by epileptic activity in an EZ, perturbing the ionic/metabolic homeostasis in a distal area. This perturbation then results in distal areas also going into seizure. For a detailed analysis of the epileptor dynamics see [7, 8, 21].

Structural network connectivity, time delays and epileptogenic areas

Patient-specific connectivity networks were inferred from white-matter tractography obtained via diffusion MRI [6, 22]. In patient-specific networks P1—P2, a 84-area (Desikan-Killiany) parcellation of brain areas was used, while for P3-P5 a finer 162-area parcellation was employed.

The inferred patient-specific connectivity matrices W were further truncated and normalized. To attenuate potentially large artifacts in the estimation of white-matter connectivity, entries larger than the 95% percentile (across all entries) were set to this threshold value. After that, all entries were normalized by this threshold such that 0 ≤ W_ij ≤ 1. Connectivity matrices for patients P1—P5 are shown in S1 Fig in S1 File.

Time delays τ_ij in the model were computed from the length of the estimated white-matter fiber tracts connecting brain areas related to nodes i and j [22]. The maximum length in the 5 different patient-specific connectivity data corresponded to 200 mm. We set the speed to 60 mm per time unit in the model. Given our choice of setting the time unit in the model to 0.02 s (see Model simulations), a speed of 3000 mm / s (i.e. 60 mm × 1/0.02s) was obtained, with a maximum time delay of 0.067 s for the largest fiber length of 200 mm. This is slower than expected when compared to only axonal transmission delays, and it should be understood here as including also delays in synaptic activity and rising times in neuronal population responses, for example. The time delays for each patient-specific network are shown in S1 Fig in S1 File. We also explored a range of much slower delays, but found no qualitative differences in the spreading dynamics in these epileptor network models. We elaborate on related issues in the Discussion section.

The following nodes and brain areas were identified as the epileptogenic zones (EZs)

1. P1: {61, 64} for ctx-rh-lingual and ctx-rh-parahippocampal, respectively;
2. P2: {48, 60, 81} for right-amygdala, ctx-rh-lateral orbito frontal and ctx-rh-temporal pole, respectively;
3. P3: {135, 154, 156} for right-rhinal-cortex, right-hippocampus-anterior and right-amygdala, respectively;
4. P4: {35, 51, 53, 54, 73, 74, 75} for left-temporal-pole, left-collateral-sulcus, left-parahippocampal-cortex, left-rhinal-cortex, left-hippocampus-anterior, left-hippocampus-posterior, left-amygdala, respectively;
5. P5: {83, 87, 116, 126, 135, 154, 156} for right-orbito-frontal-cortex, right-F3-pars-opercularis, right-temporal-pole, right-T2-anterior, right-rhinal-cortex, right-hippocampus-anterior, right-amygdala, respectively.

Model simulations

The epileptor network model (Eqs 1–5) were implemented as a system of stochastic differential equations in the Itô calculus sense, where ξ_i(t) and η_i(t) thus denote Wiener processes. For numerical integration we used the Heun’s method, with the Wiener process increments implemented as i.i.d. Gaussian random variables with zero mean and variance σ² ⋅ Δt. The above five-dimensional node system was implemented as a 6-dimensional system by rewriting Eq 6 as corresponding differential equation. We wrote Python code to implement the numerical simulations and used an integration step of 0.05. To match typical oscillatory properties of focal epileptical seizures and for plotting purposes, we assumed a time unit in the model to correspond to 0.02 s. Given the simulation step size of 0.05, this resulted in a sampling rate of 1000 simulation steps per second.

The epileptor network was initialized with an initial conditions near identified stable equilibrium point of the single node epileptor model. The excitability parameter was set to the same nonepileptic value for all nodes (including all EZ nodes). Simulations were run for 20,000 time steps to let transients decay and the network to settle into a steady-state. At that point, the excitability for the selected EZ node was set to x₀ = −1.6. Then, a seizure may or may have not developed in this active EZ node after some varying time interval depending on initial conditions and parameters. After a seizure terminated, we enforced a postictal (refractory) period by uncoupling the corresponding node where the seizure terminated from the network. This postictal period was enforced to all nodes regardless if the node was an EZ node or not.

We take seizure spread as a collective variable (see below) and examine its behaviour by taking x_0,i in the non-EZ nodes and w as varying control parameters. In addition, we examined noise drive with different variances and the effect of the location of the selected active EZ node. We varied x_0,i in the non-EZ nodes in a non-epileptic range, as well as the global coupling strength and the noise variance (e.g. Fig 2). We applied a grid to the parameter space (x_0,i, w). For each point in the grid, we generated 30 stochastic realizations. The same approach was then applied to the case where a different node was selected as the active EZ node.

Epileptor network observations, seizure spread as an order parameter, and critical dynamics

For plotting and visualization purposes (e.g. Fig 1), each node in the full 6-dimensional epileptor node is observed via the difference of the variables x_2,i and x_1,i, which are related to the slow and fast ictal oscillations, respectively.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Seizure spread in the epileptor-network model.

Three qualitatively different behaviors are observed for different global connectivity strength (w) and surrounding-node excitability (x₀) parameter values. Figures show simulations of the epileptor network model for patient P1. Node 64 was identified as the target active EZ node (x₀ = −1.6). The standard deviation of the noise input level was set to σ = 0.05. (A) The active EZ node (in red) was “inhibited” by the surrounding nodes in the network and, as a consequence, did not transition into a seizure. Global connectivity strength and excitability in the surrounding nodes were set to w = 1.6 and x₀ = −2.3, respectively. (B) The seizure started in the target EZ node, but there was no spread in the network (w = 0.9, x₀ = −2.2). (C) The seizure started in the EZ node near time 50 s and spread to all of the network nodes (w = 0.9, x₀ = −2.1). Spread onset times varied across different nodes. Partial spread to a subset of nodes can also occur (see Figs 2 and 8). D A zoomed in view of the activity shown in panel C. After a seizure terminated in any given node, the node went into a postictal refractory state implemented by uncoupling the node from the network and, in the case of the EZ node, the excitability parameter x₀ was switched back to the same non-epileptic level of the surround. Interictal spikes are observed prior to the seizure onset. The decrease in the spike rate after the seizure terminates is due to both a “refractory” effect of the intrinsic node dynamics and to the enforced uncoupling.

https://doi.org/10.1371/journal.pone.0272902.g001

In order to quantify seizure spread and define a corresponding collective state variable, we ignored the details of the temporal dynamics in each node of the simulated epileptor network, and considered only discrete states, i.e. whether a seizure had spread or not to a given node, together with the temporal ordering of the spread, i.e. which nodes spread first and so on. As seizure activity showed large amplitude deviations from normal activity, specific thresholds were easily identified to ensure accurate detection of seizure onset and termination. A seizure was detected to spread to a given node when the activity in the corresponding state variable g (Eq 6 rewritten as a differential equation) crossed a specified threshold. We chose the variable g, instead of x_1,i or z_i, because it presented a smoother activity.

Thus, the state of each node was then represented as a binary variable (1 for seizure and 0 otherwise). The collective state variable or order parameter was then defined as the spread size (number of nodes that entered seizure activity) across the network in any given simulation.

Given simulations of epileptor networks and the observed behavior of the defined order parameter, we built empirical phase diagrams as a function of control parameters, here excitability in the surround nodes and global connectivity strength. Qualitative changes in the order parameter, e.g. spread vs no spread, were used to construct phase diagrams and phase transition curves, i.e. curves tat separate distinct phases.

We refer to the region near these phase transition curves as near-criticality regions. As stated in the Introduction, here we use the terms critical and criticality in a broader sense than in the strict use in statistical mechanics. Briefly, in that context, critical point and criticality refer strictly to second-order phase transitions [23, 24]. In second-order equilibrium phase transitions, the order parameter (in many cases given by the 1st-order derivative of the free energy with respect to an external field) changes continuously as control parameters are varied around a critical point (i.e. there are no jumps in the phase transition), while response functions (second-order derivatives of the free energy wrt control parameters; e.g. magnetic susceptibility in Ising models) and correlation length diverge as power-law functions at the critical point in the thermodynamic limit N → ∞. In nonequilibrium systems [24], critical phenomena has mostly been characterized as absorbing state phase transitions [25], self-organized criticality [26], and branching processes [27] where the distribution of event (avalanche) sizes and durations in the system appear as power-law scale-invariant functions. Here, in the case of seizure spread in epileptor network models, the exact nature of the nonequilibrium phase transitions in the defined order parameter (spread size) remains an open question. For this reason, we use the terms critical and criticality in a broader sense where they relate primarily to bifurcations in stochastic dynamical systems. In this case, a critical point refers to a specific value of a control parameter where a bifurcation happens, regardless of whether the order parameter changes continuously or jumps, and whether response functions diverge [19, 20].

Local linear stability analysis

Fixed points for the epileptor network model were computed by solving for the zeros of the deterministic system of differential equations (Eqs 1–6) with the Matlab function fsolve (Mathworks, Inc). We first determined the fixed points for a single-node epileptor. Given our choice of parameters and initial conditions for the numerical search, we found three fixed points: one typically stable, and the other two unstable. When solving for fixed points of the epileptor network, we started with initial conditions near to the stable fixed point found for the single-node epileptor model.

For local linear stability analyses, we write the Jacobian matrix for the 6-dimensional N-node epileptor network model as: (7) where each component is itself an N × N matrix. The non-zero components are , , , , , , , , , , , , , , J_g,g = −γ I, and is given by (8) In the above, I is the N × N identity matrix, is a diagonal matrix with the i^th diagonal element equal to ), and δ_i,j is the Kronecker delta function.

As explained in the Results section, we consider two instantiations of the above Jacobian matrix. In the first case, the EZ node is included in the analysis and the Jacobian remains a 6N × 6N matrix. In this case, we set x_0,i = −1.6 for the EZ node, and a non-epileptic x₀ value homogeneously for all of the remaining nodes, and proceed with the linear stability analysis.

In the second case, the EZ node is replaced by an “external” constant input to its target surrounding nodes and the Jacobian is a 6(N − 1) × 6(N − 1) matrix obtained by removing the rows and columns related to the EZ node. The active EZ node contribution to the nodes in this network was replaced by a constant input computed by averaging x_EZ,1 during the seizure period.

This constant input to the nonepileptic nodes was mediated via the same diffusive couplings via the permittivity variable z as in the usual epileptor network. For example, assuming the i^th-node as the active EZ node, its contributions to the variable z_j in the jth nonepileptic node was given by . We note that this constant input does not appear in the Jacobian matrix itself, but it changes the fixed point at which the Jacobian is evaluated (The fixed points of the 6(N − 1) × 6(N − 1) system). The constant input shifts the fixed points of the system towards instability.

The combined use of the eigenvalues of both Jacobians lead to the best prediction of the phase diagrams for seizure spread in the epileptor networks. We also note that time delays were ignored, i.e. we set τ_ij = 0 in the derivation of the Jacobian matrices. Proix et al. (2014) have argued that because the interaction between nodes is via the slow permittivity variable z, axonal propagation delays play no major role on bifurcations of equilibrium points and thus on the seizure spread dynamics. Each Jacobian matrix was evaluated at the fixed point of the corresponding epileptor network system.

We also considered local stability analysis and related predictions based on a 2-dimensional reduction of the epileptor model given by [6, 9] (9) (10) where τ₀ = 6667 and I₁ = 3.1.

Prediction of the temporal order of node-recruitment during seizure spread

When a seizure spread is predicted, i.e. the real part of the maximal eigenvalue of the Jacobian matrix satisfies λ_R > 0, we used the corresponding eigenvector to predict the order of spread or node-recruitment during the seizure. Specifically, we computed the magnitude of the leading eigenvector components corresponding to the variable z_i, for nodes i = 1, 2, …, N in the surrounding, i.e. excluding the active EZ node. Nodes with larger component magnitudes were predicted to be recruited earlier during the seizure than nodes with smaller component magnitudes.

Epileptor network model: Seizure simulations

Focal epileptic seizures were simulated with the epileptor network model. The intrinsic dynamics in each network node is formulated as a 6-dimensional ODE system, which captures the generation of interictal spikes, slow and fast oscillations during the seizure itself, and a much slower process (permittivity variable) that spontaneously brings the network in and out of seizures (Fig 1; Methods).

The networks were instantiated with 5 different patient-specific (inter-node) connectivity matrices, axonal propagation time-delays, and corresponding identified epileptogenic zones (Methods). In two of the cases (P1 and P2), connectivity matrices were obtained with an 84-node parcellation of the brain areas, and 162-node parcellation in the remaining (P3-P5). Corresponding connectivity weight matrices together with descriptive graph theoretic measures are given in the S1 Fig and S1 Table in S1 File. A list of the EZ nodes and corresponding brain area labels is given in the Methods section. While the connectivity matrices are symmetric and non-negative, diffusive coupling via the permittivity variables in different nodes allows not only for excitatory but also inhibitory and restraint effects in the network (Fig 1; Methods).

Each seizure was simulated by starting the entire network in a non-epileptic excitability level (x₀ parameter in the model; Methods). After a specified time to allow network dynamics transients to settle down, the excitability level for a single selected target EZ node was set to an epileptic level (x₀ = −1.6). This target EZ node may or may have not then transitioned into seizure after varying times, and led to partial or full seizure spread, depending on the level of surrounding node excitability (homogeneous x₀ across nodes), global connectivity strength (parameter w), and on the different stochastic realizations for a given input noise level (Fig 1). We simulated thousands of stochastic realizations under variations of the above parameters and EZ node selections.

Seizure spreading dynamics behave akin to phase transitions

As stated above, we defined seizure spread as a collective variable given by the spread size (or equivalently fraction of seizing nodes in the network) for any given patient-specific epileptor network simulation. To characterize qualitative different states or phases of this collective variable, we computed phase diagrams as a function of control parameters given by the surround excitability, global connectivity strength, noise levels and location of the active or target EZ node. Phase diagram examples for three different patient-specific networks, active EZ node, and parameter settings are shown in Fig 2. (See S1 File for more examples from all 5 patients and the respective active EZ nodes).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Seizure-spread phase diagrams.

Seizure spread is examined as a collective state variable, i.e. spread size or number of seizing nodes, under variations of global connectivity strength (w) and homogeneous surrounding node excitability (x₀). Top (A,B,C) The phase diagrams were obtained by simulating the epileptor network (patient P3, node 135 corresponds to the active EZ, with x₀ = −1.6) for three different noise strengths (σ = 0.05, σ = 0.10, σ = 0.15). Blue, green and yellow colors represent no seizure, no spread and seizure spread, respectively, obtained from 30 realizations of the stochastic dynamics for a given parameter setting (each point in the diagram). The colorbar indicates spread size. (D,E,F) Fluctuations across stochastic realizations are observed in the near-criticality region between the no-spread and spread phases as indicated by the standard deviation of the spread size (30 different realizations). The gray background corresponds to regions of no observed fluctuations. (G) Phase diagram predicted by the local stability analysis applied to the deterministic 2D-reduced epileptor network model (Methods). (H) Phase diagram predicted by the analysis of the deterministic original epileptor network model. Prediction of the different phases was based on the leading eigenvalues of two different Jacobian matrices, one that included the active EZ node and another that replaced the EZ node with a constant input (see Methods and Fig 3 for a schematic description). Although the true location and shape of the transition curves differ, the predictions capture well the qualitative aspects of the phase diagram. Middle (A-H) The panels show another example from patient P1 with node 64 as the active EZ. Bottom (A-H) Example from patient P2 with node 60 as the active EZ, showing distinct fluctuation regions in contrast to the above examples.

https://doi.org/10.1371/journal.pone.0272902.g002

Three main phases with clear phase transition curves appear: (a) no spread, i.e. the seizure remains localized in the target EZ node and does not spread to the surrounding network nodes; (b) the seizure spread across the network with various different spread sizes; in particular, partial spread to a subset of network nodes is observed near a small near-criticality region surrounding the phase transition curve; (c) no seizure, i.e. a seizure does not occur even in the active EZ area, in other words, the surround network activity appears to inhibit seizure onset in the target EZ area.

The existence of these three phases and their qualitative features were general across all patient-specific networks, and across different active EZ areas in the same patient (S2-S8 Figs in S1 File). These findings emphasize the role of excitability in the surrounding network and of global connectivity strength, as well as the role of inhibitory and restraint effects enabled by diffusive coupling, on seizure spread. Our results extend the 2-node epileptor network analysis in [6] to patient-specific networks with much larger number of nodes and under broader parameter and noise level variation.

Phase diagrams of seizure spread show large stochastic fluctuations in a narrow near-criticality region

By running simulations with different random noise realizations and noise levels, but fixed excitability and global connectivity strength and target EZ node, we examined the existence of stochastic fluctuations in the defined collective variable for seizure spread. Additive noise in the epileptor network model reflects the combined effect of various types of neural noise, including thermal noise in ion channels, failure in neurotransmitter vesicle release in synapses, and many other unaccounted factors.

The existence of stochastic fluctuations in the collective variable for the same levels of surround excitability and global connectivity strength has important practical implications, for instance in the predictability of seizure spread. Large stochastic fluctuations of the collective variable across different seizures under the same parameter setting make predictability harder even when one has a good estimate of surrounding neural excitability levels and global connectivity strength.

We observed stochastic fluctuations only near the phase transition curve separating the no spread and spread phases (Fig 2D–2F). Fluctuations were confined to the small near-criticality region surrounding this curve. Away from this region, seizure spread dynamics appeared deterministic for the examined noise levels and number of stochastic realizations. Similar results were obtained for other patient-specific networks (S2-S8 Figs in S1 File). In a small number of cases, unconnected fluctuation regions appeared. They consisted of two side-by-side regions with a no-fluctuation zone in between (e.g. Fig 2D–2F, bottom panel, P2, node 60 as the active EZ; S4 Fig in S1 File: P3, node 154 as the active EZ). This phenomenon can be understood as follows. In those cases, a couple of nodes were strongly connected to the active EZ node. As the the parameter variation moves toward traversing the phase transition boundary, say by increasing both {w, x₀}, initially only those strongly connected surrounding nodes go into seizure and with fluctuations across realizations. Then, with a further slight increase, these nodes always go into seizure, so there are no fluctuations. With a further increase, the seizure starts to spread more broadly and with fluctuations.

The nature of the fluctuations in this near-criticality region may have important implications for prediction and control of seizure spread. In particular, if focal seizure dynamics in actual brains follow similar dependencies on excitability and global connectivity strength, we expect prediction of seizure spread to be an easy task when the brain is away from this near-criticality region, but challenging otherwise.

Local linear stability analysis predicts the phase diagrams and the temporal order for node-recruitment during seizure spread

We examined whether local linear stability analysis of epileptor network models could predict the phase diagrams obtained from their simulations (Methods). In addition to the epileptor network models, we considered also a 2-dimensional reduced model version [6].

Fig 2G and 2H shows examples for three patient-specific network connectivity and a specific target EZ node. See S2-S8 Figs in S1 File for all the other target EZ nodes and patient-specific networks. Overall, despite differences in the exact locations of the phase transition curves, these predictions captured very well the qualitative structure of the phase diagrams based both on the full-dimensional (6D) model and an its 2D reduced version. Nevertheless, the local linear stability analysis appeared to overestimates the spread region in the phase diagrams, in particular near transition curve between no-spread and spread. This difference between simulation-based results and linear stability analysis prediction may result in part from the finite number (30 in this case) of stochastic realizations. At these two phase boundaries, the probability of spread events can be very small such that a very large number of realizations might be needed to observe such events.

Briefly, these phase diagram predictions were obtained as follows. We considered two different Jacobian matrices, i.e. local linearizations around a given equilibrium point of the epileptor network dynamics. One of the Jacobian matrices was derived from the network that included the active target EZ node, while the other was derived from a network where this active EZ node was replaced by a constant input to its connected surrounding nodes (see Fig 3 and Methods for details). Based on the leading eigenvalues of these two matrices, we predicted the corresponding phase of collective variable. We predicted seizure spread when the two leading eigenvalues indicated an unstable fixed point, i.e. when Re(λ) > 0 for both of the leading eigenvalues of the two matrices. (We note that while the patient-specific structural connectivity matrix is symmetric, the Jacobian matrices are not since they include also asymmetric block diagonals that reflect the intrinsic node dynamics.) The other phases of the collective variable were predicted by different combinations of the sign of the real part of these two eigenvalues (Fig 3).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Schematics for the linear stability analysis of the deterministic epileptor network.

Prediction of phase diagrams is based on the leading eigenvalues of two different Jacobian matrices evaluated at a given fixed point, and for a given global connectivity strength (w), homogeneous surrounding-node excitability (x₀), patient-specific connectivity and given active EZ node. (A,B) In this case, the Jacobian matrix is derived from a network that includes the active EZ node. The term λ_R denotes the real part of the corresponding largest eigenvalue of this Jacobian matrix. The fixed point of the system is stable for λ_R < 0, implying no seizure will happen, and unstable (seizure) otherwise. Based on this analysis alone one cannot predict the phase diagram for the seizure spread. (C,D) In this case, the Jacobian matrix (and corresponding leading eigenvalue ), are derived from an epileptor network where the active EZ node has been removed and replaced by a constant external input to nodes connected to it. This external input is related to the average seizure activity in the active EZ node (Methods). This analysis alone also cannot predict the phase diagram, in particular the phase where seizure in the active EZ node is “inhibited” by the surrounding nodes. (E,F) The joint consideration of the signs of the leading eigenvalues of the two Jacobian matrices leads to qualitatively correct phase diagrams for seizure spread as illustrated in Fig 2. Diagrams corresponds to patient P3 and node 135 as the active EZ area (x₀ = −1.6). Time delays are ignored in this linear stability analysis (Methods).

https://doi.org/10.1371/journal.pone.0272902.g003

Beyond predicting the phase diagrams for different patient-specific networks, we also assessed how well linear stability analysis predicted the temporal order in which each node in the network is recruited into seizure during the seizure spread. For that purpose, when the stability analysis indicated seizure spread (unstable fixed point), we ranked the magnitude of the (complex) components of the leading eigenvector (dominant mode) of the Jacobian matrix. (More specifically, we ranked the magnitudes for the eigenvector component related to the permittivity variable z in different nodes |v_z,i|; see Methods.) The components of the leading eigenvector indicate the directions (nodes) along which the unstable dynamics diverges. Furthermore, nodes with larger corresponding magnitudes should be recruited first. These predictions agree well with the temporal order observed in the simulations as shown in Fig 4 for a specific example and in Fig 5 for a statistical summary across patient-specific connectivity matrices and all variations in global connectivity strength, surrounding node excitability, target EZ node, and noise levels.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Relationship between the leading eigenvector of the Jacobian matrix and the temporal order of seizure spread (seizure onset times) across different nodes.

The magnitude of the components of the leading eigenvector of the Jacobian matrix from the epileptor network are denoted by |v_z,i|, with z, i indicating that the component relates to the permittivity variable z in the ith node. Similarly for , but now for the leading eigenvector computed from the Jacobian matrix where the active EZ node is replaced by a constant input to nodes connected to it (see Fig 3). (A,B) panels are based on the 2D-reduced epileptor network model while (C,D) panels are based on the original epileptor network model. Note that these magnitudes predict roughly wells the temporal order of seizure spread across the network, in particular the magnitudes based on the Jacobian matrix of the network that includes the active EZ node. The examples were obtained from a single stochastic realization for patient P1, node 64 as the active EZ (x₀ = −1.6), surrounding-node excitability x₀ = −2.13, global connectivity strength w = 1.1, and noise level at σ = 0.05. See Fig 5 for a statistical summary for different patient-specific connectivity matrices, different active EZ locations, noise levels and stochastic realizations.

https://doi.org/10.1371/journal.pone.0272902.g004

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Prediction of temporal order of seizure spread based on linear stability analysis: Statistical summary.

Each row corresponds to a different patient-specific connectivity (P1—P5 from top to bottom). Left column: the panels show the onset-time rank (rank 1 corresponds to a node which was first recruited into the seizure, and so on) and the corresponding seizure spread rank predicted by the ranked magnitudes of the leading eigenvector of the Jacobian matrix for a given epileptor network. As before, v_z,i indicates the eigenvector component related to the permitivity variable z in the ith node. Each black dot corresponds to a network with a particular active EZ node, global connectivity strength and surrounding-node excitability {w, x₀} parameter set, noise level and stochastic realization. Three noise levels were examined (σ = 0.05, σ = 0.1, and σ = 0.15). The range of (w, x₀) values were set such that seizures spread to the full network (e.g. within the yellow region in Fig 2). Thirty stochastic realizations were simulated for each parameter set, noise level and active EZ node. In each case, x₀ = −1.6 for the active EZ node. The overlapping red dots correspond to the empirical 95% confidence regions. As seen for all patient-specific connectivity matrices, true and predicted ranks tend to concentrate around the diagonal line. Middle column: the panels show the empirical probability distribution for the (integer) absolute error in the prediction of rank of recruitment during seizure spread. Right column: the panels show the mean absolute error of rank prediction (black) and the corresponding 95% confidence region (red). Prediction errors are very small for the nodes that are recruited first and last during a seizure, with the error increasing toward nodes that are recruited in between. For reference, a uniform random prediction of recruitment rank would show a mean prediction with a “U” shape, i.e. the opposite of these results.

https://doi.org/10.1371/journal.pone.0272902.g005

Further analysis revealed a roughly linear (in logarithmic scale) relationship between the magnitudes of the components of the leading eigenvector and the connectivity weights from the activate EZ node to the surrounding nodes (Fig 6A–6D). Therefore, while the prediction of the phase of the collective variable requires examination of the eigenvalues of the Jacobian matrices, prediction of the temporal ordering and spread onset times for different surrounding nodes can also be easily read out from the patient-specific connectivity matrices (Fig 6E). We also examined the relationship between spread onset time and different graph theoretical measures computed from the patient-specific connectivity matrices, including the shortest path length between the surrounding nodes and the active EZ node in a given seizure, as well as corresponding betweeness centrality and clustering coefficients. Of those, only the shortest path length showed a strong linear relationship (Fig 6F). A statistical summary across patient-specific connectivity matrices, active EZ locations, (E,w) parameter variations, and stochastic realizations is given in Fig 7. The 95% confidence regions, shown in red color, indicate that the good statistical linear relationships, observed in Fig 6 for a specific case, are overall preserved.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Relationship between the magnitude of the leading eigenvector components and the patient-specific connectivity matrix.

The magnitudes of the leading eigenvector components |v_z,i| are plotted versus the interaction weight W_i,EZ from the EZ area to the corresponding ith surrounding node. As before, v_z,i indicates the eigenvector component related to the permitivity variable z in the ith node. (A,B) panels are based on the 2D-reduced epileptor network model with the active EZ node included (A) or replaced by a constant input (B). (C,D) Similarly, but for the original epileptor network model. E Relationship between the seizure onset times for each network nodes and the corresponding interaction weight (W_i,EZ) from a given active EZ node. F Similarly, but with respect to the shortest path between a given surrounding node and the corresponding active EZ node. The panels show an example from a single stochastic realization from patient-specific network P1, with node 64 as the active EZ (x₀ = −1.6), surrounding node excitability set to x₀ = −2.13, global strength connectivity w = 1.1, and noise level with σ = 0.05. The panels show strong relationships among the leading eigenvector components, the patient-specific connectivity, seizure recruitment times, and shortest path length from the active EZ node to surrounding nodes.

https://doi.org/10.1371/journal.pone.0272902.g006

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Relationships among magnitudes of leading eigenvector components, patient-specific connectivity matrices, seizure onset times and path lengths: Statistical summary.

Same conventions as in Fig 6. Each row corresponds to a different patient-specific connectivity (P1, P2 and so on, from top to bottom). As in Fig 5, each dot corresponds to a network with a particular active EZ node, global connectivity strength and surrounding node excitability {w, x₀} parameter set, noise level and stochastic realization. Three noise levels were examined (σ = 0.05, σ = 0.1, and σ = 0.15). The range of (w, x₀) values were set such that seizures spread to the full network (e.g. yellow region in Fig 2). Thirty stochastic realizations were simulated for each parameter set, noise level and active EZ node. In each case, x₀ = −1.6 for the active EZ node. Because different values of (w, x₀) shift the seizure onset times in the target EZ node and corresponding spread onset times in the surrounding nodes, in each realization we corrected for these shifts by subtracting the averaged onset times over all nodes in that realization from onset times of each node that gives mean-corrected onset times.

https://doi.org/10.1371/journal.pone.0272902.g007

The above analysis of predictability of spread rank across nodes during full spread brings the question of how much variability there is in the spread path across different stochastic realizations of the same epileptor network under fixed parameters. We found that, typically, ranks changed by only 1 or 2 across different realizations Thus, overall, there was little variability in the actual spread path across different realizations for the noise levels used here. (The choice of input noise variance levels was constrained in this study by the requirement that the surrounding nodes do not go spontaneously into seizure states. In other words, the seizure always starts at the seizure onset zone.). S9 and S10 Figs in S1 File shows this analysis in detail for two specific examples.

Finally, we examined the power of eigenvector centrality, rather than the leading eigenvector of the Jacobian matrix, to predict spread rank across nodes. Eigenvector centrality has been used in previous studies to predict different seizure stages and propagation (e.g. [28, 29]). In Network Science, the eigenvector centrality, computed on estimated adjacency matrices, is commonly used as a complementary quantity to examine network properties. One way to estimate these adjacency matrices is to first compute a functional connectivity matrix, e.g. a matrix estimated via pairwise crosscorrelation or spectral coherence functions. Here, we assessed this approach by first estimating functional connectivity based on extrema of cross-correlation functions computed from the simulated epileptor network time series according to [30]. Statistically significant pairwise correlations contributed to edges to the derived (binary, symmetric) adjacency matrices. For this particular case of epileptor networks, we observed no successful prediction based on eigenvector centrality. Furthermore, eigenvector centrality, computed directly from the actual structural (white-matter) patient-specific connectivity matrices, lead to similar results. (This is understandable since the structural matrix does not specify by itself the location of the epileptogenic node.) S11 and S12 Figs in S1 File show the analysis and additional details for two example cases. We emphasize that these findings seem specific to epileptor network models and explored parameter regions. As stated above, the usefulness of eigenvector centrality analysis has been demonstrated before in the case of recorded electric potentials during seizures.

Seizure spread fluctuations in near-criticality regions show bimodal distribution

We examined in more detail the nature of the collective variable fluctuations in the near-criticality regions, i.e. the narrow region surrounding the transition curve between the no-spread and seizure spread phases (Fig 2). Seizure spread size showed a bimodal distribution, with modes located at the two extremes: no or small partial spread and full spread across the network (Fig 8A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Relationship between statistical fluctuations in the near-criticality region and the recruitment rank for surrounding nodes based on linear stability analysis.

Left column We examined the spread size (number of recruited surrounding nodes into the seizure; first column) when the epileptor network dynamics is in the near-criticality region (see Fig 2). Each row corresponds to a given patient-specific connectivity matrix. The empirical probability distributions show two clear modes of seizure spread in this region, one at small or no spread and the other at full network spread. These distributions are based on all of the {w, x₀} parameter set choices that resulted in fluctuations in the phase diagrams, different active EZ nodes, three noise levels (σ = 0.05, σ = 0.1, and σ = 0.15) and 30 stochastic realizations for each case. Right column When a small spread happens, for a given active EZ node the same small set of nodes tends to be recruited across different seizures. We checked whether the recruited nodes and their recruitment probability corresponded to high predicted ranks (e.g. rank 1, 2, etc) based on the ranked magnitudes of the leading eigenvector components computed from the Jacobian matrix (with the active EZ node included in the network). For example, the panel A-2 shows that the node with highest probability of recruitment (black bar) in the small spread cases also had the highest mean rank (rank close to 1, red curve). Results are plotted for one of the active EZ nodes of patients that are from top to bottom P1:61, P2:60, P3:154, P4:51, P5:154. These results show that, conditioned on knowing whether one is dealing with a small or full spread, predictability based on linear stability analysis features is still possible in the near-criticality fluctuation region, although noisier than in the full spread phase. In the high fluctuation region, linear stability analysis will fail to predict the extent of spread, but if there is spread it can predict the temporal order of seizure recruitment.

https://doi.org/10.1371/journal.pone.0272902.g008

We note that the absence of events in the intermediate spread size range is likely to reflect in part the range of explored noise levels as well as the number of stochastic realizations generated for each parameter set. In other words, the absence of events in intermediate spread size values might reflect their very low but finite probability, which could increase with higher noise levels. In any case, these results indicate a strong bimodal structure.

Finally, we verified that nodes from the same small subsets were recruited during small partial seizure spread in different realizations. Furthermore, the magnitudes of the leading eigenvector components |v_z,i| corresponding to these nodes tended to be ranked much higher than the magnitudes of those components corresponding to nodes that were never recruited in the small partial seizure spread.

These results indicate that, conditioned on knowing whether one is dealing with a small partial seizure spread or a full spread event, linear stability analysis related quantities and network connectivity still provide some predictability in these near-criticality regions. A remaining problem for future examination is whether information in preictal activity in patient-specific epileptor network models can help improve the prediction of the seizure spreading dynamics, especially near-criticallity regions where stochastic fluctuations appear to lead to substantial variability in the extent of seizure spread across the surrounding areas.