Clinically localized seizure focus maybe not exactly the position of abating seizures: a computational evidence

By modeling the brain as a network, the challenge of abating seizure can be recast as a problem of network control. In the premise of bringing network under control, the minimum number of nodes prerequisite for controlling seizures are thus a natural aim of interest, which is still an outstanding issue. Here, we use the network structural control theory to guide the selection for the optimal control nodes with the aim of fully abating seizures. Firstly, we construct the dynamical complex network of pathological seizure by estimating the synchronicity and directionality of information flows over time between EEG signals from 10 patients with focal epilepsy. Then, based on the controllability and observability principles of complex systems, the minimum key nodes which are effective to fully control the network seizure behaviors are obtained. Results show that the calculated control nodes are distinct with the focus zones from clinic report. This suggests that the full control of epileptic network may not only related to the focus zone, the other non-focus nodes could also play important roles. This finding is validated by using the spatiotemporal neural network model connected with our modeled dynamical adjacent matrix. It successfully reproduces the original EEG signals which can be effectively abated by applying pulse stimulation on the identified key nodes or resecting them, while the partial effects can be obtained when functioning onto the clinically identified focus zones of 60% patients. Interestingly, for another 30% patients lesser nodes than clinic reports are need to fully control seizures. In addition, our work facilitates to identify the evolution paths of information flows, so the non-clinic focus zones identified through controllability principle can be supposed to be the potential seizure foci. In sum, our work propose a general methodology or strategy for seizures focus localizations that could comprehensively consider the time-evolving information flow and the analysis of controllability mechanisms driven by the real seizures data, as well as the computational validation. This may promote to develop the canonical computational framework with the core intent of providing support for clinical treatment decisions.


Introduction
Epilepsy is a chronic disorder which is considered a change in brain state (interictal, preictal, ictal and postictal) [1]. Applications include localization [2], detection [3] and prediction [4] of epileptic activities have been broadly explored and viewed as the most challenging issues.
In principle, the aberrant brain activity of neurological diseases [5,6] can be restored through electrical stimulation. It is currently believed that electric stimuli can cause specific state-dependent transitions [7,8], which have been considered as the potential supplemental approach to abate seizures [9][10][11].
Neocortical focal seizures are thought to be arised from the spatially localised pathological tissue of brain, i.e., epileptogenic zone (EZ) [12][13][14]. Surgical removal [15,16] and electrical stimulation [17] of EZ are often the final treatment options. However, they, in general, achieve only partially successful results. This may due to that the most widely used techniques for localizing the EZ in patients with focal epilepsy is only the visual identification of raw data, lacking the spatio-temporal correlations of different seizure regions which might mask the true EZ and fail to prevent future seizures [18].
If the distinct brain regions are represented as nodes and the statistical dependencies between regions are represented as edges, the brain can be modelled as a network. Hence, the challenge of suppressing seizures can be recast as a problem of finding key nodes which can regulate the whole seizure network [19]. Note that, seizures are characterized by rapid state transitions [20] and the information flow propagation through the brain may depend on the evolutions of the coherent state between brain regions [21]. This suggests that the seizures network should be modelled as a timeevolving dynamics effective network which can capture the dynamic influences over changing patterns of network coherence [22].
However, the graph measures (e.g., degree and betweenness centrality) found to correlate with epileptogenic focus [23,24] didn't give the controllability analysis of the seizures network, thus providing weak understanding of the epileptogenic network mechanisms. By means of the graph theory, [25,26] revealed the alterations of functional connectivity of the left and right temporal lobe epilepsy and the focal epilepsy, respectively. But they paid less attention on the dynamic evolutions of information flow over time. Based on mathematical model, state-feedback-based controllers [27] with the assumption that all states are known over time was introduced to suppress epileptic seizures. But it is hard to validate the effect of the controller in clinic due to that some of the state may not be measured once their physical meaning is not fully established according to current knowledge on this model. Therefore, there still lacks a canonical computational framework for seizures focus localizations that could comprehensively consider the timeevolving information flow and the analysis of controllability mechanisms driven by the real seizures data, as well as the model validation. Hence, our aim is to propose a general methodology or strategy for seizures focus localizations that can indicate the key regions with the intent of providing support for clinical treatment decisions.
Evidences have shown that epileptic seizures are generated by the collective activity of extended neuronal networks. They are often referred to as pathologically hypersynchronized states [28,29]. Synchronization measures have also been applied to predict and localize the epileptic activity [30][31][32]. This means seizure control may be induced by modulating synchronization level. However, the synchronization between different brain areas, in particular, the synchronous coupling direction between signals, cannot be measured directly. To better understand seizure onset and propagation, as well as propose an effective seizure interfering strategy, it is needed to identify the coherent evolution and transferring process of the dynamic signal within the epileptogenic network by applying an appropriate analysis method.
Based on the above discussion, firstly, the timeevolving dynamic effective network will be constructed by employing a methodology to investigate the synchronicity and directionality (leading role) over time between signals from different brain areas. Then, based on the constructed network, we will consider the network mechanisms of seizures, with which to identify the key nodes of influencing the network state. In the premise of bringing system under control, the minimum number of nodes, prerequisite for controlling seizures are thus a natural aim of interest. Also, the macroscopic neural population model network dynamics will be used to compare with the clinical EEG data. Finally, the effect of stimulation on seizure abatement will be computationally explored.

Clinical data
In this paper, the SEEG data were recorded from more than 100 channels of ten patients who had been admitted to the hospital with refractory focal epilepsy at Sanbo Brain Hospital of Capital Medical University in Beijing. This study protocol was approved by the ethics committee of Sanbo Brain Hospital of Capital Medical University and the subject was written informed consent. Due to the hospitals regulations, the data cannot be made publicly available. Note that, in this paper, we mainly use the seizure data of the first patient (see Table 1) to perform the analysis, validation and visualization of our results. Part of the time series of the signals from 10 channels of the first patient are shown in Fig. 4. Detailed description was exhibited previously in Yang et al. [33,34]. The 10 channels are denoted by F08, M09, C09, D03, D04, G05, H10, J07, K01, K02, where F: the middle frontal gyrus-medial frontal lobe, M: the middle frontal gyrus-outer side of frontal lobe, C: the front of the forehead back triangle area-frontal island cover-buckle belt back front lower part, D: the posterior central anterior sulcus, G: the back of the middle frontal gyrus-the inner side of the frontal lobe, H : the anterior central gyrus-cingulate gyrus, J : the front central back to upper middle, K : the central posterior middle-cingulate gyrus, respectively. The clinical EEG showed that two clinical seizures of epilepsy were recorded during the monitoring, which were mainly manifested as complex movements. At the beginning, the F08, M09 and C09 displayed low-amplitude and high-frequency discharges. After 0.5-1s later, the discharges spread to D(03-04), then to G05, H10, J07 and finally to K(01-02), which successively showed the transitions from the low amplitude and fast activity to spike-slow wave discharges.

Definition of event
For a time series x(n) (n = 1, 2, ..., T ) with total T steps, the events can be defined as the local maxima by the following law [35,36]: where the event times are denoted by t k (k = 1, ..., n e ) with n e being the total event number; M and h are the parameters which have different values for the EEG data from epileptic animal model and epileptic patient.

Coupling direction and synchronization measurement
Consider two time series x i (n) and x j (n) with equal T total steps, where the event times are t i r and t j s (r = 1, ..., n i ; s = 1, ..., n j ), respectively. Allowing a time delay τ (smaller than half the minimum interevent distance of these two time series), we can define a quasi-simultaneous appearances of events, which can be described by Eq. (5). Then, e τ (x i |x j ) (e τ (x j |x i )) is used to denote the number of times an event appears in x i (x j ) shortly after it appears in x j (x i ). Their symmetrical and anti-symmetrical combinations expressed by Eqs. (2) and (3) can be used to describe the event total synchronization level and event delay pattern, respectively. Thus, Q τ and q τ are normalized to [0,1] and [−1, 1], respectively. Here we suppose that the delay behavior suggests the directions of information flow between x i and x j signals, i.e., as a tool to identify the leading role between them. That is we can identify the driver-response relationship between them. If q τ in Eq. (3) is larger than 0, x i can be considered as the reason which drive the propagations of pathological information to x j . where To present the on-line changes in event synchronization and event delay pattern (leading role), we give the time resolved variants of Q τ (n) and q τ (n), which can be obtained by simply modifying the conditions of Eq. (4), i.e., where the n is no bigger than T . Hence, where Q τ i j (n) increases one step only if an event is found both in x i (n) and x j (n) within the window τ , q τ i j (n) takes one step up (down) every time an event in x i (n) shortly ( i.e., also within the window τ ) precedes (after) one in x j (n).
To further observe the instantaneous changes in event synchronization and event delay pattern ( directionality, i.e., the leading role of information flow between two signals) at time n, a window technique was applied to calculate them averaged over the slide window n = 5120 time steps, expressed as follows:

Network measurement
Epilepsy has been understood to be the network disorder, but how to capture the evolutions of pathological information over time and localize the focus remain a hard issue. Thus, we need to reconstruct the epileptic network. To measure the network synchronization, the above equations can be modified as follows, √ n i n j (13) where i = j and i, j ∈ {1, 2, ..., N }. The directionality of synchronous information flow, i.e., drive-response relationship of any two signals within network can be represented by the directed connectivity matrix: where γ =1000 is the amplifying parameter, Here our analysis assumes a loop free graph, thus a ii = 0, for 1≤ i ≤ N .

Controllability and observability
Due to the ubiquity of complex networks, it is highly desirable to be able to apply proper control to guide the network dynamics toward states with the best performance. In fact, the brain system dynamics can be described by the generic state-space form: where x(t) = (x 1 (t), ..., x N (t)) T ∈ N captures the state of a system network of N nodes at time t.
Here, each node corresponds to every contact of EEG channel which records the electrode activity level of neuronal population from separate local brain region. u(t) = (u 1 (t), ..., u M (t)) T ∈ M is the control input which captures the influence of the environment (e.g., external stimulation pulse input). Hence, the system's internal state depends not only the time t and the system's internal state x(t), but also the external input u(t). The field vector is described as f (t, x(t), u(t)). y(t) = (y 1 (t), ..., y K (t)) T ∈ K is the output. To observe the system we need to monitor a set of variables of y(t) ∈ K . They also depend on the time t, the system's internal state x(t) and the external input u(t), which are represented as their function h(t, x(t), u(t)). Even though most real systems are driven by nonlinear processes, the controllability [37,38] of nonlinear systems is in many aspects structurally similar to that of linear systems [39]. Thus, the study of this nonlinear system can be reduced to a canonical linear timeinvariant dynamic system: (17) A N ×N is the adjacency matrix and N is the number of nodes. M is the number of driver nodes, B is the N × M input matrix, and u(t) is the time-dependent input. C is the K × N output matrix. According to Kalman's controllability rank condition [40,41], the above system can be controlled from any initial state to any desired state in finite time, if the N × N M controllability matrix has full rank, that is rank(J 1 ) = N . In addition, it is needed to estimate and reconstruct the system's internal state from its outputs, i.e., the observability [42,43] of complex system. Of course, the simultaneous measurement of all internal variables can completely observe the system's state, but, in practice, the complete description of a complex system can only be limited to a subset of variables (termed sensors in [42]). According to the previous method [42], here we will also adopt a graphical approach to determine the sensors that are necessary to reconstruct the full internal state of a complex system. The observability criteria can be formulated algebraically for linear dynamical systems, in particular, the system (16) is observable if the Jacobian matrix has full rank, that is rank(J 2 ) = N

Neural field model
The model we used here is modified from Wendling et al. [44]. Researches show that the collective behavior of neural populations leads to epileptic activities and produces bursts of spikes on the EEG. This model, followed the lumped-parameter approach in which the neural populations are modeled as nonlinear oscillators, is particularly suited to model the SEEG signals.
The single-compartment model includes two subsets of neurons, namely the main pyramidal cells and the inhibitory interneurons. The corresponding block diagram representation is shown in Fig. 1. h e and h i , viewed as the static gains of filters, describe the impulse Note that, each linear transfer function h e or h i can introduce a second-order ordinary differential equation. Hence, according to the block diagram, the following set of 3 second-order differential equations that describe the state of every population can be established: where a = 100s −1 and b = 50s −1 are the membrane average time constant and dendritic tree average time delays, respectively; A and B = 35 mV are the average excitatory and inhibitory synaptic gain, respectively.
), e 0 = 5s −1 is the maximum firing rate, ν 0 = 6 mV is the post-synaptic potential corresponding to the half of e 0 , r = 0.56 mV −1 is the steepness of the sigmoid. In the model, the noise current is modelled as where μ is the bias current, σ is the noise intensity and ξ(t) is the Gaussian white noise with zero mean and unit variance. K i j is introduced to define the coupling gain constant from populations j to population i, where a impulse response h d , viewed as a filter, is introduced to model the delay connection from population i. Hence, the two-compartment coupled model is governed by 4 second-order differential equations describing the inter-population behavior: where x 4 is the output of the EPSP transfer function h d , a d = 33s −1 is the average time delay on efferent connection from a population.

Stimulation setup
It is shown that external stimulus and electromagnetic induction as well as the noise can disturb the neuron population activity [45][46][47]. In particular, transcranial random noise stimulation (tRNS) [48][49][50][51] involving the application of random noise oscillations to the selected brain regions was shown to be a painless and unharmful neuromodulation method. It can increase the neuronal excitability via stochastic resonance and hence modulate the cortical plasticity. Due to the advantage of noise stimulation, in this paper, we employ the normally distributed noise stimulation with cathodic phase pulses (cNDNS) to computationally explore the seizure abatement effect. During the simulation, normrnd(μ, σ, 1, time/dt), i.e., a Matlab function with μ = 75 and σ = 25 is used, where time is the simulation duration, dt is the integration step size. The reason why choose cNDNs stimulation is also because that in Eq. (20) or (21) p(t) globally measures the influence of neighboring or distant populations which was assumed to be Gaussian white noise. The Gaussian white noise had been computationally evidenced to play roles in the seizure generations. The cathodic phase pulses stimulation of cNDNs is intent to counteract the effect between two different populations.
The standard fourth-order Runge-Kutta integration scheme under the MATLAB (MathWorks, USA) simulating environment was used with Eqs. (20)-(21). Initial conditions were set to zero in all simulations, and an integration step size of ∼ 4 ms (corresponding to the clinic sampling frequency: 256 Hz) was used. Simulation duration is 100s for Fig. 2a, 20s for Fig. 2c, d, 136.25s for Fig. 10.

Synchronization and leading role between two coupled EEG signals
The above mentioned statistical analysis for estimating synchronization and delay pattern [35] was directly applied into the analysis for EEG data without the verification of effectiveness. However, this is more or less a shortcoming. Hence, in this paper we first verify the effectiveness of this method using a coupled neural system, as shown in Fig. 1. Figure 2b shows that the firing of the system typically depends on the average excitatory synaptic gain A. It is seen that when A is larger than A ≈ 4.5 mV the system can generate spikes from background state, the dominant frequency of which can successively increase with A progressively increasing. As illustrations, in the middle and upper panels of Fig. 2a, the system shows the spiking state and background state with A = 4.85 mV and 3.25 mV, denoted by system 1 and 2, respectively. However, as shown in the lower panel of Fig. 2a, the coupling function (K 12 = 200) from system 1 makes Further, we explore the synchronization level and the leading role between the two unidirectional coupling systems using the above statistic method for measuring synchronization (Q τ ) and time delay pattern (q τ ). It is clearly observed that when K 12 =200 the two coupled systems have high level of synchronization (Fig. 2d, Q τ ≈ 0.9 ) and the system 1's leading role obviously dominates (Fig. 2c, q τ ≈ 0.9). In particular, this method can consistently, steadily and subtly capture the increasing evolutions of time delay pattern (i.e., the leading role of system 1) and synchronization by gradually increasing the K 12 from 0 to 300. However, for the surrogate data which is obtained by disordering the discrete time series of system 2 (lower panel in Fig. 2a, it is seen from the green lines in Fig. 2c, d that as expected the drive-response rela-tionship with system 1 disappears and desynchronizes (a background 'zero' level) between them. Therefore, these results verifies the reliability of the above mentioned statistical analysis [35] in determining the real synchronous level and drive-response relationship of two signals.
Then two original signals from channels D03 and D04 (Fig. 3a) are used to analyze their synchronous level and delay patterns (i.e., leading role). In particular, we aim to observe their evolutions over time, i.e., time-resolved event synchronization (Q τ (n)) and delay pattern (q τ (n)), where n denotes datapoint (total 40,000 datapoints with sampling rate 256 Hz used in this paper, the seizure onset is labeled at n=27880 denoted by vertical dashed lines). It is seen from Fig. 3b that with time going the synchronization between the signals of channels D03 and D04 gradually increases. Along with that, as shown in Fig. 3d, the leading role (event delay patterns) of D04 for D3 also increasingly dominates. Remarkably, after the seizure starts, their synchronization level and the leading role of D04 are transiently, quickly and significantly increased. This suggests the quick synchronous information flow from D04 to D03, which is consistent with the previous clinical results that pointed out that increasing synchronization may promote seizure termination. The corresponding instantaneous time-resolved event synchronization (d Q τ (n)) and delay pattern (dq τ (n)) can be observed in Fig. 3c, e, where there are evident humps at the beginning of seizures. Then, the synchronization increasing become weak. And the leading role of D04 also weakens a little, i.e., the slight decreasing in q τ (n) which corresponds to the dq τ (n) < 0.
These results are credible when compared to the pairs of D04 signal and the surrogate signal of D03, where the surrogate data ( the usefulness of such surrogates was discussed in more detail in [35]) are obtained by shifting the signals 2000 data points (∼7.8 s) to the left without changing the individual properties of each signal. The result of surrogate pairs corresponds to the cyan lines which reach to lower level and even the background 'zero' level for the synchronization level and delay patterns.  Table 1).

Network analysis based on the clinical EEG signals
Next, we aim to explore the network evolutions of event synchronization and delay patterns by reconstructing the directional network based on the analysis of clinical EEG signals. Figure 4 displays the original EEG data with 40000 datapoints from 10 channels, i.e., F08, M09, C09, D03, D04, G05, H10, J07, K01, K02, respectively. These ten channels are mainly divided into four groups (see Fig. 7) according to the clinical observations, i.e., group I: F08, M09, C09; group II: D03 and D04; group III: G05, H10, J07; group IV: K01 and K02. Epileptiform discharges are observed to spread from group I to group IV. Figure 5a shows the time-resolved event delay pattern evolutions between G05 and the other 9 channels. Original EEG data recorded from 10 channels of the first patient (see Table 1 The leading roles of G05 to other 9 channels can be determined if q τ (n) > 0. Note that, to avoid the hidden wrong outcomes induced by the noise disturbance, we set an event delay threshold q τ (n) = 10 that can define the significant leading role of G05 when q τ (n) > 10 (the same hereinafter). In particular, after seizures start (denoted by the vertical dashed lines), we can observe the evident driving roles of G05 to M09, C09, D04, K01 and K02, respectively, in seizure spreading. Figure 5b give the corresponding instantaneous changing rate of q τ (n), i.e., dq τ (n). Figure 5c shows the corresponding time-resolved event synchronization evolutions of each pair of signals consisting of G05 and each other channel, as well as the whole network. As a whole, the network presents the higher synchronization level while seizure occurs. As shown in Fig. 5d, the instantaneous rate is also high, meaning the very active synchronous activity.
The other similar scenarios can be found in Fig. 6, where we provide rich drive-response relationships between each channel and every other channels. The significant leading roles are in bold. We can observe the larger delay patterns in Fig. 6b, d, e, g, h, especially in Fig. 6e where q τ (n) > 50. Based on the analysis of event delay patterns of the whole network, as shown in Fig. 7, we can particularly identify the intra-group drive-response links, such as the F08 → M09 → C09 in group I, D04 → D03 in group II, J 07 → G05 and J 07 → H 10 in group III, and K 01 → K 02 in group IV, which may be difficult to detect in clinic but facilitates to determine the leading role node in each group. However, the inter-group drive-response relationships or leading roles conforming to the clinical observation has not clearly captured. In fact, inter-group interactions is progressed in a scattered and staggered manner. For example, G05 (see Fig. 5a) or J07 (Fig. 6g) have the significant leading roles (event delay) for about 5 or 6 other channels in other groups, respectively. They thus may be the potential common hidden driving source. However, observations in detail on the timing reaching to the threshold q τ (n) = 10 show that there exit evident relative biases in the five event delay patterns. As  Fig. 7b, c, this may suggest the presence of different spatiotemporal propagating delays from G05 or J07 to these five channels, respectively. In contrast, the other channels only display the sporadic leading roles. We should remark that Fig. 7b, c clearly show the presence of a time delay of one signal with respect to the other, but do not necessarily prove a spatial and temporal disparities of signal propagations, although it might suggest it. These results may assist us to design the effective network control strategy.

Seizure focus of directed network
In this section, we aim to reconstruct the directed seizure network based on the above analysis of the time-resolved event delay patterns. Then, based on this constructed directed network we try to seek the potential seizure focus nodes using the in-degree and outdegree distributions of networks. The channels with the maximum out-degree are regarded as the theoretical epileptogenic focus. Figures 8b, c give the reconstructed directed network before (Fig. 8b) and after (Fig. 8c) seizures. Clinical stereo-EEG report shows that in the inter-ictal period and ictal period F08/D03 and D03 are the possible epileptogenic focuses, corresponding to the node 1 and node 4, respectively. However, our results, shown in left panels of Fig. 8b, c, suggest that node 8, i.e., J07, with largest out-degree and no in-degree before seizures, may be the primary focus. Detailed observation shows that before seizure node 8 has the leading role for node 1 and ignites it. After seizure, the leading role of node 8 to node 1 turns to node 6 (i.e., G05) and the leading roles of node 6 to others (i.e., nodes 2, 3, 4, 5, 9 and 10) are doubled than that before seizure. In sum, nodes 6 and 8 traverse almost all of the network nodes. Therefore, they potentially play the vital role to the epileptic network seizures. Due to the complexity of epileptic network, these results may broaden the understanding of seizure origins from the network view. Our another desire is to identify the key nodes in directed network by utilizing the theory of complex network that can effectively control the network seizures and can also sense the full internal state of complex system. That is, we will explore the controllability and observability of directed network, respectively.

Controllability
Note that, due to the similarity of exploration for the network before and after seizures, we only focus on the post-seizure case. The post-seizure intra-network connectivity can be denoted by the following adjacent matrix, where a i, j , i, j = 1, 2, ..., 10, calculated by Eq. (14), indicates the leading role of channel j relative to i. As shown in the schematic diagram of Fig. 8a, we wish to know whether it would be possible to apply a certain number of control signals to typical nodes which can drive the network system toward some desirable state, i.e., the network controllability. Of course, the whole network would be controlled with high likelihood if each node receives one control signal. But this will cost high. Thus, in the premise of bringing system under control, the minimum number of nodes, N d , prerequisite for controlling seizures are thus a natural aim of interest. The network states before (b) and after (c) seizures can be controlled (rightward arrow, controllability) to reach to some desired final state from initial state by inputs stimulation perturbations, u = (u 1 (t), u 2 (t), u 3 (t), u 4 (t)), applied on unmatched nodes (e.g., nodes 1, 4, 6 and 8). In turn, the initial state can also be inferred from the current state (leftward arrow, observability, see also Fig. 9). The networks in the left of (b) and (c) also show the heterogeneous degree distribution. The node size denotes the out-degree of the network. The yellow nodes have only the indegree. The width of the arrows represent the significance of leading role (q τ (n) >10). The right of (b) and (c) show all the four matching paths (green arrows) orderly connected by pink arrows. Other links are marked by light gray Liu et al. [37] proposed the maximum-matching set of network nodes (i.e., the nodes with incoming links but not sharing the starting or ending links) which can be used to assess structural controllability of the built network. In other word, we can also identify the minimum number of unmatched nodes (viewed as driver nodes), N d , (the number of maximummatching set is N M = N − N d ) using the structural controllability theorem and the minimum inputs theorem. That is, N d nodes are required to be con-trolled for the system to satisfy the full-rank condition (18).
In our constructed network, according to the minimum inputs theorem, N d can be calculated as N d = 4. And as an example the nodes 1, 4, 6 and 8 are identified as a driving nodal (unmatched nodal) set. Thus, in Eq. (18) M = N d = 4 and the input matrix B is given as follow: It is easy to be calculated that rank(J 1 ) = N = 10 is full rank and the network is structurally controllable. Note that, the driving nodes are not unique, but their minimum number is. This suggests that the controllability of epileptic network depends not only on the clinical or computational seizure focus but also on the method of network analysis. In particular, the control channels may be not unique. In fact, following the above structural controllability theorem, the unmatched nodal sets {1, 5, 7, 8}, {1, 2, 6, 8}, {1, 2, 7, 8} can also satisfy the full-rank condition of controllability. Similar to [38], corresponding to matching nodes, we can determine the maximum-matching links (i.e., the green arrows in Fig. 8b, c that do not share the starting or ending nodes. The number of maximummatching links is equal to the number of maximummatching nodes, i.e., N M . We can next give the matching paths that is the subset of the maximum-matching links. As shown in the right panels of Fig. 8b, c, before and after seizures the network can be divided into four independent matching paths, respectively. Each matching path is either (i) an "isolated" node (node 8) or (ii) starts from an unmatched node and ending at a matched node without outgoing link belonging to the maximum-matching set. Before seizure, these matching paths shown in Fig. 8b are 1→7, 8, 6→5→2→3, 4→10→9, respectively. After seizure, these matching pathes shown in Fig. 8c are 1→5→2→3, 4→10, 6→9, 8→7, respectively.
Another issue of significant practical interest in network control is to reduce the N d so that the network can be harnessed by less control signals. In particular, when N d =1, only one node is unmatched and the network can be optimally controlled by mere one input control signal, i.e., which can make the network meet the full-rank condition [37,38]. However, the matching paths determined above are abjoint. A way to reduce N d is to make unmatched nodes to be matched. This can be implemented by enhancing beforehand the controllability of the network via adding links between two independent matching paths because the topological changes can alter the network controllability according to [37]. In addition, to make the network fully controllable under only one external controller, i.e., N d =1, we need to identify a minimum number N − N M − 1 of additional links. In particular, we first randomly order the found four matching paths and then link the ending nodes of each matching path to the starting nodes of the matching paths next to it in order, as illustrated in right panels of Fig. 8b, c. Then, the network can be fully and structurally controlled with a single controller applied on the starting node of the first matching path (here, the driver can be the node 1), i.e., satisfying the full-rank condition with a single input (see Supplemental Material [38]). Note that, the obtained network in our paper doesn't include close matching paths where any one node can be viewed as the driver. These results may be instructive to the epileptic network reconstruction and modification.

Observability
A complex system can be quantitatively described and observable by estimating and reconstructing the system's internal state, respectively. However, this description is constrained by our estimating ability from experimentally accessible outputs. This is mainly because in practice experimental access is limited to only a subset of variables, or sensors, even though the simultaneous measurement of all internal variables can offer a complete description of a system's state.
Here, based on the reconstructed network in Fig. 8, we try to search for the potential sensor nodes that are necessary to estimate the full internal state of the epileptic network system. Because the epileptic network can capture the pathological information flow in inferring the state of individual variables, we give the inference diagram in Fig. 9. Note that, the inferred directed diagram for observability analysis has the reverse direction ( Fig. 9) in contrast to the network for controllability analysis (Fig. 8). We suppose the system is estimated by the following equations: Fig. 9 (Color online) Graphical approach of structural observability of epileptic network before (a) and after (b) seizure, corresponding to Fig. 8 but having the opposite directions. The inferred balance equations of seizure network is shown in Eq. (24). Inference diagram is constructed by drawing a directed link (x i → x j ) if x j appears in the right-hand side of x i 's balance equation. The red nodes having no incoming links is termed as root sensor nodes (observation group) whose measurements allow us to reconstruct the state of all other variables or to ensure observability of the whole system where the variables before the semicolon in each equation correspond to the pre-ictal case while the variables after the semicolon correspond to the post-ictal case. The inferred directed link x i → x j in the inference diagram is determined if x j appears in x i 's differential equation. This means one can get the information on x j by observe x i . (•) in x 8 's differential equation means no information on any other variables can be captured by monitoring x 8 . Same to the (•) in other differential equations.
In the Fig. 9a, all the nodes except for node 6 are sensors while in the Fig. 9b all the nodes except for nodes 1 and 8 are sensors. However, it is particularly interesting to identify the minimum set of sensors from whose observability one can infer the states of any other nodes. In particular, we define the sensor nodes having no incoming edges as the root sensor nodes shown in  (26) In addition, it is found that even though these root sensor nodes are not identical before and after seizures, they belong to the same group I, III and IV, respectively, as shown in Fig. 7. The information on nodes 4 and 5 in group II can be captured by monitoring the nodes in other groups. Here, node 4 is the clinically identified focus. These results suggest that in terms of the epileptic network observability the root sensor nodes may be not always the clinically localized focus node.

Reproduction of the clinical results using a spatially extended neural network model with the stimulation modulation of cNDNS
Computational model is a promising method to investigate the epileptic brain network dynamics [52,53]. Finally, we try to reproduce the clinical results using the single-compartment neural field model of Eq. (20).
To computationally reconstruct the clinical epileptic network, we need to extend the above model to model network (Eq. 21) by coupling different populations, denoting the various cortical columns or brain structures. During the simulation, the a i j (t) depending on time in Eq. (14) obtained by the EEG analysis is used to describe the coupling connectivity of network with time going. Presently, we take n = 10 to reconstruct the epileptic network model representing the 10 spatial locations corresponding to the clinical 10 contacts' EEG recordings. That is, during the 10-compartment epileptic network, population i (∈ {1, ..., 10}) receive the outputs of EPSP transfer function h d of the other 9 populations, i.e., j∈{1,...,10}; j =i K i j x j 4 . Throughout a seizure, the directed effective networks can represent the distinct regimes of neural connectivity. Figure 10 gives the effective connectivity of directed network during pre-ictal (Fig. 10a), trans-ictal (Fig. 10b) and post-ictal (Fig. 10c) across ten nodes. It is interestingly shown that before and after seizures the system network shows cluster connections, while during seizure transitions the network system is scatteredly connected. This suggests that the seizure transitions from pre-ictal to post-ictal periods is along with the evolution and recombination of information flows.
The simulated results can be found in Fig. 11. Compared Fig. 11a with Fig. 4, we can see the qualitative agreement between original EEG signals and the simulations. Due to the advantage of painless and unharmful neuromodulation of cNDNS, i.e., normally distributed noise stimulation with cathodic phase pulses (here supposing the brain has the noise with anodic phase pulses), we then test the effect of cNDNS on the (a) (b) (c) (d) Fig. 11 Time series for the simulated ten channels corresponding to Fig. 4a. Simulated EEG a without stimulation perturbation; b with stimulation perturbations applied on channels F08, D03, G05, J07; c with stimulation perturbations applied on channels C09, K01; d by resecting channels F08, D03, G05, J07. The integration step size of 4ms (sampling frequency: 256 Hz) is used. A in Eq. (21) is equal to linspace (4.35 mV, 4.75 mV) for first two channels and 3 mV for other channels seizure abatement. In particular, we consider the effect of cNDNS applied on the typically selected contacts, e.g., the nodes 1, 4, 6 and 8 (i.e., F08, D03, G05 and J07), respectively, which are the driving nodes identified from the controllability analysis at last section. As shown in Fig. 11b, the epileptic discharges are almost completely abated. Whereas, cNDNS is feeble in abating seizures when applied on the non-driving nodes, e.g., nodes 3 and 9, i.e., C09 and K01, as shown in Fig. 11c. In addition, we also get rid of the nodes F08, D03, G05 and J07 from the network reconstruction to simulate the clinical resection of seizure focus, which shows the similar seizures abatement effect with cNDNS. This evidences the effectiveness of graph theory in identifying control nodes of seizures. These computational results also verify the dynamical network characteristics of epileptic seizures.

Method reliability
Note that, the focal cortical epilepsy is of very heterogeneous nature and in clinical practice sampling of the brain is usually very under sampled spatially and differ greatly patient by patient. Thus, it will be very important to test the method by application on mul-tiple dataset consisting of several patients intracranial EEG data and present the results of the testing to prove the usefulness of the method in this field and its clinical applicability. For this purpose, in Table 1, we provide the results about SEEG data of another 9 focal epilepsy patients with more than 100 seizure contacts. Importantly, the results similar to Fig. 11 are obtained. In fact, it is noted that in clinic each patient has more than 10 electrode channels with more than 100 contacts. Nevertheless, as shown in Table 1, we randomly selected about 10 contacts from each patient to construct the effective network. It is obviously shown that the calculated focal locations are largely distinct from the results of SEEG report. However, the shared focal locations can be observed between them which may be the network nodes needed to be paid more attentions in practice. In addition, although the calculation shows that for most patients (60%) to full control seizures needs performing inputs on more nodes, it is interestingly observed that for the patients 3, 6 and 9(30%), less control nodes are needed to modulate network seizures. This results may provide new insights into the clinical treatment decisions

Conclusion and discussion
This paper describes an application of a graph theory and system identification for finding an epileptogenic nodes of epileptic seizure substrate to inform the epilepsy surgery or targeting the electric stimulation of the brain to suppress epilepsy seizures. Firstly, the reliability of a previously proposed methodology employed to estimate the synchronicity and directionality of information flows over time between EEG signals, is numerically assessed with a coupled mass neural model. Then 10 channels' EEG signals from a patient with focal epilepsy are used to reconstruct the dynamical complex network of pathological seizure. Based on the controllability and observability principles of complex systems, we can focus on the key nodes which is effective to control the network seizure behaviors and the key ones that can allow us to estimate the state of all other variables. Results show that the full control of epileptic network may not only related to the focus zone, and the other non-focus nodes could also play important roles. In addition, we use the spatiotemporal neural network model connected by our modeled dynamical adjacent matrix to successfully reproduce the original EEG signals which can be effectively abated by applying the normal distribution noise stimulation with cathodic phase pulses (cNDNs) on the identified key nodes or resecting them. More importantly, to validate the reliability of the proposed method, we provide the results about SEEG data of another 9 focal epilepsy patients with more than 100 seizure contacts.
In conclusion, the results obtained above evidence the effectiveness of graph theory in identifying control nodes of seizures and also verify the dynamical network characteristics of epileptic seizures. In particular, we computationally found that applying stimulations on the clinically localized seizure focus just partially abates seizures. That is not only the focus zone of clinic report but also the other non-focus nodes may contribute to the full control of the epileptic network. This suggests that some focus zones of clinic report may be not the real focus, or at least should not be the clinic surgical positions of abating seizures. This also hints that the clinically localized focus may be driven by a hidden source, which facilitates to the search for the potential recessive, dominant or extended foci. These results largely enrich the clinical results and provide new insights into the seizure resection and electronic stimulation therapies.
It is noted that each patient has more than 10 electrode channels and each channel has about 10 contacts located on the different depth of the electrode. Thus, in future, the multi-layer and deep-layer network across more nodes will be considered with the aim of informing more exact clinical decisions. We also should note that the key nodes obtained by using our proposed strategy are not unique. In particular, during the simulation, the parameters used to reconstruct the seizure networks of the different patients are varied and patient-specific. In view of this, in the following work the uniform computational framework solving these issue will be further explored. In addition, in Eq. (20) or (21) p(t) globally measures the influence of neighboring or distant populations. This influence was assumed to be Gaussian white noise which had been evidenced to play roles in the seizure generations. The principle using the normal distribution noise stimulation with cathodic phase pulses (cNDNs) to abate seizures is just to counteract the effect between two different populations. At last, it has to mention that the practical confirmation of the theoretical seizure abating capabilities with the methods presented in this manuscript has not been performed. This is what we will focus on in the future.

Final remark
Seizure focus localization is the key to control seizures. But how to accurately identify the focus which can effectively control seizures is still an open challenge. The controllability of complex system using graph models is a great vehicle and shown to be an effective tool in exploring of biologic complex systems as well as in controlling of such. Even in epilepsy, there has been multiple approaches using a graph theory to better describe epileptogenesis and target the therapies (either surgery or electrical brain stimulation). However, there still lacks a canonical computational framework or research strategy for seizures focus localizations that could comprehensively consider the time-evolving information flow and the analysis of controllability mechanisms driven by the real seizures data, as well as the the computational validation. Based on latest research in epilepsy and graph theory applications, our work proposes a general methodology or strategy for seizures focus localizations that can indicate the key players (network nodes) in generating seizures, which may also provide a good way to go and explore to inform clinical decisions.