Optimal Network Interventions to Control the Spreading of Oscillations

Oscillations are a prominent feature of neuronal activity and are associated with a variety of phenomena in brain tissue, both healthy and unhealthy. Characterizing how oscillations spread through regions of the brain is of particular interest when studying countermeasures to pathological brain synchronizations. This paper models neuronal activity using networks of interconnected excitatory-inhibitory pairs with linear threshold dynamics, and presents strategies to design networks with desired robustness properties. In particular, we develop a dynamical description of the brain through a network where the state of each node models the firing rate of a region of neurons and where edges capture the structural connectivity between the regions. We characterize the presence of oscillations and study conditions on their spreading. We also discuss strategies to optimally design networks which are robust to oscillation spreading. We demonstrate our results with numerical simulations.


I. INTRODUCTION
O SCILLATIONS in the brain have been associated with numerous cognitive processes [1], such as attention [2], memory [3], and the elaboration of sensory information [4], as well as undesirable phenomena, such as tremors [5] and seizures [6] that are linked to several neurological disorders (Parkinson's Disease, schizophrenia and epilepsy). As an example, epilepsy is characterized by seizures, which are a temporary change in the neuronal activity of the brain, marked by excessive or synchronous neuronal activity. This pathological brain activity leads to physical symptoms such as loss of awareness and convulsions [7], which can prove harmful to the patient experiencing them. A recent study [8] revealed the existence of pathological activity in brains of clinically asymptomatic subjects. In other words, epilepticlike oscillations are observed in the brains of both healthy and non-healthy patients, however only in the latter class do localized oscillations spread throughout different regions of the brain into a symptomatic epileptic event. It is hypothesized in [8], [9] that healthy and non-healthy brains differ with respect to their structural robustness to the spreading of localized oscillations. Modification of the network structure to increase structural robustness is a potential avenue for treatment of neurological disorders. This modification can be done via electrical stimulation, as is the case with deep brain stimulation [10], [11], or surgical resection, where connections between certain brain regions are severed in order to prevent the spread of harmful oscillations [12]. This motivates the importance of characterizing the spatiotemporal dynamics of oscillations onset and propagation in the context of pathological activity in the brain [8], [13], [14], as well as developing design principles for modifying the network structure optimally.
In this paper we introduce a mathematical model of the brain based on a network with linear threshold dynamics and use it to characterize how oscillations originating in one region of the brain might, or might not, spread to neighboring regions. The type of oscillations we consider represent seizures and other typically harmful brain activity. We assume that both healthy and unhealthy brains exhibit localized oscillations, but symptomatic pathologies arise as a consequence of the oscillations spreading from one region of the network to another. Our goal is to study how the structure of networks affect their robustness to oscillations spreading, and to develop strategies for the optimal design of robust networks. We validate our results through numerical simulations.

A. RELATED WORK
Modeling brain dynamics through networked systems has a rich history (we refer the interested reader to, e.g., the books [15], [16]). In particular, mean-field models have proven especially powerful when modeling oscillatory trajectories and synchronization in brain activity [17]. Several models have been proposed, e.g., the Jansen-Rit model [18], the Kuramoto model [19], and the Wilson-Cowan model [20]. Although all of the mentioned models have been the subject of numerous studies, we here focus on those involving the Wilson-Cowan model, e.g., see [21]- [24]. We model the dynamics of the networks after a modified Wilson-Cowan model, where the original sigmoidal activation function is replaced by its piecewise-linear approximation. This model, termed Linear Threshold Network (LTN), has received considerable attention [25]- [27]. An unbounded version of LTNs, termed rectified linear units (ReLU), is widely studied in the context of machine learning [28], [29]. In the context of neuroscience, it has been used to model diverse brain states, including goal-driven selective attention [30], [31] and epilepsy [32]. In particular, recent works [30], [31], [33] have studied structural conditions for oscillations in networks of EI pairs by directly linking the lack of stable equilibria to the presence of oscillations. Here, we take a different approach by studying how the dynamical properties of each EI pair in the network changes as a consequence of its coupling to neighboring EI pairs. This allows us to characterize whether a particular node in the network takes part in the oscillation, and to use these conditions to design networks with desired localized patterns of oscillations. Although the literature on network design and topology control of dynamic systems is vast, e.g., see [34]- [36] and the references therein, to our knowledge this is the first work to introduce methods for network design for controlling oscillations in LTNs.

B. PAPER CONTRIBUTION
The main contribution of this paper is to characterize conditions for the spreading of oscillations in brain networks and to formulate and solve optimization problems for the design of networks that are robust to oscillation spreading. In particular, we model the excitatory and inhibitory activity of a small brain tissue (micro-domain) through a Wilson-Cowan model [20] characterized by a piece-wise activation function. We refer to this atomical unit as an Excitatory-Inhibitory (EI) pair, the structural properties of which have been well characterized, see e.g., [30]- [32]. We build networks of EI pairs in order to model the complex interactions among domains of the human brain. Our goal is to exploit the known properties of the single EI pairs to infer global properties of the brain network.
Once formal conditions on the spreading of oscillations are derived, we develop and solve a series of optimization problems through which we can efficiently compute conditions to isolate localized oscillations from the rest of the network. We show how these optimization problems are computationally efficient and practically effective. We conclude the discussion with extensive numerical simulations on synthetically generated networks. Motivated by our previous work [32], in this paper we tie the discussion of oscillations in the brain to the the study of epilepsy; however our results can be also applied to other oscillation-related problems in network control.

II. PRELIMINARIES
This section presents our notation and discusses basic concepts on linear-threshold network models.

A. NOTATION
Let R, R ≥0 , and R ≤0 denote the set of real, nonnegative real, and nonpositive real numbers, respectively. Given a vector x ∈ R n , we use x i to refer to its ith component. Given a matrix A ∈ R n×n , we use A i j to to refer to its i jth component. For

B. PLANAR LINEAR-THRESHOLD NETWORKS
Consider a neural network consisting of N dynamically coupled neuronal populations. The dynamics of the network evolve according to the linear-threshold model where x ∈ [0, m 1 ] × · · · × [0, m N ] is a vector whose ith component is the average firing rate of the ith population, W ∈ R N×N is the weighted adjacency matrix of the network, m is the saturation vector, and u ∈ R N is the vector whose ith component is the average background input to the ith population. We assume that the network satisfies Dale's Law [37], which states that every node in the network is either excitatory (E) or inhibitory (I). If the ith node is excitatory, then the activity of the ith node enhances the firing rate of its neighbors, and the ith column of W is non-negative. On the other hand, if the ith node is inhibitory, then the activity of the ith node suppresses the firing rates of its neighbors, and the ith column of W is non-positive.
To study epileptic events, we are primarily interested in understanding oscillatory trajectories of linear-threshold networks. Because trajectories of (1) remain in the bounded region [0, m 1 ] × · · · × [0, m N ], we say that a solution x(t ) to (1) is oscillatory if x(t ) does not converge to a fixed point as t → ∞. This notion of oscillation includes both periodic and chaotic oscillations, which is a common practice in computational neuroscience [33], [38].
While previous work [33] identifies structural conditions for the existence of oscillations in general LTN models, we focus here on oscillations that arise in LTN networks composed of coupled excitatory-inhibitory (EI) pairs, as motivated by [30], [32]. For an EI pair, shown schematically in Fig. 1(a), the linear threshold model has the form where a, b, c, d > 0. Following Dale's Law, x E (resp. x I ) is the excitatory (resp. inhibitory) node. The fixed points and attractors of the EI pair have been fully characterized [30], [32], [33] for all possible parameter values. Here we pay attention to the following dynamical features of EI pairs: the origin as a stable fixed point (corresponding to healthy neurological behavior) and stable limit cycles (corresponding to a pathological oscillation characteristic of an epileptic seizure).  (2). The dynamics has a unique stable limit cycle if and only if d + 1 < a − 1 and

Lemma II.1 (Stable equilibria in EI pairs [32, special case of Theorem 2]): Consider the EI pair (2). If the input satisfies
Finally, we give conditions on time-varying inputs to (2) such that the corresponding solution is oscillatory. The result follows as a straightforward consequence of the analysis in [32, Section III].

Lemma II.3 (Sufficient conditions on inputs giving rise to oscillations): Consider the EI pair (2). Let
Then, i) U is nonempty if and only if , then all solutions of (2), except for those corresponding to equilibria, are oscillatory.

III. INTERCONNECTIONS OF EXCITATORY-INHIBITORY PAIRS
We consider complex networks resulting from the interconnection of EI pairs. Consider a network of N nodes, where each node corresponds to a single EI pair (so the state of the network is actually 2 N-dimensional). For i ∈ {1, . . . , N}, the parameters of the ith EI pair in the network are We interconnect the individual EI pairs to form a coupled network, with state space The dynamics of the network are given by (1), where is a weighted adjacency matrix which characterizes the connections between the excitatory nodes of each pair in the network, A II ∈ R N×N ≥0 models connections between the inhibitory nodes of each pair, and connections from excitatory to inhibitory nodes and from inhibitory to excitatory nodes are given by A EI ∈ R N ≥0 and A IE ∈ R N×N ≥0 , respectively. Fig. 3 illustrates a network built by coupling EI pairs in this manner.
We seek to characterize the oscillations in the network using knowledge of the parameters of the individual pairs,  as well as the interconnections A EE , A EI , A IE , A II . We say that the ith node in the coupled network is oscillatory if for all solutions x(t ) of the interconnected system, x i (t ) does not converge to a constant as t → ∞. We say that the ith node is inactive if for all solutions, x i (t ) → 0 as t → ∞. While other behaviors are of course possible, such as converging to a nonzero constant, our focus on these particular ones is driven by their relevance for neurological applications.
Determining whether a particular node in the network is oscillatory or inactive is, in general, nontrivial because of the complex effect of the interconnection on the dynamics of individual nodes. In particular, note that, if the parameters of the ith node satisfy the conditions in Lemma II.1, then the ith node is not necessarily inactive in the coupled network. Likewise if the parameters of the ith node satisfy the conditions of Lemma II.2, this does not necessarily mean that the ith node is necessarily oscillatory in the coupled network. We illustrate these observations in the following example. Although both cases are interesting in their own right, our focus in this paper is on characterizing the robustness of the dynamical properties of the nodes, rather than the emergence of new features through network interconnection. In panel (a), we study a simple network made up of two EI pairs (specific network parameters are reported in (a3)). When taken separately, each EI pair is oscillatory, cf. (a1). In panel (a2), the two EI pairs are interconnected through a simple excitatory-to-excitatory interconnection. As a result of this reciprocal excitation, the two nodes in the interconnected network become saturated.
In panel (b), we consider a network of two EI pairs (specific network parameters are reported in (b3)), where the first node has a globally asymptotically stable limit cycle and the second node has a globally asymptotically stable fixed point at the origin, cf. (b1). After interconnecting the two EI pairs in panel (b2), the second pair becomes oscillatory as a consequence of the incoming activity from the other EI pair.
Given these observations, this paper has two goals. The first goal is to characterize conditions under which properties of  Panel (a) shows a family of seven EI pairs. Colored nodes depict an active EI pair, i.e., a pair which exhibits an oscillatory trajectory (as described in Fig. 1). The hollow nodes represent a inactive EI pair, i.e., a pair whose trajectory is constant and zero. Panel (b) shows the interconnection between the EI pairs according to the adjacency matrix A. As a consequence of this interconnection, node 7 becomes active (the oscillation, we say, spreads to inactive pair 7). The goal of the problems described in Problem 2 is to compute optimal perturbations to the matrix A characterizing how oscillations might (or might not) spread from active to inactive pairs. individual pairs in the network, such as being oscillatory or inactive, are preserved when the pairs are interconnected with one another. The second goal is to develop an approach to modify the network parameters so that given sets of desired nodes are either inactive or oscillatory. We formalize both of these problems next.
Problem 1: What are the conditions on the node parameters W i , u i , and m i , i ∈ {1, . . . , N}, and the interconnection parameters A EE , A EI , A IE , and A II that determine when the dynamical properties of the ith node are preserved after interconnection?
Problem 2: Consider a network (c.f. Fig. 4) whose interconnection is described byÂ EE ,Â EI ,Â IE , andÂ II and let I oscillatory , I inactive ⊂ {1, . . . , N} be disjoint sets of nodes. How should the interconnection structure be modified so that in the resulting network every node in I oscillatory is oscillatory and every node in I inactive is inactive?
Both problem formulations touch upon the spatio-temporal spreading (or lack of thereof) of oscillations in networks and are of general interest. Here, we are particularly motivated by the spreading of microseizures (localized pathological activity in the brain) to clinical seizures (diffused pathological activity in the brain). In this framework, one can interpret Problem 1 as studying whether a given brain network is prone to the insurgence of a clinical seizure as a consequence of a microseizure. Consequently, Problem 2 can be viewed as the study on how interventions on the brain might target specific areas of interest.

IV. SUFFICIENT CONDITIONS FOR PRESERVATION OF DYNAMICAL PROPERTIES OF SUBSYSTEMS
In this section we give a solution to Problem 1 by deriving sufficient conditions which characterize when dynamical properties of subsystems in (7) are preserved after interconnection, i.e., whether the subsystem in question is inactive or oscillatory. The fundamental idea is, for each given node in the network, to (i) view the combined inputs and neighboring signals as a disturbance, and (ii) derive conditions on the interconnection structure to ensure that each node is robust to disturbances. We begin by introducing notation that will be useful for the proofs of the main technical results. The dynamics of the ith node in the coupled network arė where (ũ E i (t ),ũ I i (t )) incorporates the combined input to the ith node from its neighbors: We first present sufficient conditions for the ith node in the network to be inactive. As a consequence of Lemma II.1, the ith node taken individually is inactive when (u E i , u I i ) ≤ 0. The following result gives conditions which ensure that (ũ E i (x),ũ I i (x)) ≤ 0 for all x ∈ X , so the inactivity of the ith node is robust with respect to all inputs the node receives from neighboring nodes.
Theorem IV.1 (Sufficient conditions for robust inactivity): Assume that for i ∈ {1, . . . , N}, Then, the ith node in the coupled system (7) is inactive. Proof: Let x(t ) be a solution to (1). Note that by (10), Since the input to the ith EI pair isũ i (t ) = (ũ E i (x (t )),ũ I i (x(t ))), the result follows by Lemma II.1. Note that the condition in (10) holds only if (u E i , u I i ) ≤ 0. In fact it is not possible for the ith node to be inactive with respect to the interconnected network, unless the ith node individually is inactive.
We now move on to discussing sufficient conditions for the ith node in the network to be oscillatory after interconnection. Our technical approach relies on the following result which formalizes the following observation: the ith node taken individually is oscillatory when (u E i , u I i ) ∈ U i , where U i is the set (4) for the parameters corresponding to the ith EI pair, so the oscillation of the ith node persists after interconnection if

Lemma IV.2 (Robustness of oscillations in coupled networks):
Then, the ith node in the coupled system (7) is oscillatory. Proof: Let x(t ) be a solution to (1). Then by (11), u i (x(t )) ∈ U i for all t ≥ 0. Sinceũ i (x(t )) is the input to the ith node, the result follows by Lemma II.3.
Note that the condition (11) holds only if (u E i , u I i ) ∈ U i , meaning that Lemma IV.2 characterizes only the case where a pair that is oscillatory when viewed individually remains oscillatory after being interconnected in a network. However, as shown in Example III.1, it is possible for nodes (u E i , u I i ) / ∈ U i to become oscillatory after interconnection, though we do not consider such cases in this paper.
In general, the condition in Lemma IV.2 is not easy to check computationally. However, as we show next, in the special case where there are no interconnections from the inhibitory neurons in each subsystem (i.e., A IE = A II = 0) the condition (11) can be expressed in terms of affine constraints on the adjacency matrices.
. In such case, the ith node in the coupled system (7) is oscillatory.
A EI i j m I j ,

FIGURE 6. This plot shows the spreading of oscillations in two 35 × 35 grids of EI pairs, as discussed in Section V-B. Panel (a) shows four consecutive snapshots of the network activity for excitatory (top) and inhibitory (bottom) nodes in each EI pair of a randomly generated grid. We notice how oscillations originate from a subset of pairs, and spread throughout the nominal network. Panel (b) shows the network modified using (14)
, so that all pairs i in I inactive are inactive. We find that 13 out of the 297 edges to nodes in I inactive from nodes not in I inactive have been modified (i.e., changes were made by (14) to less than 5% of the links between EI pairs).

so, using (3), the condition (11) holds if and only if
A EI i j m I j ≥ 0, We obtain (12) by rearranging the above equations.
As the following result shows, in the presence of inhibitory coupling, it is still possible to derive checkable conditions in the form of affine constraints on the entries of the adjacency matrices given by (13). The conditions are valid for all interconnection structures, and are only sufficient for condition (11) to hold. However, in the special case where A IE = A II = 0, the equations (13) reduce to (12) and are both sufficient and necessary for (11) to hold.
Theorem IV.4 (Affine conditions for robust oscillations with inhibitory coupling): Let i ∈ {1, . . . , N}, and suppose that . Then, condition (11) holds and the ith node in the coupled system is oscillatory. Proof: Begin by observing that, and similarly that, It follows immediately that if (13) holds, then (11) holds as well.

Remark IV.5 (Comparison with the literature):
The results presented here differ from the characterization in [30], [33], [39] of oscillations in linear threshold networks in two ways: first, previous results use the lack of stable equilibria as a proxy for the existence of oscillations, whereas we take a slightly different approach by showing when conditions for the existence of oscillations in the ith node are robust with respect to all inputs from neighboring nodes. Second, previous results simply ensure the existence of oscillations, while the results here allow us to determine whether a given node participates in the oscillation, or remains inactive.

V. NETWORK DESIGN USING SUFFICIENT CONDITIONS
In this section we apply the results of Section IV to address Problem 2, and determine how to modify the structure of a given network to control the spread of oscillations. The approach we take is to determine the network interconnection structure as the solution to an optimization problem, where the constraints of the optimization problem correspond to the conditions in Theorems IV.1 and IV.4. We begin by presenting the proposed optimization problems, and then complement them with numerical examples.

A. NETWORK OPTIMIZATION PROBLEMS
Let I oscillatory and I inactive ⊂ {1, . . . , N} be disjoint sets, where I oscillatory denotes the indices of the nodes in the network we desire to be oscillatory and I inactive denotes the indices of the nodes in the network we desire to be inactive. Assume we are given a nominal network of coupled EI pairs whose interconnection is determined byÂ = (Â EE ,Â EI ,Â IE ,Â II ). The problem we address here is to modify the network parameters, so that the modified network given by A = (A EE , A EI , A IE , A II ) has oscillating and inactive nodes as determined by I oscillatory and I inactive , respectively.
We consider two scenarios. In the first, the weights of each of the interconnection matrices can be varied continuously, as is the case with deep brain stimulation, where the magnitude of the change to the network increases with electrical stimulation. In the second scenario, we no longer have finegrained control over the weights of each of the interconnection matrices, and the network structure can only be modified by removing edges entirely. This scenario corresponds to surgical resection, where connections between brain regions are severed to control the spread of oscillations.

1) CONTINUOUS MODIFICATION OF NETWORK WEIGHTS
To find the modified network in the first scenario, we introduce the following optimization problem: The interpretation of the optimization problem (14) is that it finds the minimal modification to the parameters of the nominal network that ensure they satisfy the sufficient conditions derived in Section IV. Because these conditions are affine with respect to A EE , A EI , A IE , and A II , (14) is a quadratic program and can be solved efficiently using standard convex optimization solvers. The solution, A = (A EE , A EI , A IE , A II ), gives the network with the desired properties.

2) MODIFICATION OF NETWORK BY REMOVING EDGES
We now consider the second scenario, where the network can only be modified by removing edges. To solve this problem, we introduce the following optimization problem: Here, denotes component-wise multiplication. In the above problem, S k ∈ {0, 1} N×N , where k ∈ {EE, EI, IE, II} is a matrix where S k i j = 0 if the edge between i and j is severed, and S k i j = 1 if it is preserved. The interpretation of the problem is that it finds the minimum number of edges that need to be severed in order for the modified network to have the desired properties. The modified network is given by (14), the problem (15) is not convex, but rather a mixed-integer program (MIP), which in general is NP-complete. However, modern MIP solvers [40] employ a number of heuristic techniques which allows them to solve relatively large problems efficiently.

3) FEASIBILITY AND OPTIMALITY OF (14) AND (15)
We now give conditions for (14) and (15) to be feasible. Intuitively, both problems are feasible if, in the absence of any interconnection, every i ∈ I inactive is inactive, and every i ∈ I oscillatory is oscillatory. This intuition is formalized in the following result.
Proposition V.1 (Feasibility of (14) and (15)): The problems (14) and (15) are feasible if and only if for all i ∈ I inactive , u i satisfies the hypothesis of Lemma II.1, and for all i ∈ I oscillatory , W i , u i , m i satisfy the hypotheses of Lemma II.2.
Proof: We begin with the forward direction. Suppose that A k = 0 for k ∈ {EE, EI, IE, II}. Then it follows immediately that these matrices satisfy the constraints in (14). Similarly, if S k i j = 0 for all 1 ≤ i ≤ j, then it follows immediately that the constraints of (15) are satisfied.
To show the reverse direction, suppose that A k for k ∈ {EE, EI, IE, II} is feasible for (14). If i ∈ I inactive , then (ũ E i (x),ũ I i (x)) solve (10) for all x ∈ X , which implies that (13) for all x ∈ X , which implies that (u E i , u I i ) ∈ U i . The same argument can be applied for the constraints of (15).
Finally, we show that the optimization problem (14) (14). Then for all i ∈ I oscillatory , the ith node of coupled network determined by A will be oscillatory, and for all i ∈ I inactive , the ith node will be inactive; b) Let A = (A EE , A EI , A IE , A II ) solve (15). Then for all i ∈ I oscillatory , the ith node of coupled network determined by A will be oscillatory, and for all i ∈ I inactive , the ith node will be inactive.

B. SIMULATIONS
We illustrate here how the solutions to the optimization problems above produce network designs that accomplish the desired controlled spread of oscillations.  Fig. 5 shows the results. Note that the nominal network does not satisfy the desired properties since there are nodes i ∈ I inactive which oscillate, and i ∈ I oscillatory do not oscillate and instead saturate at the upper limit. We modify the network designs both by continuously modifying the network weights using (14), and by severing edges using (15). Note that with both designs, the nodes in i ∈ I inactive are inactive, and nodes in i ∈ I oscillatory are robustly oscillatory.

2) EXAMPLE WITH SPATIAL PROPAGATION
For our second example, consider a network of N = 1230 nodes. The first 1225 nodes in the network arranged in a 35 × 35 grid and satisfy the conditions in Lemma II.1. The remaining 5 nodes, whose indices we denote by I driver , satisfy the conditions in Lemma II.2. The network has only excitatory-excitatory coupling, where each node in the grid is coupled to the nodes in the cells above, below, to the left and to the right, and the driver nodes are randomly coupled to nodes in in the grid (Â EE denotes the adjacency matrix of the network). As shown in Fig. 6(a), the oscillations from the 5 driver nodes spread throughout the network.
We now mark a region in the grid where we do not want the oscillations to spread. Let I inactive denote the indices of these nodes. To determine how to modify the network interconnection, we use the optimization problem (14) to determine the interconnection weights A EE of the modified network. As shown in Fig. 6(b), the oscillations in the modified network spread from the driver nodes to other regions in the network while avoiding the region determines by the nodes in I inactive .

VI. CONCLUSION
We have studied the spreading of oscillations in complex networks of interconnected EI pairs. Such networks, modeled here with a piecewise-linear activation function, are an expressive modeling tool for oscillatory networks, with meaningful connections to brain dynamics. Motivated by the link between the spatial spread of oscillations and seizures in the brain, we have identified formal conditions on the network interconnection structure that determine which regions oscillate and which remain inactive. We have built on this understanding to propose strategies that mitigate the spread of oscillations among brain regions. These strategies formalize network design objectives by means of optimization programs with attractive numerical properties. The simulation results show that the proposed approach may be effective in practice. Future work will address the study of the emergence of oscillations through network interconnection and in particular characterize conditions for the nodes in networks of coupled EI pairs to be oscillatory without assuming that they admit oscillations individually. Further, we hope to compute tighter bounds on the sufficiency conditions for robust oscillations and to further validate these results through data-generated models of brain networks from human subjects.