Evolving Network Model that Almost Regenerates Epileptic Data

Increasingly, many complex systems are being modeled as networks since the framework of nodes representing the basic elements of the system and the interconnections of the network representing the interaction between the elements fits well for a theoretical study. When the complex systems are large-dimensional dynamical systems, the network framework comprises many interacting subsystems of smaller dimension each of which constitutes a node. As a particular example, the whole or a portion of the entire neuronal activity in the human brain can be regarded as the consequential dynamics of interacting subsystems, where the dynamics of a subsystem is generated by a group of neurons. The enormously interconnected subsystems in the brain generate a wide variety of dynamical patterns, synchronised activities and rhythms.

Epilepsy is a disorder that affects the nerve cell activity which in turn intermittently causes seizures. During such seizures, patients could experience abnormal sensations including loss of consciousness. Clinical and theoretical research have shown that underpinning the cause of this disorder has not proved to be easy, in particular the predicament of resolving the question as to whether the disorder manifests due to the nature of the interconnections in the brain or the pathologicity of a portion of the brain tissue itself remains.

Since epilepsy is one of the most common neurological disorders with an estimate of more than 50 million individuals being affected [26], there is a strong need both for curative treatment and as well for the analysis of the underlying structure and dynamics that could bring about seizures. In this paper, our attention is on focal epilepsy where the origin of the seizures are circumscribed to certain regions of the brain called the seizure focus [4, 3, 7, 19, 27, 32, 36], and the aim of this paper is to foster theoretical investigation into the connection strengths of the underlying nodes in such a seizure focus in comparison to the other nodes in the network.

Different models have been proposed to understand different aspects of focal epilepsy [1, 13, 38]. Mathematically speaking, dynamical system models of focal epilepsy studied in the literature are mainly of two broad categories: (i). Models that consider noise to be inducing seizures in a node [15, 31, 34] (ii). Models that consider seizures to be caused or terminated by a bifurcation, a change in the intrinsic feature of the (autonomous) dynamics at the node [5, 24, 28].

While modeling the neuronal activity of the human brain, there are far too many numerous parameters which also dynamically evolve on separate spaces and scales, and it is unrealistic to observe or encapsulate all these aspects in a static network or an autonomous dynamical system (see Appendix A). Since the brain evolution does not depend on its own internal states as it responds to external stimuli and the interconnections in the brain are likely to change in response to stimuli, the human brain can be safely said to be a nonautonomous system and in the language of network dynamics, this is effectively an evolving network (see Section 2).

Besides evidence that network connections play a role in initiating seizures (e.g., [2, 4, 20, 32]), some authors also provide evidence that stronger interconnections play a role in epileptic seizures (e.g., [11, 25]). Also, a commonly found feature of biological systems is adaptation: an increase in exogenous disturbance beyond a threshold activates a change in the physiological or biochemical state, succeeded by an adaptation of the system that facilitates the gradual relaxation of these states toward a basal, pre-disturbance level. Based on all these, we present a novel phenomenological model where seizures could be activated in a node, if certain interconnections change to have larger weights resulting in a larger exogenous drive or disturbance or energy transmitted from other nodes into it, and upon experiencing seizure, adjustment of interconnections takes place eventually to subside the seizures. Since the interconnections are allowed to change, the model dynamics emerges from an evolving network.

In this paper, we propose a patient and data-specific, macroscopic, functional network model that can act as a “surrogate” representation of the electroencephalography (EEG) or electrocorticography (ECoG) time-series. The network model is new in the sense that it satisfies the following simultaneously: the network connections are directed and weighted; the network connections are time-varying; there is a discrete-time dynamical system at each node whose dynamics is constantly influenced by the rest of the network’s structure and dynamics; the ability of the network to generate a time-series that sharply resembles the original data is based on a notion of robustness of the dynamics of the network to its initial conditions what we call as the “state-forgetting” property.

The essence of the state-forgetting property, roughly put is that the current state of the brain does not have an “ distinct imprint” of the brain’s past right since from its “genesis” or “origin” (rigorous definitions are provided later in the paper). Effectively, this means that the brain tends to forget its past history as time progresses. Such a notion has been adapted from the previously proposed notions for studying static (artificial neural) networks with an exogenous input in the field of machine learning viz., echo state property [12, 22], liquid state machines, [21], backpropagation-decorrelation learning,[33] and others. However, the network model that we represent is very different from that of the artificial neural networks and moreover the nodes in our model correspond to different regions of the brain.

Our methodology is as follows: We present an evolving network model which exhibits state-forgetting property for all possible inter-connections, and a particular sequence of interconnection strengths (matrices) is determined from a given EEG/ECoG time-series so that the model can almost regenerate the time-series. The weights of the interconnections are chosen from a set of possibilities based on a heuristic algorithm that infers a causal relationship between nodes – a (directed) synchrony measure between two nodes $i$ and $j$ represents the conductance (ease at which energy flows) from node $j$ into node $i$, and weight of the connection from node $j$ into node $i$ has some degree of proportionality with the synchrony measure.

The state forgetting property ensures the synthesized evolving network has in some sense captured the original time series as an “attractor” (see [23] and Proposition C.1), and hence the network acts as surrogate representation of the original time series. In other words, since any initial condition presented to a state-forgetting network has no effect on its asymptotic dynamics, and that the asymptotic dynamics resembles the time-series, we assume that the evolving network model has captured the essence of the time-series data.

To indicate the potential utilities of the model, we analyse the constructed evolving network from a patient’s clinical data. By comparing the relative strength of the interconnections, we identify nodes that are outgoing hubs of the network as a theoretically determined seizure focus. The nodes in such a focus have a strong location-wise correlation with the set of nodes in a clinically determined seizure focus. We also, discuss the efficacy of surgical removal of the determined focus in subsiding seizures.

The remainder of the paper is organised as follows. We introduce the notion of state-forgetting evolving networks in Section 2; formal definitions are made in Appendix B. In Section 3, we present the evolving network model in detail and state a result concerning the state-forgetting property; the complete mathematical proof of this result is presented in Appendix C. In Section 4, we present a method to obtain the weight matrices to obtain a time-varying network from a clinical time series. In Section 5, we present some results on the network inference followed by a brief discussion and conclusion in Section 6.

The aim of this section is to give an intuitive description of the state forgetting property of evolving networks; the formal and rigorous definitions of this property are found in Appendix B. We employ discrete-time systems [6] since they are more amenable to computer simulations as the computed orbits or solutions unlike in the case of differential equations do not depend on numerical methods for computing the solutions and their related accuracy. Also, discrete-time systems can exhibit complex dynamical behaviour (e.g., [6]) even in one-dimension whereas a continuous time-system needs a dimension of at least three.

Using Fig. 1, we present a less formal and simplistic description of an evolving network with three nodes $i=1,2,3$. Let $W(n)$ be a $3\times 3$ matrix defined for each integer $n$ with $W_{ij}(n)$ being the strength of the incoming connection from node $j$ to node $i$; here $W_{ij}(n)$ is the element in the $i^{\mbox{th}}$ row and $j^{\mbox{th}}$ column of $W(n)$. There is a discrete-time dynamical system at each node in the evolving network represented in Fig. 1, and the $i^{\mbox{th}}$ node dynamics is given by

where $h_{i}$ is some function (not dependent on $n$) that calculates $x_{i}(n+1)$ using information from the dynamical state in the three nodes at time $n$ and the respective incoming strengths into node $i$. By specifying $h_{1}$, $h_{2}$ and $h_{3}$ together with the sequence of weight matrices $W(n)$, the evolving network is specified completely. By denoting $x(n):=(x_{1}(n),x_{2}(n),x_{3}(n))$, we can represent the evolving network dynamics also as a nonautonomous system $x(n+1)=g_{n}(x(n))$, where the variation of $g_{n}$ with $n$ represents the inherent plasticity of the system (change in $W(n)$ above). We can also allow $g_{n}$’s to represent external stimuli to the network and other obfuscated factors.

Since an evolving network is essentially a nonautonomous system, we recall some relevant results of a nonautonomous system. Of interest to us is the notion of an entire-solution of a nonautonomous system. An entire solution (or orbit) of a nonautonomous system $\{g_{n}\}$ is a sequence $\{\vartheta_{n}\}$ that satisfies $\vartheta_{n+1}=g_{n}(\vartheta_{n})$ for all $n\in\mathbb{Z}$. It is possible to show under natural conditions that there exist at least one entire solution [22]. In the particular case, where the evolving network has a unique entire solution, we say that the network (or the nonautonomous system) has the state-forgetting property (formally defined in Definition B.2).

Now, we comment on the rather intriguing consequence of an evolving network $\{g_{n}\}$ having multiple entire-solutions. Suppose $\{\vartheta_{n}\}$ and $\{\theta_{n}\}$ are two distinct entire-solutions of a nonautonomous dynamical system, then for some time $n=n_{0}$, we have $\vartheta_{n_{0}}\not=\theta_{n_{0}}$ i.e, the system could have evolved into two possible states at time $n_{0}$. Since $\vartheta_{n_{0}}=g_{n_{0}-1}(\vartheta_{n_{0}-1})$ and $\theta_{n_{0}}=g_{n_{0}-1}(\theta_{n_{0}-1})$, it is necessary to have $\vartheta_{n_{0}-1}\not=\theta_{n_{0}-1}$. Similarly, $\vartheta_{n_{0}-1}\not=\theta_{n_{0}-1}$ would necessitate $\vartheta_{n_{0}-2}\not=\theta_{n_{0}-2}$ and so on. Thus, in general, at any time-instant $n_{0}$, $\vartheta_{n_{0}}\not=\theta_{n_{0}}$ implies that there are two entire-solutions that have had distinct infinitely long pasts. In other words, the dynamics of such systems with distinct entire-solutions “depends” not only on the inherent plasticity or the exogenous stimuli represented by the maps $\{g_{n}\}$, the distinct solutions carry a distinct and non-vanishing imprint of their past.

From the arguments in the preceding paragraph, for nonautonomous systems like the human brain, having more than one entire solution translates to saying that the current state of the brain is not just the result of the changing interconnections in the brain or the exogenous input but to also say that its past states have a non-vanishing effect on the asymptotic dynamics.

Since it is hard to envisage or believe that seizures in current time had something to do with the state of the brain right from its “genesis” (genesis is used as a metaphor to refer to the earliest factual time instant from which the current brain state has evolved from, and this time could be matter of debate), we assume the more plausible – seizures at a certain time are a consequence of how the network structure has evolved/changed till that time or a result of a stream of exogenous inputs to the brain, but certainly not a consequence of the state at which the brain was in during its genesis. This also forms the basis for employing echo-state networks in artificial neural network training and machine learning[12, 22].

Whenever a nonautonomous system (technically, defined on a compact state space) has the state-forgetting property, it also has state-space contraction towards its entire-solution in the backward direction (see Proposition C.1). In fact, a stronger result [23] holds: if $\{x_{n}\}$ is the unique entire solution, then for every $y_{0}\in X$ the successive iterates $y_{n+1}=g_{n}(x_{n})$ would satisfy $|y_{n}-x_{n}|\to 0$ as $n\to\infty$. In other words, $\{x_{n}\}$ “attracts” all other state-space iterates towards it (component-wise) as time tends to infinity. This convergence due to the attractor property of the state-forgetting property is what makes our evolving network model generate a time-series that actually sharply resembles (mathematically, converges to) the data from which the network was conceived. Further discussion on attractivity is beyond the scope of this paper. The larger point is that for an evolving network to be of utility as a surrogate representation of the time-series data, we need the state-forgetting property.

Before, we actually describe the evolving network model in (23), we list and describe its features: (i). The evolving network model has $N$-nodes for all time. (ii). Placed at each node is a certain bi-parametric one-dimensional discrete-time dynamical system which is autonomous when there are no incoming connections from other nodes to it, and nonautonomous otherwise. We refer to the dynamics of the system at a node, as the dynamics of the node or nodal dynamics. (iii). Among the two parameters that specify this one-dimensional dynamical system is a stability parameter. The nodal dynamics in autonomous mode has a globally stable fixed point at $0$ whenever the stability parameter is in a certain range (stable-regime); a globally stable fixed point attracts all solutions, i.e, the time series values measured in the node are close to $0$ and tend towards to it. (iv). When the stability parameter is chosen not to lie in the stable regime, the autonomous nodal dynamics can render a range of different behaviours including oscillatory behavior, complex behavior like chaos [6] depending upon the value of the stability parameter. (v). When the nodal dynamics is nonautonomous, the dynamics of the node is perturbed at each time-instant from the dynamical states of other nodes. We call the net perturbation from other nodes as the drive to the node under consideration. The exogenous drive from other nodes has its signature on the nodal dynamics, and if sufficiently strong, is also capable of setting the nodal dynamics far away from the fixed point $0$. (vi). We term a seizure in the node to be a dynamical excursion away from the fixed point, and that the seizures in a node are brought about (and terminated) in two distinct scenarios: (via). changing the stability parameter of the node in the autonomous mode of operation (vib). by the nature of the exogenous drive from the other nodes in the network in the nonautonomous mode. (vii). In this paper, we consider the case of nodal dynamics slipping into seizures as in scenario (vib) above, i.e, we consider each node to have its stability parameter to lie in the stable regime, and the exogenous drive to cause seizures. The drive is time-varying on two accounts: first, being a function of the states of the other nodes it is time-varying as other nodes have a dynamical system and not static; second, the drive is time-varying since it is a function of the interconnection strengths of the nodes.

Our model is also motivated on the lines of previous research [14, 30, 35, 37] showing that the hyper-interconnected regions in the brain are prone to seizures (see Section 3). Hence, we expect the drive to an individual node is likely to be having large magnitude if it has large weighted incoming connections.

To describe a family of evolving networks (as in Definition B.1), we first sketchily describe the nonautonomous dynamics at each node $i$ (the map $h_{i}$ described in Section 2 or in Definition B.1) via the update equation:

where the various entities are listed and briefly explained:

• •

  The iteration of the map $f_{a,b}$ on an initial condition gives the autonomous dynamics of the node. Some Neuron models like that of Fitzhugh-Nagumo (e.g., [9, 29]) permit an individual neuron to have a regenerative self-excitation via a feedback term. Based on this principle, we consider a highly simplistic discrete time model of the node dynamics where feedback from the node dynamics can cause self-excitation, and set

  $$f_{a,b}(x):=a\>x^{3}-b\>x,$$

  (19)

  where parameters $0<a<1$ and $0<b<3$ (for a graph of $f_{a,b}$, see Fig. 2(a)). Explanation on how a parameter can cause self-excitation dynamics is made in Section 3.1.

• •

  Adopting the principle that the individual drive from a connecting node has a larger magnitude if the incoming connection weight and the magnitude of the state value of the connecting node are together larger, we set the total sum of the drives from the other $(N-1)$ nodes into the node $i$ to be $=$

  $$-\>\sum_{j\not=i}\log\Big{(}\phi(W_{ij}(n))\cdot(d+x_{j}(n))\Big{)},$$

  (20)

  where $W_{ij}(n)$ is the weight of the incoming connection from the $j^{\mbox{th}}$ node to the $i^{\mbox{th}}$ node; $d>0$ is an offset parameter to enable a positive argument for the $\log$ function; $\phi:[0,\infty)\to[0,\infty)$ is an increasing function and is defined by

  $$\phi(\star):=\frac{1}{d}\;\;e^{\sqrt{\log(1+\star)}}.$$

  (21)

  The function $\phi$ although increasing, has a monotonically decreasing slope (see Fig. 2(c)). The net effect of the function $\phi$ can be thought of as the node’s mechanism to have a relatively stronger inhibition to a stronger incoming connection – node’s internal mechanism to resist seizures by inhibiting stronger exogenous drive. Lastly, the negative sign in (20) is only for mathematical convenience and has no physical significance. When the drive is zero, the dynamics of the node $i$ is in autonomous mode.

• •

  The saturation function henceforth would be the mapping $\sigma:\mathbb{R}\to\mathbb{R}$ which is linear in a certain adjustable range, and its slope becomes zero asymptotically (see Fig. 2(b)) defined by:

  $$\sigma(y)=\begin{cases}y&\text{: $y\in[-p,p]$}\\ \frac{1}{r}\log(\frac{1}{r}+(y-p))+k&\text{: $y>p$}\\ -\frac{1}{r}\log(\frac{1}{r}-(y+p))-k&\text{: $y<p$},\end{cases}$$

  (22)

  where $p>0$ and $\sigma$ is the identify function on $[-p,p]$; $r>1$ determines the asymptotic rate at which the slope of $\sigma$ decreases and $k=p-\frac{1}{r}\log(\frac{1}{r})$. (see Fig. 2(b)); $\sigma$ is linear in the region $(-p,p)$ with a slope $1$, and the slope of $\sigma(y)$ approaches $0$ as $y\to\infty$ monotonically according to $\frac{1}{1+r(y-p)}$.

Putting together the above into (18), the update equation at the $i^{\mbox{th}}$ node is modeled by:

Note that in the model in (23), the weights $W_{ij}(n)$ are unspecified, but they actually can be solved for given a ECoG/EEG time-series $\{x_{n}\}$ (topic of Section 4). In Section 3.1 and 3.2, we describe the dynamics of a node when it operates in autonomous node and in the nonautonomous mode. In Section 3.3, we state sufficient conditions on the parameters $a$,$b$ and $d$, to ensure state-forgetting property of any evolving network arising out of (23). However, to understand the model, a few remarks are in order concerning the role of the offset parameter $d$ and the saturation function $\sigma$ in (23):

From (20), it follows that the individual drive from the node $j$ into node $i$, at time $n$ is given by $u_{ij}(n):=-\log\Big{(}\phi(W_{ij}(n))\cdot(d+x_{j}(n))\Big{)}$. To make the role of the weights meaningful, a smaller weight $W_{ij}(n)$ should ensure smaller individual drive from node $j$, i.e., we must have $u_{ij}(n)$ close to zero whenever $W_{ij}(n)$ is close to zero. We note that this would be the case whenever the offset parameter $d$ is sufficiently large. To observe this, we remark on the behavior of $u_{ij}(n)$ as the weight $W_{ij}(n)$ gets close to zero. By definition of $\phi$ in (21), when $W_{ij}(n)$ is close to zero, it follows that $\phi(W_{ij}(n))$ is close to $\frac{1}{d}$. When $\phi(W_{ij}(n))$ is close to $\frac{1}{d}$, $u_{ij}(n)$ is close to $\log(\frac{d+x_{j}(n)}{d})$. For practical purposes we restrict $x_{i}(n)\in[-1,1]$ (see Section 4), and hence if $d$ is much larger than $x_{j}(n)$, then $\frac{d+x_{j}(n)}{d}$ is close to $1$, and $\log(\frac{d+x_{j}(n)}{d})$ is close to zero, and thus $u_{ij}(n)$ is very close to zero. However, it is also to be noted that if $d$ is made unreasonably large, the signature of $x_{j}(n)$ on $u_{ij}(n)$ gets weaker, so $d$ is chosen in an appropriate range, although with a large value. Physiologically, the parameter $d$ has no direct relation, but serves as a mathematical plug to overcome the difficulty in modeling the drive directly from the dynamical states of a connecting node that take both positive and negative values. The offset $d$ has another role: its value could determine the state-forgetting property of an evolving network (see Theorem 3.3.1).

The particular choice of the saturation function $\sigma$ is not critical to the model; for the convenience of mathematical proofs, we need $\sigma$ to be differentiable and unbounded and hence its form. On the other hand, for practical reasons, choosing such a $\sigma$, helps in the simulation stage to troubleshoot if each of the parameters $a$,$b$ and $d$ are chosen in an appropriate range lest the network lose the state-forgetting property. While the network does not have the state-forgetting property, the dynamics can be unbounded in the absence of a saturation function since the range of $f_{a,b}$ is unbounded and the numerical values on a computer could indicate an error.

Henceforth, for ease of jargon we call the $i^{\mbox{th}}$ component of a $N$-dimensional time series data $\{x(n)\}$, i.e., $\{x_{i}(n)\}$ the data in the $i^{\mbox{th}}$ channel. Given $\{x(n)\}$, we present a method for solving the weights $W_{ij}(n)$ in(23) based on inferring a causal relationship between the nodes, which we call a synchrony measure.

To explain the methodology of solving weights from (23), we rewrite (23) by rearranging as

where $\sigma^{-1}$ is the inverse of the saturation function $\sigma$. Clearly, given a $N$-dimensional time-series $\{x(n)\}$ of EEG or ECoG clinical data, and having fixed $a$, $b$ in the model, the right hand side of (24) can be computed easily for each $n$ from the $i^{\mbox{th}}$ channel data. Let $r_{i}(n)$ denote the right hand side term in (24).

For a fixed $i$, we have a single equation

where there are $N-1$ unknown quantities $\left(W_{ij}(n)\right)$ as $j$ is varied from $1$ to $N$ with $j\not=i$. Hence, multiple solutions of $W_{ij}(n)$ could exist since we have more unknowns than equations. We regard seizure spread in the network due to sustained energy transmitted from a node that is seizing to the ones which is not. With this principle, we hypothesize $W_{ij}(n)$ to be an increasing parametric function of a certain synchrony measure $\rho_{ij}(n)$ between the two distinct nodes $i$ and $j$. Physiologically, the synchrony measure between two nodes $i$ and $j$ represents the ease at which energy dissipates (a measure of conductance) from node $j$ into node $i$, and hence is directed – nodes $i$ and $j$ could have different interconnections, and energy dissipation between nodes is not symmetric. To get a measure of such conductance, we let the synchrony measure $\rho_{ij}(n)$ depend on (i) a sustained influence from node $j$ into node $i$ measured by a leaky integrator placed between the two nodes (ii) the similarity between the average amplitude or power in the two channels up to time $n$.

The algorithm for computation of directed synchrony $\rho_{ij}(n)$ between two channels is described in Section 4.1. We now explain the idea/reasoning behind some of the steps of the algorithm. The discussion can be read in conjunction with the various steps of the algorithm.

The very first step of the algorithm involves normalisation of the amplitude of data for setting the amplitude of the time-series in any channel lies between $[-1,1]$ (see (S1) of Section 4.1).

In a succeeding stage of the algorithm, we compute the average power at each instant in a window of length $\tau$ in each channel (see definition of $P_{i}(n)$ in (S6) of Section 4.1). It turns out that the time interval in which there is sustained influence of one node on the other by the synchrony measure is determined by $\tau$. We comment on the choice of $\tau$ in Section 5.1

The directed synchrony $\rho_{ij}(n)$ depends on the product of two factors (see (S7) of Section 4.1). One of the factors is the similarity between the powers in the $i^{\mbox{th}}$ and $j^{\mbox{th}}$ channels which we quantify by a similarity index $(1-|P_{i}(n)-P_{j}(n)|)$. Note that due to normalisation in (S1), the average power $P_{i}(n)$ in any window is between $0$ and $1$, and hence $0\leq(1-|P_{i}(n)-P_{j}(n)|)\leq 1$, and the similarity index attains $1$ if and only if $P_{i}(n)$ and $P_{j}(n)$ are identical.

The other factor that contributes to $\rho_{ij}(n)$ is $Q_{ij}(n)$, the output of a (directed) nonlinear leaky integrator placed between two distinct nodes and whose states are evaluated recursively (for details see (S8) to (S12) of Section 4.1)). The purpose of such a leaky integrator is to model the sustained influence of the signal strength in the node $j$ on node $i$ – the word sustained points to the fact that the integrator leaks, and unless the influence is not sustained over a period of time, the influence diminishes. Any such leaky integrator in (S11) of Section 4.1 is a nonautonomous dynamical system in its own right, and hence its future states are obtained by iteration on the current state. In essence, a leaky integrator sums up (or integrates) its previous state with what is fed to it typically as an input, but also leaks a small fraction of the sum while integrating. Here, we set the leaky integrator dynamics for the connection from node $j$ to $i$ as

where the cumulative average power $\bar{\bar{P}}_{j}(n)$ in the $j^{\mbox{th}}$ channel is calculated as in (S9) of the algorithm. This cumulative average fractional power is determined by first calculating the fractional power in the $j^{\mbox{th}}$ channel (see (S8) of the algorithm). The cumulative average of this fractional power is calculated as in (S9) of the algorithm; this cumulative average is calculated iteratively using the elementary idea: to calculate the average of a collection of $M$ numbers $(s_{1},\ldots,s_{M})$, if the average of $M-1$ numbers is given to be $A_{M-1}$, then the average of $M$ numbers $A_{M}$ can be updated by $A_{M}=(M-1)A_{M-1}+s_{M}$.

Once the synchrony $\rho_{ij}(n)$ from the time series from the algorithm in Section 4.1 is computed, then for convenience, we proceed to deduce the intermediate quantity $V_{ij}(n):=\phi(W_{ij}(n))$ from which the actual weight $W_{ij}(n)$ can be readily obtained (see (28)).

We hypothesise $V_{ij}(n)$ to be an increasing parametric function of $\rho_{ij}(n)$ and is in the form

where the parameter $\alpha_{i}(n)$ can be solved as explained below. Substituting $V_{ij}(n)=\frac{1}{d}e^{\alpha_{i}(n)\cdot\rho_{ij}(n)}$ in Eq. (25), we have the following straightforward deduction to get an expression for $\alpha_{i}(n)$:

Once $V_{ij}(n)$ is found, the weight $W_{ij}(n)$ can be easily found since $\phi$ is invertible. From (21), we get

We infer cortical connectivity from a 58 year old male patient diagnosed at Mayo Clinic. The patient had intractable focal seizures (seizures that fail to come under control due to treatment) with secondary generalization (seizures that start in one area and then spread to both sides of the brain). The patient underwent surgical treatment and pre-surgical evaluations were performed using a combination of history, physical exam, neuroimaging, and electrographic recordings. For delineating the epileptogenic tissues, electrocorticography (ECoG) recording was performed at a sampling frequency of 500 Hz for four days continuously with grid, strip, and depth electrodes implanted intracranially on the frontal, temporal and inter-hemispheric cortical regions. A total of nine seizures were captured during the recording period and complete seizure freedom was achieved (ILAE class I) following the resection of pathological tissues in the left medial frontal lobe and medial temporal structures. All recordings were performed conforming to the ethical guidelines and under the protocols monitored by the Institutional Review Board (IRB) according to National Institutes of Health (NIH) guidelines. The clinical data and ECoG recordings for this patient are publicly available on the IEEG portal (https://www.ieeg.org) with dataset ID ‘Study 038’.

In our analysis, we consider only the electrodes in the $(8\times 6)$ grid placed on the frontal lobe where both the seizure onset and its spread were captured. The electrodes which circumscribed the epileptogenic tissues are shaded in black on the schematic shown in Fig. 3(a). A small segment of ECoG recording is sufficient to demonstrate the efficacy of our modelling framework. Therefore, we select an ECoG segment comprising of both pre-ictal and ictal activities to explore how the network evolves during seizures. We remove the noise artefacts and the power-line interference from the raw ECoG signals by applying a bandpass filter between $1-70$Hz and a notch filter at $60$Hz. The ECoG signals from a few exemplary electrodes (labeled with a prefix LG) are shown in Fig. 3(b).

Next, we infer an evolving network synthesized from the ECoG data. We consider the cortical area under each ECoG electrode as a network node and apply algorithm detailed in Section 4.1 to get the directed and weighted edges or interconnections between the nodes so that (23) is satisfied. In particular with parameters $a=1$, $b=0.9$, $d=500$, $\tau=5000$ we apply the algorithm to obtain the weights of the interconnections between the nodes at each time sample of the data. To cross validate that the weights obtained eventually render a state-forgetting evolving network, we adopt the following procedure. With an arbitrary initial condition $\hat{x}(n_{0})$ and the obtained weights at time $n_{0}$, we use (23) to obtain $\hat{x}(n_{0}+1)$. Repeating this, we obtain $\hat{x}(n_{0}+2)$, $\hat{x}(n_{0}+3)$ and so on. The time-series in a few channels obtained from this procedure is shown in Fig. 3(d). Apparently, the time-series in Fig. 3(b) and Fig. 3(d) are similar. Actually, the similarity is true barring a few hundred samples at the beginning (see Fig. 3(c)). This convergence of the model output from an arbitrary initial condition to the given time-series verifies state-forgetting property and validates our modelling framework for the study of network evolution during seizure onset and seizures.

As a passing remark to the reader, the model’s ability to sharply reconstruct the given data as indicated in Fig. 3(c) (plots of the first 2000 samples of an ECoG channel data, corresponding model output and their difference) is robust even for a different choice of parameters in the algorithm. In fact, Theorem 3.3.1 tells us that such a reconstruction is possible as long as the parameters are in a certain range.

With an objective of identifying the nodes crucial for the seizure genesis, we study the network properties before and after the seizure onset. Consequently, we divide the ECoG recordings into two segments: pre-ictal (before seizure onset) and ictal (after seizure onset).

For each segment of this data, we compute an ensemble average of node strengths as explained next. It may be recalled from Section 3 that the drive from node $j$ to $i$ at time $n$ is given by $u_{ij}(n):=-\log\Big{(}\phi(W_{ij}(n))\cdot(d+x_{j}(n))\Big{)}$. For each node $i$, we compute its mean incoming drive strength of node over a time period $T$ as $\displaystyle{\frac{1}{|T|}\sum_{n\in T}\frac{1}{N-1}\sum_{j=1,j\not=i}^{N}|u_{ij}(n)|,}$ where $N$ is the number of nodes, and $|T|$ is the number of sample points in $T$. In a similar vein, we compute the mean outgoing drive strength of node $i$ over a time period $T$ as $\displaystyle{\frac{1}{|T|}\sum_{n\in T}\frac{1}{N-1}\sum_{j=1,j\not=i}^{N}|u_{ji}(n)|.}$

Figure 4 shows the mean incoming and outgoing node strengths for each node during the pre-ictal and ictal periods. Fig. 4 (a) indicates a strong correlation of the mean incoming node strengths with the clinical seizure onset zone or seizure focus (depicted by the electrodes shaded in black) during the pre-ictal period, while Fig. 4 (b) indicates a more uniform spread of mean incoming node strengths across all nodes during the ictal period. The correlation of the mean incoming node strengths with the clinical seizure onset zone during the pre-ictal period can be understood as that the external drives from other nodes play a supportive role in initiating seizures (a hypothesis of our model as well). Since during the ictal period all channels seem to have similar average power due to seizures, all nodes are likely to receive a similar quantum of external drive. Hence, the mean incoming drive strength appears to be uniform.

In contrast, Fig. 4 (c, d) shows that the distinctively higher mean outgoing node strength (during both ictal and pre-ictal period) is much strongly “correlated” with the clinical seizure onset zone than that shown by the corresponding mean incoming node strengths. Therefore, from Fig. 4 (c, d), it is evident that the nodes whose node strengths are shown as columns in red (i.e., nodes ${9,10,11,12,13,18,19,20,27,28,29,35}$) are the outgoing hubs. We refer to this set of nodes as the theoretical seizure onset zone.

Incidentally in this experiment there are 12 nodes in the theoretical seizure onset zone, similar in cardinality to the set of the nodes at the site of clinical resection ( i.e., nodes $(9,10,11,12,17,18,19,20,25,26,27,28)$). We refer to the nodes at the site of clinical resection as the clinical seizure onset zone. It is evident that the theoretical seizure onset zone and the clinical seizure onset zone point to very similar regions in the brain. Thus, the model has the potential utility of identifying an approximate seizure focus zone prior to surgery by analysing some pre-ictal and ictal data.

Although the theoretical and clinical onset zones differ slightly, from a theoretical standpoint of view, we can hypothesise the outcome of the surgical removal of these two zones and compare these outcomes if we have a suitable objective measure. To explain this, let $\mathfrak{R}$ be the indices of the subset of the nodes that are omitted from a network of $N$ nodes in a simulation experiment, i.e, the nodes in $\mathfrak{R}$ are made non-functional by setting $W_{ij}(n)=0$ for all $n$ if either $i$ or $j$ belongs to $\mathfrak{R}$. Using such modified weight matrices in (23) we can obtain a time-series $\{y(n)\}$ of the depleted network can be calculated by starting with any initial condition that satisfies $y_{i}(n_{0})=0$ whenever $i\in\mathfrak{R}$. In other words, the signal has to have zero amplitude in all the channels corresponding to the nodes in $\mathfrak{R}$ since they are assumed to not participate in the network dynamics.

We quantify the efficacy of the removal of the nodes that belong to $\mathfrak{R}$ by calculating the ratio of the average power in those channels outside $\mathfrak{R}$ prior to removal of nodes and after it by:

where $|T|$ denotes the number of samples over a time-interval $T$ and $|\mathfrak{R}|$ denotes the cardinality of the set $\mathfrak{R}$. Similarly the average power in the signal generated from the depleted network can be calculated by $\displaystyle{P_{\mbox{depleted}}=\frac{1}{|T|}\sum_{n\in T}\>\>\left(\frac{1}{N-|\mathfrak{R}|}\sum_{j\in\mathfrak{R}}(y_{j}(n))^{2}\right)},$ and we calculate the efficacy of removing the nodes $\mathfrak{R}$ by the ratio $G(\mathfrak{R})=P_{\mbox{original}}/P_{\mbox{depleted}}.$

We tabulate the values of $G(\mathfrak{R})$ below, where $\mathfrak{R}$ is determined by the following choices: (i) the clinical seizure onset zone (ii) theoretical seizure onset zone (outgoing hubs) (iii) an arbitrary set of 12 nodes (iv) 12 nodes randomly drawn (average of $G(\mathfrak{R})$ is tabulated in this case). Now, we discuss the potential application of the model in a clinical study when surgery is part of the treatment. The aim of a clinical surgery would be to a minimise the removal of tissue from the brain along with subsiding seizures post surgery. In our study above, the removal of a set of nodes $\mathfrak{R}$ from the evolving network is intended to simulate a surgical removal of the tissue beneath the nodes in $\mathfrak{R}$. A time series $\{y(n)\}$ from the depleted network is intended to simulate the resultant clinical ECoG time series observed after surgery, and the quantity $G(\mathfrak{R})$ is meant to measure the efficacy of the surgery. Thus prior to an actual surgery, the theoretical efficacy of surgical removal of nodes can be studied by using the corresponding depleted network. From the tabulated values of $G(\mathfrak{R})$, we observe that the removal of nodes from the outgoing hubs that form the theoretical onset zone maximises the efficacy of a surgery when (the tissue beneath) a set with 12 nodes is assumed to be removed in each simulation. In fact the large value of $G(\mathfrak{R})$ obtained after the removal of the nodes in the theoretical onset zone translates to complete subsiding of seizures in the resultant time-series (not shown here) generated from the depleted network. When a small subset of the outgoing hubs were retained in the network, $G(\mathfrak{R})$ was still observed to be smaller and the seizures continued to persist from the resultant time-series (not shown here) generated from the depleted networks. From this ability of outgoing hubs to subside seizures in the rest of the network upon their omission, we also refer the outgoing hubs to as spreaders. For a reader interested in finding an optimal set of nodes for subsiding seizures, we note that one does not have to necessarily aim to get a very large value of $G(\mathfrak{R})$, but any $G(\mathfrak{R})$ that yields a seizure free time-series is adequate.

We remark that in our simulations, the mean incoming (outgoing) node strengths and the mean incoming (outgoing) weights of a node exhibit similar variation across nodes, but we employ the former as it accounts for more accurate interaction between two nodes since it subsumes both the effect of the weight of the interconnection and the signal strength of the node that connects.

In this paper, we have presented a modeling approach to transform data into an evolving network that has the state forgetting property or robustness to initial conditions. Although, this principle of transformation can be applied to data originating from any source, the particular evolving network model concerned an epileptic network. Although time-series can be modeled with reasonable accuracy by static networks of very large dimension such as recurrent neural networks, the individual nodes in such case have no physical significance and they do not correspond to the regions in the brain. The interconnection-strengths between the nodes were solved via a causal relationship inferred from the sustained power similarity between two nodes.

Since the model can produce a sharp reconstruction of the time-series, there is flexibility to study the interconnections in the evolving network in conjunction with the value of the time series at any time point, or over a time-interval. We have observed from modeling clinical data of an epileptic patient that certain nodes turn into outgoing hubs during the pre-ictal period and persist through the ictal period. These outgoing hubs also have a location in the brain similar to that of the clinically identified seizure onset zone (seizure focus). This hints that the (pathological) seizure causing areas in the brain could be having strong outgoing connections during seizure onset and during seizures.

Further, to quantify the influence of such outgoing hubs on the network dynamics, we have presented a result quantifying the theoretical efficacy of removal of such nodes in subsiding seizures. Of course, the result is based on data of a single patient, but then a comprehensive study of the effect of a removal of any node in an evolving network has to be backed up with a justification – since altering even a single interconnection could alter the entire network dynamics, a reliable model should show small changes in the network dynamics for small changes in the weights. In a sequel paper, we present a mathematical justification for inferring network information from a depleted network where we establish that changing the interconnections results in a continuous change in the time series whenever the evolving network is state-forgetting. Such a result is needed to justify a more detailed study of the theoretical efficacy of removal of nodes. In fact, it is possible to prove a very strong result showing that the effect of the removal of a node can be studied reliably if and only if a ‘nontrivial’ evolving network has state-forgetting property.

The author would like to thank Nishant Sinha at Nanyang Technological University for assistance in presenting plots in Section 5 and some discussions there in. The author thanks the team of the IEEG portal (https://www.ieeg.org) for providing access to the ECoG data. Access to the image(s) and/or data was supported by Award number U24NS06930 from the U.S. National Institute of Neurological Disorders and Stroke. The content of this publication/presentation is solely the responsibility of the author, and does not necessarily represent the official views of the U.S. National Institute for Neurological Disorders and Stroke or the U.S. National Institutes of Health. The author thanks an anonymous referee for pointing out several typographical errors.