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 spacing (0.38 mm), slice number (34), field of view (240 × 240 mm²), 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.

DPARSF (date processing assistant for restingstate 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 (MathWorks, 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).

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].

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_ij = 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:
(1) $m=\frac{1}{N-1}\sum\limits_{i=1}^{N} K_{i}$ 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.

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 (Wang et al., 2017). 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.

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 N × N matrix {d_{i j}}, where the matrix of initial efficiency $\text{ }\!\!\{\!\!\text{ }{{e}_{ij}}\text{ }\!\!\}\!\!\text{ = }\!\!\{\!\!\text{ }\tfrac{1}{{{\text{d}}_{ij}}}\text{ }\!\!\}\!\!\text{ }$. 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)：
(2) ${{C}_{i}}=\text{a}{{L}_{i}}(0)(a\ge 1),i=1,2,...,N$
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:
(3) $E=\frac{1}{N(N-1)}\sum\limits_{\text{i}\ne j,i,j=1}^{N}{{{e}_{ij}}}$
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:
(4) ${{e}_{ij}}(t+1)=\begin{cases} {{e}_{ij}}(0)\frac{{{C}_{i}}}{{{L}_{i}}(t)},{{L}_{i}}(t)>{{C}_{i}}& \\ e{}_{ij}(0),{{L}_{i}}(t)\le {{C}_{i}} \end{cases}$
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.

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:
(5) ${{E}_{\operatorname{glob}}}=\frac{1}{N(N-1)}\sum\limits_{i\ne j,i,j=1}^{N}{{{\varepsilon }_{\text{ij}}}}$
and the local efficiency by:
(6) ${{E}_{\text{loc}}}=\frac{1}{N}\sum\limits_{i=1}^{N}{E({{G}_{\text{i}}})}$
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:
(7) ${{D}_{i}}=\frac{1}{N-1}\sum\limits_{i,j=1}^{N-1}{{{E}_{ij}}(0)}-{{E}_{ij}}(t)$
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.

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.

Figure 1.

Decreased efficiency of nodes after overload among GTCS subjects and controls. The horizontal axis label refers to the rest nodes after overload, the vertical axis label refers to the decreased efficiency among two groups, controls (green line) and GTCS (blue line).

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.

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.

Figure 3.

Relationship between the global efficiency, local efficiency and the tolerance parameter among GTCS and controls. The ordinate (E) gives the efficiency (global and local) among GTCS and control groups and the abscissa (a) represents the tolerance parameter, the red line indicates controls, and the blue line gives the response under GTCS.