Distributed control in a mean-field cortical network model : Implications for seizure suppression

Brain electrical stimulation (BES) has long been suggested as a means of controlling pathological brain activity. In epilepsy, control of a spatially localized source, the seizure focus, may normalize neuronal dynamics. Consequently, most BES research has been directed at controlling small, local, neuronal populations. At a higher level, pathological seizure activity can be viewed as a network event that may begin without a clear spatial focus or in multiple sites and spread rapidly through a distributed cortical network. In this paper, we begin to address the implications of local control in a network scenario. To do so, we explore the efficacy of local BES when deployed over a larger-scale neuronal network, for instance, using a grid of stimulating electrodes on the cortex. By introducing a mean-field model of neuronal interactions we are able to identify limitations in network controllability based on physiological constraints that suggest the need for more nuanced network control strategies.


I. INTRODUCTION
Pathological oscillations are characteristic of many neurological disorders. A striking example of this is epilepsy, the condition of chronic unprovoked seizures, which affects over 50 million people worldwide [1]. For many of these patients, seizures remain poorly controlled despite maximal medical management. In these cases, invasive treatment options are a viable alternative, including resective surgery, in which the brain region implicated for seizure genesis is removed [2]. When the epileptogenic zone includes eloquent cortex (e.g., motor or speech cortex), alternative invasive treatments are considered, including brain electrical stimulation (BES) [3]. The most common BES method is vagus nerve stimulation, thought to affect brain regions (e.g., the thalamus) that might, in theory, increase cortical inhibition and thereby lessen or modulate seizures [4]. Seizures may also be aborted by direct electrical stimulation of the cortex. For example, cortical afterdischarges (seizure-like activity elicited by direct cortical electrical stimulation) may be arrested by applying brief bursts of pulse stimulation to the cortex [5] at a variety of frequencies [6,7].
At the microscopic scale of individual neurons and small neural populations, the impact of electrical stimulation on seizure-like activities has been thoroughly studied [8,9]. But the effect of BES on a network of macroscopic brain regions remains incompletely understood. Characterizing the impact of cortical stimulation at the macroscopic level will be important in determining optimal stimulation sites, transducers, and parameters for terminating seizures.
Animal models of epilepsy have allowed researchers to explore different BES methods that currently pose unacceptable risks to human subjects [9,10]. In addition, mathematical models provide an alternative means to safely investigate the effects of new BES paradigms [11]. Typical studies have focused on the mechanisms of control in local neuronal populations, for instance, cortical columns [12,13], single cells [14], and small neuronal populations [15,16]. While these works have offered insights into local controllability, none directly address the notion of distributed control of the seizure activity propagating through a large network. Such an understanding is particularly important to treating epilepsy in which some seizures can be understood as network events, beginning in a circumscribed region but rapidly recruiting other brain areas through a cascade of spreading activity. An optimal control strategy may thus target multiple spatially localized regions while being sensitive to the network structure.
This paper offers an initial theoretical exploration of control issues in cortical networks. We are motivated by existing technologies, such as electrocorticography (ECoG), in which grids of electrodes are placed on the cortical surface [see Fig. 1(a)]. Animal experiments and medical technologies such as deep-brain stimulation show that voltage stimulation can suppress activity in neural tissue [9]. Thus, it is plausible that such a grid could be used in a distributed control scheme. However, would localized suppression at a few spatially disparate regions be sufficient to stabilize the entire network? In order to investigate this question at the spatial scale relevant to ECoG, we use a mesoscale model consisting of coupled mean-field neuronal oscillators [17]. The propagation of seizure activity was previously modeled in a similar context in Ref. [18], where a controller was proposed to modulate a localized epileptogenic focus. Here, we examine distributed control at multiple sites. The model is simple enough to be analytically tractable but maintains some important features from neurophysiology. By characterizing the feasibility of control in this model, we can identify limitations that would persist in more complex biophysical models.
To this end, we study two types of cortical networks: those coupled by diffusion and by synaptic transmission. The former has been employed in mean-field modeling (e.g., [19]) and permits analytic characterizations, although the connection to neural physiology is less clear. Here, we use diffusive coupling to gain basic analytic insight before moving to the case of synaptic coupling, which is the commonly accepted modality by which most corticocortical transmission occurs, although this improved realism limits rigorous analysis. The controllability of the network is different in the two cases. In the diffusive network, a sufficiently dense grid of local control nodes can always suppress pathological activity. Conversely, when coupling is synaptic, distributed control is not sufficient to completely suppress pathological activity but may be able to firewall activity within a section of cortex. Consequently, the analysis demonstrates fundamental limitations that suggest the development of new BES network paradigms.
The outline of this paper is as follows: In Sec. II we introduce the mean-field cortical network model. In Secs. III-V we describe the controllability of such a network in the case of diffusive and synaptic connectivity. For the latter we consider both nearest-neighbor and small-world network topologies. Conclusions and an outline of future work are provided in Sec. VI.

A. Cortical network model
The Wilson-Cowan model [20] is used to describe oscillatory behavior in a single cortical macrocolumn, i.e., a cylindrical portion of cortical tissue ∼0.5-3 mm in diameter [21]. It is a mesoscale model and offers a compromise between the biophysical detail of single neuron models of the Hodgkin-Huxley type [22] and the ability to describe macroscopic brain activity, such as field potentials, in large populations and over large cortical distances.
Here, we construct a network of N Wilson-Cowan columns, given bẏ where (e j ,i j ) are, respectively, the activity in excitatory and inhibitory cell populations of the j th column. The function F is of the standard sigmoidal form The biophysical interpretation and typical values for the parameters in Eq. (1) are given in Table I. The network interconnectivity arises through the functions C e (·) and C i (·), whereē and the connections between columns are only excitatory. The term u j (t) is the exogenous input, or control, while w(t) is a white noise process that models random background activity in each column. If column j is actuated, for instance, by a surface electrode, then b e j ,b i j are nonzero. Again, we are motivated by existing electrode grid configurations, e.g., Fig. 1, in which only a subset of cortical columns are directly stimulated. For most of the simulations we will assume that the control uniformly deactivates both excitatory and inhibitory populations, i.e., b e j = b i j < 0. In Sec. VI we will consider the issue of imperfect control, where these coefficients may be uncertain due to, for example, differential effects of stimulation depending on neuron orientation [8,23].
In the absence of both interconnections, i.e., when C e (·) = C i (·) = 0, and control, i.e., when b e j = b i j = 0, a single column exhibits two dynamic regimes of interest, depending on the level of external input P j (t),Q j (t) to the excitatory and inhibitory populations, respectively. Specifically, when P j (t),Q j (t), and w(t) are sufficiently small, the column displays small deviations about a stable equilibrium [note that (e j ,i j ) = (0,0) is a stable equilibrium when P (t),Q(t) = 0]. The maximum values of the excitatory, inhibitory response functions 1,1 r e ,r i The absolute refractory period of the excitatory, inhibitory subpopulations 1,1 a e ,a i The value of the maximum slope of the logistic curve for the excitatory, inhibitory subpopulation The position of maximum slope of the logistic curve for the excitatory, inhibitory subpopulation 4, 3.7 Such behavior is representative of normal neuronal population activity. Conversely, when P j (t) is large, the column displays oscillations in the (e j ,i j ) plane. Such an oscillation is more indicative of a pathologically entrained neuronal ensemble and will be used herein as a surrogate for epileptogenic activity.

B. Control in a single column
Previous studies in BES, both experimental and computational, have typically focused on the problem of control in an isolated neuronal population [3,8,24]. Emphasis has been placed on understanding the specific mechanisms of electrical stimulation, including the sensitivity of cells to electric fields [25] and their orientation with respect to the cortical surface [23]. From a control perspective, effort has been directed at addressing practical considerations such as minimization of delivered current and charge balance [3,26].
There remain many open questions with respect to singlecell and small-population control. Our focus here is on how such control, once achieved, may affect activity at the network level. Hence, we assume that there exists a controller design that successfully regulates the neuronal activity of a single cortical column in real time, i.e., a local control of the form where g(·) maps the excitatory neural activity into y j , the observed brain voltage signal. For the purposes of our simulations, we will assume, without loss of generality, that g(e j ) = e j and h(y j ) = ky j , i.e., simple proportional control based on (1) and (2). Such a scheme, though simplistic, is sufficient to ensure suppression of an individual cortical column in a more detailed model [12] and, more importantly, enables us to assess the implications of local control in a network of interacting elements. In reality, the specific local control law would likely rely on desynchronization of the underlying population for achieving suppression [16,27,28].

III. CONTROL IN DIFFUSIVELY COUPLED NETWORKS
The simplest model of coupling between neuronal populations is standard diffusion. Here, the terms C e (·) and C i (·) in (1) and (2) take the form where N denotes the set of neighboring columns and k d is a diffusive constant. In the absence of input, the controllability of the network about the origin can be assessed directly by linearizing (1) and (2) about the point (ē,ī) = (0,0). The more interesting case occurs when P j (t) is nonzero and the network exhibits pathological oscillation. Consider a simple two-column network, reciprocally coupled as shown in Fig. 1(b), and let P 1 (t) be such that both columns exhibit limit cycle behavior. Assume that control is applied locally to the second column in order to suppress activity there, i.e., b 1 = 0,b 2 = 1, so that FIG. 2. Two columns with diffusive coupling. External input to column 1 drives oscillations in both columns. Control applied to column 2 at t = 400 ms eliminates oscillations in both columns. Proposition 1. For k d sufficiently large, local control in column 2 suppresses limit cycles in the connected two-column network.
See the Appendix for the proof to Proposition 1. The intuition behind this result follows from the fact that the diffusive term, k d e 1 , acts as a proportional feedback for column 1.
To illustrate control in the two-column network, consider the standard parametrization [20] with w(t) being a Gaussian random process of variance 0.1. Figure 2 plots e 1 ,e 2 from the resulting simulation when control at column 2 is turned on at t = 400 ms. As shown, this control induces stability about a low-activity fixed point in both columns. In general, for certain parameter ranges, a column may exhibit multiple attractor states [29]. However, this phenomenon will not occur for the parametrizations considered in this paper. Proposition 1 suggests that, in a diffusively coupled network of more than two nodes, a sufficiently dense grid of control electrodes can suppress the generation and propagation of pathological activity (consider the limiting case of a control actuator placed at every other node in the network). To demonstrate this property, we consider a simulation of a network consisting of 225 columns arranged in a 15 × 15 sheet. Here, we assume that a column encompasses a volume that is roughly 6 mm in diameter, or the equivalent of roughly four macrocolumns, which we reason is approximately the region affected by a typical intracranial electrode. In this setting, our network spans a surface area of 9 cm × 9 cm, which is approximately the size of a typical intracranial electrode grid. On this network, we will consider 25 (5 × 5) equally spaced actuators, equivalent to roughly a 1-cm interelectrode spacing, again, typical of actual intracranial grids. Thus, our configuration equates to a grid density (the ratio of controlled to uncontrolled columns) of around 1 : 10. to a column in the lower right quadrant [P j (t) = 1.25 at the column], which kindles and entrains the rest of the network. Nodes at the network boundary receive input from neighboring nodes but do not provide outgoing excitation to other nodes. Figure 3(b) shows the same network with control applied using the grid of 25 actuators.
Localized control at these sparse columns is sufficient to stabilize the network. Figure 3(c) shows a Monte Carlo analysis (n = 100 simulations) of the entrainment (as defined by the sum of excitatory activity over all columns) of the network in both the uncontrolled and controlled cases. During each simulation trial, a new realization of the background noise w(t) is generated, and the model parameters are randomly and uniformly chosen to within 5% of the values in Eq. (8). Over time, the controlled network exhibits much lower entrainment than the uncontrolled network.
While distributed control here seems effective, the diffusive coupling model has an important deficiency. It assumes passive inhibition, i.e., quiescence in column 2 suppresses column 1. Such inhibition may occur in the brain via ephaptic means or through gap junctional connections. The dominant mode of transmission in the cortex is, however, synaptic coupling, in which the presynaptic area must be excited in order to modulate activity in the postsynaptic area. In the next section we explore the limitations of networked control in this more realistic scenario.

IV. CONTROL IN SYNAPTICALLY COUPLED NETWORKS: NEAREST-NEIGHBOR TOPOLOGY
Consider a simple representation of synaptic coupling given by where N again denotes the set of neighboring columns and k s is the coupling strength [30]. Here, excitation in a neighboring column induces activity in both e and i populations. The following result highlights a fundamental limitation in the controllability of such a network Proposition 2. Consider the network of two columns illustrated in Fig. 1(b) and described by (1) and (2), coupled reciprocally via synaptic coupling (9), where P 1 (t),Q 1 (t) is such that limit cycles are displayed in both columns in the absence of control. Unlike the diffusive model, the system is not controllable through control applied to the second column.
See the Appendix for the proof of Proposition 2. Figure 4 provides an illustration of this property, showing the same two-column network as in Fig. 2, except with (9) and k s = 2. Although column 2 is locally suppressed, column 1 continues to display a limit cycle.
Note that the exact size of an epileptic focus is, in general, not well known and may be larger than a column. Proposition 2 is intended only to illustrate that nonlocal control, offset from pathological tissue, cannot yield complete suppression of activity in the network. Although complete suppression is not possible, distributed control may still limit, or firewall, pathological activity. As an example, consider Fig. 5(a), which shows propagation of pathological activity in an uncontrolled, synaptically coupled network. The parametrization here is the same as in Fig. 3, except with P (t) = 0.8,k s = 1.5 for all columns. Again, excess drive is initially applied to a single column in the lower right quadrant [P j (t) = 1.25 at the column], which is sufficient to kindle and entrain the remainder of the network. When local control is applied at 25 equally spaced columns FIG. 4. Two columns with synaptic coupling. External input to column 1 drives oscillations in both columns. Control applied to column 2 (at t = 400 ms) fails to control column 1. [ Fig. 5(b)], the pathological activity stays firewalled within the region of kindling.
We note from Fig. 5(c) that, over time, the pathological activity may still begin to accumulate within the network. Indeed, depending on the strength of coupling, the firewalling effect may be eventually overcome by recurrent excitatory activity. Nevertheless, containing or slowing propagation may still be valuable as a tool to detect or localize seizure activity, providing added time for clinical intervention with pharmacological or other therapeutic agents. Moreover, slowing seizure propagation may allow intrinsic mechanisms within epileptic foci to "run their course" and terminate endogenously without affecting larger areas of the brain. Remark 1. In the two-column network (Fig. 4) one may surmise a strategy in which column 2 is placed in an antiphase oscillation to that of column 1, providing a sort of destructive interference. While this may work in suppressing column 1, such an antiphase oscillation requires activation of column 2, which may entrain pathological activity in other columns activated by column 2 in a larger network.

V. CONTROL IN SYNAPTICALLY COUPLED NETWORKS: SMALL-WORLD TOPOLOGY
The nearest-neighbor connectivity topology, while simple, may not accurately reflect the longer-range anatomical con-FIG. 6. Control efficacy in small-world topologies: (a) Control efficacy when, with probability p, one nearest-neighbor edge is replaced with long range connection. (b) Propagation and control efficacy in a small-world network topology with higher connection density as compared with purely nearest-neighbor coupling. A more dense (i.e., controlled nodes to uncontrolled nodes) control grid is required to recover comparable propagation suppression performance. Lines indicate mean of n = 100 simulations, with parameters fixed and different instantiations of small-world topology and noise.
nections within the cortex [31,32]. To understand the effects of these longer-range connections, we perform simulation studies in a small-world network topology [33]. Here, each column may, with some probability p, make a long range connection to a non-neighboring column. To maintain the overall connectivity density, when a long range connection is made, one of the nearest-neighbor connections is randomly and equiprobably removed. Figure 6(a) demonstrates the control efficacy of this scenario for various values of p. As shown, the performance of the control scheme is approximately retained over all cases. The result can be understood in terms of the uniformity of our control grid. Since we use equally spaced actuators that cover the entire network, any small-world edge emanating from the focus is effectively firewalled by an actuator located near the termination of the connection. If our grid was not uniform and only covered the region of the network containing the focus, then the presence of small-world edges would cause the activity to escape containment. This result suggests that, for a uniform grid, the critical factor in efficacy is not the degree of randomness but instead a ratio of the strength and density of connections between columns to the spacing or density of actuators.
To examine this issue, we also study a small-world situation in which connection density is not maintained. In this setting, when a long range connection is made, a nearest-neighbor connection is not removed. Figure 6(b) illustrates the propagation of activity in this setting for p = 0.1, where the simulation parameters are otherwise identical to that of Fig. 5. Not surprisingly, the uncontrolled propagation is faster than in the nearest-neighbor-only configuration, and the efficacy of the 5 × 5 control grid is diminished. When the grid density is increased to 6 × 6 (relative to the same 15 × 15 network), the control efficacy improves substantially, with the propagation approximating that of the 5 × 5 case in Fig. 5(a). The results suggest that random long range connections exacerbate the speed of activity propagation in this model only when they increase the overall connection density, and that this effect can be offset by a compensatory increase in grid density.
We can thus interpret the effect of our distributed control scheme as dependent on two complementary factors: (i) a spacing, or density, of actuators sufficient to contain the spread of excitatory activity and (ii) coverage by our actuator grid of those regions affected by long distance connections emanating from the seizure focus. The combination of these two factors leads to what we have termed the firewalling of seizure activity. When firewalled, the activity remains localized in the region of the seizure focus, even in the case of long range connections. Thus, given a set of known long range connections (from, say, brain mapping studies [34]), along with nearest-neighbor local connectivity, the model could be used to investigate the minimum actuator spacing or density required to achieve propagation containment. Moreover, such studies may suggest the design of nonuniform grids that place electrodes only near areas impinged upon by the seizure focus. Characterizing in detail the relationship between network connection strength and actuator grid spacing and design will be the subject of future work.

VI. LIMITATIONS DUE TO IMPERFECT LOCAL CONTROL
As mentioned above, the methodology introduced in this paper is intended to assess distributed control over a large cortical network. We do not explicitly treat the problem of local control, which is itself an active area of research. Issues such as the geometry of the control electrode with respect to the cortical surface may significantly alter the efficacy of any fixed feedback scheme of the form (5). For instance, depending on the orientation of the sulci near the electrode, the stimulation may preferentially affect the inhibitory interneurons as opposed to the excitatory pyramidal cells [8]. Similarly, the same stimulation parameters may depolarize, rather than hyperpolarize, some elements of the underlying cell population. In our model, this would amount to a change in the coefficients b e j and b i j in Eqs. (1) and (2). To investigate this issue, we examine the effect of imperfect local control on the distributed control scheme. Specifically, we repeat the simulation of Fig. 5 by randomly choosing a portion of the control network (10%, 30%, or 50% of nodes) for which control is imperfect, modeled here by setting b e j or b i j to zero, i.e., control of only the inhibitory or excitatory population.
As shown in Fig. 7, the critical factor in the efficacy of the scheme is the control to the excitatory cell population. If control is preferential to the inhibitory population (i.e., local failure causes b e j = 0), then the performance degrades as a function of the failure rate [ Fig. 7(a)]. Functionally, this is analogous to a decrease in the grid density. Moreover, by selectively suppressing inhibitory elements, the propagation in the network is actually enhanced due to an excess of excitation compared to inhibition (compare with Fig. 5). Conversely, if control is preferential to the excitatory population (i.e., local failure causes b i j = 0), performance is unaffected [compare Fig. 7(b) with Fig. 5]. The result is promising since excitatory cells, owing to their longer axons [35], are more apt to be affected by electrical stimulation.
We do not consider the case when the local electrode may depolarize the underlying population, i.e., b e j ,b i j > 0. Uncertainties in the sign of b e j ,b i j would likely necessitate specialized local controllers in place of (5). For instance, a unique situation may arise when b e j ≈ 0, b i j > 0. Here, in order to suppress excitatory activity, one may surmise a local strategy that depolarizes the inhibitory cell population.

VII. DISCUSSION AND CONCLUSIONS
Understanding the biophysical mechanisms of BES remains an active and important area of research. While elucidating these mechanisms, it will be important to consider the implications of local BES when deployed over broader networks. This paper has introduced a modeling framework for this purpose, a platform on which to bridge analysis at the level of small neuronal populations with that over larger brain areas.
We show that local control may be insufficient to completely suppress pathological activity, which may be important when such activity is diffuse in space or does not have a predictable focus. This is due to a basic limitation in the controllability of the underlying dynamical system. Nevertheless, distributed control has the capacity to firewall or delay propagation, which may prove useful for real-time seizure detection or localization algorithms, potentially enabling intervention with other therapeutics. The controller would only need to remain on until such time as these other interventions were made (or such time as the seizure dynamics ceased endogenously).
The size of our model network has been chosen to match the scale of typical intracranial electrode grids. Based on our analysis, we expect that the qualitative behavior of the distributed control scheme does not change with increased network size as long as the interelectrode distances are maintained. Nevertheless, in future work, the model may be paired with anatomical data from brain mapping studies to examine larger networks and account for additional complexities in cortical connectivity, for instance, long range excitation from column a onto only inhibitory cells in column b. When combined with more principled methods for differentially targeting inhibitory or excitatory populations, this additional detail may suggest entirely new schemes for spatially containing propagation. Such schemes have the potential to overcome the limitations illustrated herein.
Evaluation in more detailed computational models and experiments must be completed before any possible implementation. However, the model shown herein can be used to examine theoretical yet fundamental issues such as control grid density, precision of electrode placement, and parameter configuration. By rapidly synthesizing potential control strategies, the model can serve as a precursor to important, but challenging, experimental studies in the future. ACKNOWLEDGMENTS M.A.K. and S.C. acknowledge support through Career Awards at the Scientific Interface from the Burroughs Wellcome Fund. S.C. and E.N.B. are supported in part by NIH DP1-OD003646.

APPENDIX
Proof of Proposition 1. Since column 2 is controlled, it follows that e 2 (t) can be made arbitrarily small, and thuṡ e 1 = −e 1 + (k e − r e e 1 )F e [(c 1 − k d )e 1 − c 2 i 1 + P 1 (t) + ], (A1) where P 1 (t). As shown through structural bifurcation analysis in Ref. [29], for any P 1 (t),Q 1 (t), there exists a k d > 0 such that no limit cycles exist and the system exhibits stability about a fixed point.
Proof of Proposition 2. The proof is by contradiction. Specifically, assuming the converse and using the standard definition of controllability, there exists an input u 2 (t) such that [e 1 (t),i 1 (t),e 2 (t),i 2 (t)] = 0 for t > 0. Since e 2 (t) is suppressed, column 2 sends no outgoing synaptic excitation to column 1. Then it follows that the dynamics of column 1 arė e 1 = −e 1 + (k e − r e e 1 )F e [c 1 e 1 − c 2 i 1 + P 1 (t)], (A3) which, by definition, exhibits a limit cycle. In this case, [e 1 (t),i 1 (t)] = 0, and the system cannot be controllable. Note that a full dynamical analysis of this configuration is contained in Ref. [30].