Model for cascading failures in functional networks: application to epileptic patients with generalized tonic-clonic seizures

The dynamic process of epilepsy is modeled as a cascading failure model in functional networks derived from graph theory. The aim is to test whether cascading failure identified from functional magnetic resonance imaging data could simulate epileptic discharges in 18 subjects with generalized tonic-clonic seizure and 17 demo-graphically matched healthy controls. A cascading failure model was used to simulate the neural networks underly-ing generalized tonic-clonic seizure and healthy controls by stimulation of the node with the greatest number of connections. Results showed that the efficiency of generalized tonic-clonic seizure dropped significantly when compared to controls. Particular nodes whose efficiency altered significantly showed a correlation with the symptoms of generalized tonic-clonic seizure. Results also indicated that the left middle frontal lobe may be a potential focal area in the initiation of generalized tonic-clonic seizure.


Introduction
Generalized tonic-clonic seizure (GTCS) is a subtype of idiopathic generalized epilepsy (IGE), that exhibits loss of consciousness with tonic rigidity and clonic uncontrolled jerking as its main characteristics. Its definitive pathology is still unclear. Restingstate functional magnetic resonance imaging (fMRI) is a technique with high spatial resolution that in particular contributes to epilepsy research at the network level. An increasing number of studies have applied complex network theory to investigation of the pathology of GTCS Liu et al., 2017).
EEG-fMRI analysis has revealed that thalamo-cortical and default mode networks (DMNs) have a strong relationship with the generalized spike wave observed in IGE (Moeller et al., 2008). Morphometric study based on structural magnetic resonance imaging has documented significant reduction in cortical thickness in IGE (Bernhardt et al., 2009). Regional homogeneity analysis has demonstrated that altered regional synchronization of brain activity exists in GTCS during the interictal period (Yuan et al., 2011). Graph theoretical analysis has revealed disrupted structural and functional organization of the brain connectome in GTCS . Liu et al. (2017) implemented a method of functional network connectivity to study GTCS and reported that state-specific functional network connectivity disruptions occur in a model of GTCS.
However, as the collection of data during seizures is difficult and not always accurate, many studies focus on the network attributes of the static state to study the dynamic effects of seizures on brain networks. The theory of complex networks has been at the forefront of interdisciplinary research, along with "small-world" networks such as those proposed by Watts and Strogatz (1998) and the "scale-free" networks proposed by Barabasi and Albert (1999).
Cascading failure refers to the fact that the breakdown of a single node may have a profound effect on a whole network due to dynamic redistribution of activity flows. This progression is similar to the disease process, therefore, here it is hypothesized that a cascading failure model could simulate the dynamic process of GTCS seizures. A cascading failure model (CLM) proposed by Crucitti et al. (2004), is employed in this study. Kinney et al. (2005) has applied this model to North American power grids and found that the average efficiency decreased significantly when stimulating the node with the largest load. This node was significant for the network because of the massive information flow it received. As such, the node with the largest load could be fragile and the whole network might suffer an avalanche of alteration as it overloads. Concurrently, a new view has emerged that the pathology of epilepsy may potentially originate at such an overloaded node (Wang and Lu, 2012). Therefore, it was concluded that the node with largest load might be related to the pathology of GTCS.
This study examines functional connectivity with respect to the node with the largest load and that node's response to stimulation in both GTCS subjects and healthy controls. This is followed by investigation of overloaded nodes and how their efficiency was altered with respect to GTCS. The data processing assistant for resting-state fMRI software (DPARSF) is employed to initially process group GTCS and control data. An automated anotomical labeling (AAL) atlas is then used to divide the whole brain into 90 nodes and a Pearson correlation coefficient calculated for all node pairs. Subsequently, a threshold is set to structure a binary logic network that represented the functional connectivity of the whole brain. Finally, a cascading failure model is applied to stimulate the node with the largest load in GTCS and alterations of whole brain were analyzed during GTCS.

Subjects
Subjects were diagnosed with GTCS according to the criteria of the International League Against Epilepsy (ILAE) classification: (1) Typical clinical symptoms of generalized tonic-clonic seizures, including tonic extension of the limbs, followed by a clonic phase of rhythmic jerking of the extremities, and loss of consciousness during seizures; (2) Generalized polyspike-wave in the interictal scalp EEG; and (3) No focal abnormality in the structural MRI. Subjects had no history of either neurological or psychiatric disorder. All subjects provided written informed consent prior to the study which was approved by the medical ethics committee of the Second Hospital of Lanzhou University in accordance with the Declaration of Helsinki and good clinical practice guidelines.
Resting-state fMRI data was acquired with a 3-T MRI (Siemens, Germany). Functional imaging was acquired by use of a gradient echo-echo planar imaging sequence. The scan range included the whole brain and a baseline parallel to an anterior and posterior orientation. Scanning parameters were: time to repeat (2000 ms), echo time (30 ms), slice thickness (3.8 mm), slice spac-ing (0.38 mm), slice number (34), field of view (240 × 240 mm 2 ), matrix (64 × 64), fractional anisotropy (90 • ), time points (200), and scan time (six minutes forty seconds).
Subject parameters are given in Table 1. A total of 20 patients were included in this study but two patients with head motion range greater than one millimeter were excluded. No controls were excluded. There were nine female and nine male rest GTCS. Their age ranged from 7 to 46 years, (mean 26.4 years). Duration of epilepsy ranged from one to thirty years (mean duration 6.2 years). All of subjects were right-handed. There were 17 normal healthy controls, eight male and nine female. All controls were right-handed and aged from 15 to 37 years (mean age 25.8 years). The ages of both GTCS subjects and controls were normally distributed, (standard deviation 6.9 and 6.2 years, respectively). The age distributions of the two groups were not significantly different.

Data preprocessing
DPARSF (date processing assistant for resting-state fMRI) (Yan and Zang, 2010) and REST (resting-state fMRI data analysis toolkit) (Song et al., 2011) were employed to preprocess fMRI data with a software implementation in MATLAB (Math-Works, Natick, MA, USA). DPARSF processing included: format conversion from DICOM (digital imaging communication in medical) to NIFTI (neuroimaging informatics technology initiative), removal of the first 10 time points, slice timing, realignment (exclusion criteria: 1 mm and 1 • ), normalization (EPI template), smoothing (full width half maximum: 6 mm), detrending, and filtering (0.01-0.08 Hz). The processes with REST included covariable extraction and regressing out covariants (gray matter, white matter, and cerebrospinal fluid).

Functional connectivity
The functional connectivity was instructed by the REST. The AAL (automated anatomical labeling) brain atlas (Sporns et al., 2005) was employed to define the regions of interest (ROI). The whole brain was divided into 90 nodes and the time course of every node extracted. The function ROI − wise from the REST software was employed to compute the functional connectivities between the 90 nodes. A 90 × 90 matrix was obtained from this calculation. Matrix values ranged over [-1, 1].

Construction of the brain network derived from graph theory
The brain network refers to a collection of several brain areas and edges (connections) in graph theory . It consists of brain areas and connections between them, which also indicate the functional network topology of brain regions. Here, the brain network was represented by an adjacency matrix G which representing the brain network. In this matrix, G i j = 1 describes an edge between the ith and jth nodes.
The Pearson correlation coefficient matrices of GTCS subjects and controls were computed individually. Initially, the Z-score of each GTCS and control data matrix was computed to improve the normality of the distribution of their correlation coefficients. An average correlation coefficient matrix of all subjects was then taken between groups and a threshold (t) binary network constructed (t > 0) (the threshold represented a critical value indicating the presence of edges in brain regions). A Z-score matrix was employed to construct the logic network. Values in the GTCS and control matrices were set to 1 if the absolute value of the corresponding value was greater than threshold (t), otherwise matrix values were set to zero. Diagonal matrix values were set to zero to self-join. Adjacency matrices were then obtained that represented the GTCS and control networks following their computation. According to (Achard and Bullmore, 2007), the brain network was complex and sparse. When constructing the brain networks two rules were followed: One was that there should be no isolated nodes (Wang and Lu, 2012); the other that network density should range from 10% to 50% (Xue et al., 2010). Network density was obtained from: where N is the sum of nodes and K i is the degree of the ith node. The topology of the GTCS network is then analyzed to find the node with the largest load and further examine its pathological basis.

Analysis of functional connectivity of node with the largest load
To assess whether the node with the largest load was abnormal, it was initially selected as a region of interest (ROI) (Wang et al., 2012). The resting-state averaged time series from the ROI were then calculated and correlated with the time series of the remaining voxels in the whole brain to give a correlation map. To improve normality this result was converted to a z map by using Fisher's r-to-z transformation . The GTCS and control groups were examined and compared for differences by two-sample t-test. Only results significant at p ≤ 0.001 and cluster size greater than 15 voxels were considered (Ke et al., 2017). The cascading failure model is then applied to explore the networks of GTCS and control subjects.

Cascading failure model
The network was represented by a weighted and undirected graph G with N nodes, where G was given by an N × N matrix {e i j }. To define the initial efficiency e i j , the shortest path d i j was computed between the ith and jth nodes according to the adjacency matrix described above. The matrix containing the shortest paths was represented by the This model assumed that the communication between generic node pairs were connected by the most efficient path (Crucitti et al., 2004). Thus, the load L i (t) with respect to the ith node at time t was defined as the total number of the most efficient paths passing through the given node at the given time. Each node was characterized by a capacity defined as the maximum load a given node could accommodate. The capacity C i of the ith node was proportional to its initial load L i (0): where a gives the tolerance parameter of the network (every network conclued a tolerance parameter which represented the degree of protection itself). The average efficiency represented the mechanism of transportation in this model, and was given by: It was assumed that at time 0, the initial efficiency matrix was {e i j } and one node was overloaded (removed). The removal of a node changed the most efficient paths between nodes and caused a redistribution of the loads, thus created overloads at some nodes. Therefore, the efficiency at time t was obtained from the following iterative rule: where j extended to all first neighbors of i. In this rule, if the ith node was congested at time t, the efficiency of all the trajectorys passing through it would be reduced, so that the new paths which were most efficient would probably eventually alter. This may trigger a further distribution in a repetitive process unless the new paths were the most efficient, thus stabilized.

Analysis of network characteristics
When a network is overloaded by cascading failure, its topology is likely to alter as both local and the global efficiency are determined by this topology. Simultaneously, the decreased efficiency of nodes after overload provided significant insight into the topology. The global efficiency was given by: and the local efficiency by: where N represented the sum of nodes, E(G i ) represented the local efficiency of the ith node. E(G i ) was defined as the efficiency of the ith node among its all first neighbors. The decreased efficiency of nodes was given by: where E i j (0) represented the initial efficiency of the i jth connection, and E i j (t) represented the efficiency after the system had relaxed to a stationary state.

Brain network construction
The threshold was set within the range [0.2, 0.5] and increased by 0.01 at each timestep when the network density was computed. The relationship between network density (m) and threshold (t) was reported as shown in Fig. 1. In this figure, the density of the network in GTCS is shown. We found that isolated nodes only existed when the threshold was raised above 0.37 for GTCS, but it could be as low as 0.31 for the controls. No node should be isolated within the network and network density should be within the range [10%, 50%] when constructing both networks. The network of healthy controls did not meet the above criteria as the threshold was larger than 0.3, So this threshold was considered as a meaningful value. Therefore, the threshold were set equal to 0.3 when constructing both networks as network density (computed by Eqn. 1) was 17.6% in GTCS and 16.5% for controls, in the absence of isolated nodes. Consequently, a pair of unweighted and undirected adjacency matrices could represent the networks of GTCS and controls.

Analysis of functional connectivity
The total time of every node contacted by all the most efficient paths was computed and the largest load (the most times) for GTCS was found at the left middle frontal node. It was considered that this result might have a pathological basis, so this node was selected as the seed node in an examination of the functional connectivity with other brain areas in GTCS. Differences were then examined by the two-sample t-test and only results significant at p ≤ 0.001 for cluster sizes greater than 15 voxels were considered. The brain areas which altered significantly are given in Table 2.
Here, it is shown that there were four areas which decreased significantly when compared with controls, the lingual, the left middle occipital area, the left lower frontal triangle area, and left middle occipital area. Hence, the functional connectivity of the seed node (left middle frontal) demonstrated an abnormality when compared with controls.

Simulations
The left middle frontal node in GTCS indicated an abnormality in functional connectivity when compared with healthy controls. To determine how this brain area influenced the rest brain areas under the cascading failure model, the left middle frontal node was selected for overload in GTCS and controls with tolerance parameters. Results are shown in Fig. 2. Initial efficiency refers to efficiency prior to network overload. The global efficiency in GTCS and controls was 0.5354 and 0.5077, respectively, while the local efficiency in GTCS and controls was 0.7560 and 0.7523 respectively. This shows that the initial efficiency in GTCS was greater than that for controls. The global efficiency of GTCS and controls increased as the tolerance parameter was increased, and both eventually approach their initial efficiency. Similarly, the local efficiency of GTCS and controls showed positive correlation with the tolerance parameter. This showed local efficiency to be greater than global efficiency for both GTCS and controls. Only when the tolerance parameter reached 1.2 was the global efficiency returned to its initial value in GTCS. However, controls reached 1.02 when returned to its initial efficiency. Overall, the network of controls was vulnerable to reach stability compared with GTCS.

Figure 2. Relationship between threshold (t) and network density (m). Red and blue lines indicate GTCS and control subjects, respectively.
Both networks were overloaded under different tolerance parameter values, so the global efficiency in GTCS and controls would accordingly be decreased under these changes. This resulted in the mean global efficiency of every node being altered. Fig. 3 shows that alteration of controls were within the range [0, 0.05] (mean value 0.0022) and the alteration of GTCS was within the range [0.19, 0.5] (mean value 0.3095). Overall, the change in GTCS was larger than that of controls. The trend of controls was to be less variable than that of GTCS. The encephalic regions are specifically reported to exhibit greater alteration in GTCS (see also Table 3). In Table 3 it is shown that the node with the largest alteration was the right rolandic operculum, whose alteration of effiency reached 0.4702.

Discussion
The shortest path was found to be relevant to transport properties in a complex network. However some nodes in the shortest path may be heavily loaded and jammed by considerable activity. Goh et al. (2001) defined the total time taken for the most efficient paths to pass activity through a as "load" (Goh et al., 2001). Here, the node with the largest load (left middle frontal) was vital to the network because of the information it transmitted. This node could severely impact the whole network in the cascading failure model similarly to the disease process of epilepsy. Therefore, it was hypothesized that the cascading failure model studied here could simulate the dynamic process of seizure in GTCS where the node with the largest load was the potential focal area. The functional connectivity of left middle frontal node was examined along with other brain areas in GTCS. Results demonstrated abnormality in functional connectivity between it and other brain areas when compared with the controls, which suggested that this brain area might reflect a pathophysiological mechanism in GTCS. The fourth largest load in healthy controls was observed for the left middle frontal area, which indicates that this brain area was also a vital area for the controls. Based on the hypothesis, the left middle frontal area was overloaded by use of the cascading failure model in GTCS and controls.
Brain networks were constructed to model GTCS and control groups. It was found that the network density of GTCS and controls was decreased with increased threshold. This indicates that the threshold is of significance for GTCS and control networks. Furthermore, the selection of threshold may affect subsequent network behavior. Generally, the density of a brain network model should be within the range [10%, 50%], a requirement that was also the main reason for selecting the threshold. Consequently, the threshold was set to 0.3 when constructing both networks. Whatever the threshold was, it should ensure that network density is maintained within the range [10%, 50%] when constructing a brain network model. Latora and Marchiori (2001) showed that the local efficiency determined the extent of fault tolerance in a system, thus this measure showed how efficient the communication was between the first neighbors of the ith node following its removal. Alternatively, the global efficiency focused on the entire network, so a few nodes with extremely slow connections did not mean that the efficiency of the entire network was diminished (Latora and Marchiori, 2001). Here, it is noted that both the global and local efficiencies in GTCS showed a positive correlation with the tolerance parameter. This indicates that the tolerance parameter was conductive to the robustness of the GTCS network, as Crucitti et al. (2004) reported that the tolerance parameter was a realistic assumption in the design of an infrastructure network but was limited by the cost. The reason for the tolerance parameter enhancement of the robustness of GTCS was that it enhanced node capacity. Here, the capacity of every node was decided by the shortest (most efficient) paths. The larger the capacity of a node, the shorter the paths that contacted it. The shortest paths of nodes represented the functional information flow. An enhanced tolerance parameter meant the functional connectivity of nodes was strengthened. Hence, stronger functional connectivity supported greater robustness in GTCS. Studies report that patients with GTCS suffer longterm neuropsychological impairment (Vlooswijk et al., 2010).
A previous study (Ling et al., 2016) had also found a negative correlation between duration and reduced long range functional connectivity density (FCD) in GTCS. This indicates that longterm functional impairments decrease the functional connectivity of some nodes. Meanwhile, functional impairments would have a direct effect on exchange and consequently decrease robustness of the whole network. These reports indirectly provided evidence that functional connectivity had a positive correlation with robustness in GTCS. When functional connectivity was strengthened it could protect the whole network from breakdown. Nevertheless, functional connectivity could not be infinitely large because the brain network is a sparse network (the construction of the network was based on the functional connectivity). The result of global and local efficiencies in controls indicate that the networks of controls exhibited stronger robustness when compared to that of GTCS. Especially, when the network of healthy controls maintained a stationary state, GTCS produced a large drop in efficiency when the tolerance parameter ranged over [1.02, 1.12]. The reason for this might be that the left middle frontal node had the largest load in GTCS, so this node had a stronger relationship with other nodes. The networks of controls showed a higher robustness when compared with that of GTCS. This result also indicated that the left middle frontal gyrus had an essential position in the GTCS network.
The decreased efficiency of every node in GTCS was much larger than for controls. It was noticed that the trend of alteration in controls was more uniform than for alterations in GTCS. The alteration was determined by an iterative efficiency rule and the changing node load. Small changes in efficiency meant that the new node load was barely altered. This suggests that removal of the left middle frontal area has little effect on control network topology. As Crucitti et al. (2004) reports information might flow through another path as one node is overloaded (removed). This was the case for the GTCS network, where the information flow chose alternative paths to complete functional connectivity when the left middle frontal area was removed. Consequently, node load might alter greatly as the efficiency of nodes greatly reduced. In conclusion the topology of networks in GTCS altered significantly as a consequence of overloading the left middle frontal gyrus. In a previous study, Schindler et al. (2008) analyzed seizures from the EEGs of a large patient group and found that the changing functional network topology during seizures was accompanied by an initially decreased stability of the global synchronized state. The results reported here are in accordance with Schindler et al. (2008). The alteration of topology may have a deep influence on the function of nodes in GTCS.
It was found in this study that the efficiency of some brain areas in GTCS were extremely decreased compared with their initial efficiency. These brain areas included the right back cingulum, right precuneus, right rolandic operculum, right side of tongue, upper and lower right occipital areas and right insula. The right posterior cingulum and right precuneus belong to the DMN. The posterior cingulum especially is the core of a DMN (Dum and Strick, 2002;Salvador et al., 2005), an area that generated a drop of 0.3801 in its efficiency. The literature reports that the DMN may have a strong relationship with the loss of consciousness during seizure (Blumenfeld, 2012;Sporns et al., 2008). Both nodes in the DMN which decreased greatly might provide new evidence for the loss of consciousness during GTCS seizures.
It was also found that the right rolandic operculum, right side of the tongue, and upper and lower right occipital areas belong to the motor areas. These regions, which produced a large drop in their initial efficiency, enabled the functions of motor disorder. This may have a direct association with voluntary movement and generalized tonic clonic activity. Concurrently, Gotman et al. (2005) has shown that the insula is correlated with an epileptic attack, which may produce impairment with respect to the insula. So this may indirectly explain why the efficiency of right insula decreased significantly. Overall, the mean drop in efficiency was 0.3095 in GTCS. This suggests that the whole brain suffered serious damage and most functions would cease to work naturally. Such results show that the cascading failure model was significant in the study of GTCS pathology. Moreover, the left middle frontal area was chosen as the overloaded node in the cascading failure model to simulate the disease process. The pathology characteristics indicated that the left middle frontal area may be a potential focal area in the disease process of GTCS.

Conclusions
In this study, it was found that the left middle frontal area was the node in a model of GTCS with the largest loads and that this area also exhibited a greater load in control subjects. Thus, it is considered that this node might be the pathogenic brain region. For this reason the node was selected as the seed node in a model network developed to examine its functional connectivity with other brain areas. This node was found to exhibit abnormal functional connectivity when compared with controls. The left middle frontal area was pathological and it was hypothesized that the cascading failure model described in this paper was suited to simulate the dynamic process of seizure in GTCS. The left middle frontal area was oversimulated to induce GTCS in controls by employing cascading failure to simulate the disease process. Tthe efficiency of the network was altered and the function of nodes with larger drops of efficiency were analyzed. It was found that a cascading failure analysis was of great significance in studying the pathology and it is concluded that the left middle frontal may be the focal area in the disease process of GTCS.