Optimal control of excitable systems near criticality

Experiments suggest that the cerebral cortex gains several functional advantages by operating in a dynamical regime near the critical point of a phase transition. However, a long-standing criticism of this hypothesis is that critical dynamics are rather noisy, which might be detrimental to aspects of brain function that require precision. If the cortex does operate near criticality, how might it mitigate the noisy ﬂuctuations? One possibility is that other parts of the brain may act to control the ﬂuctuations and reduce cortical noise. To better understand basic aspects of controlling neural activity ﬂuctuations, here we numerically and analytically study a network of binary neurons. We determine how the efﬁcacy of controlling the population ﬁring rate depends on proximity to criticality as well as different structural properties of the network. We ﬁnd that control is most effective—errors are minimal for the widest range of target ﬁring rates—near criticality. Optimal control is slightly away from criticality for networks with heterogeneous degree distributions. Thus, while criticality is the noisiest dynamical regime, it is also the regime that is easiest to control, which may offer a way to mitigate the noise. DOI


I. INTRODUCTION
A large body of experimental evidence supports the hypothesis that the cerebral cortex operates at, or very near, the critical point of a phase transition [1][2][3][4][5].Consequently, many efforts have been made to identify and understand the functional benefits of operating in this regime [6], including maximized dynamic range [7][8][9] and maximized information transmission and capacity [10][11][12].Both experiments and models have shown that networks operating in the critical regime visit the largest variety of microscopic and macroscopic states [10,13,14].It has been hypothesized that these states could be harnessed by the network to encode and transfer information [15].In this paper, we identify a way that criticality may be beneficial in neural systems.We show that, at criticality, the activity of the system can be controlled with minimal error over the largest range of activity levels.In addition, by analytically treating the network's deviations from linear dynamics, we show that heterogeneity in the network's degree distribution reduces this controllable range considerably.Understanding the controllability of a neural system may be important for designing effective therapeutic brain stimulation, closed-loop brain-machine interface systems, and understanding basic experimental studies of optogenetic control of neural activity [16,17].
There has been a significant amount of work in quantifying the controllability of functional brain networks and networked dynamical systems.The main methodology for studying control of brain networks has been to assume an approximate linear evolution model for the node state vector s t on a network derived from functional Magnetic Resonance Imaging (fMRI) experiments (e.g., [18]) or from complete connectomes [19] of the form s t+1 = As t + u, where u is the control signal and A encodes the network connections.Using established results from control theory, one can determine the ability of the control signal to drive the state vector to a desired state.Such analyses have provided valuable insights, such as the role of specific brain regions in driving transitions between different states of brain activity [18], the functional role of individual neurons in controlling the locomotion of Caenorhabditis elegans [19], and the role of large-scale network topology on the efficiency of this type of control [20] (for a review and other applications, see [21]).Other examples of control of networked systems include pinning control of synchronization [22] and synthetic rescue in metabolic networks [23].For a review of control in complex networks, see [24].In the models that assume a linearized evolution of the node state vector, the Perron-Frobenius eigenvalue of the matrix A, λ, is typically assumed to be fixed at a value less than 1 to guarantee stability of the uncontrolled system.In contrast, here we are specifically interested in the relation between control performance and the proximity of the system to the critical state where activity is neither amplified nor damped, i.e., to λ = 1, and therefore, we consider λ as our main parameter of interest.In addition, we consider nonlinear evolution of the form s t+1 = σ (As t + u), where σ is a nonlinear function.Finally, we note that our interest in this paper is not on the ability of the control signal to drive the state vector to a target state but in minimizing the error when trying to maintain a macroscopic variable (e.g., a scalar function of the state vector) as close as possible to a desired value.
The structure of this paper is as follows.In Sec.II we present our model and the control task.In Sec.III we present the analysis of the control error for both the linear and nonlinear regimes.In Sec.IV we study the stability of the control scheme, and in Sec.V we present our conclusions.

II. MODEL
We first introduce the model without control.Following Refs.[7,13,14], we consider a network of N discrete-state excitable nodes.At each time step t, each node m may be active or inactive, corresponding to a node variable s t m with a value of 1 or 0, respectively.The nodes are coupled via a directed, weighted, and signed network with adjacency matrix A, such that A nm represents the strength of interaction from active node m to node n.These values are non-negative (nonpositive) whenever node m is excitatory (inhibitory), and networks are assumed to consist of 20% inhibitory nodes.The dynamics are probabilistic, with node n becoming active at time step t + 1 with probability σ ( N m=1 A nm s t m ), where σ is a piecewise linear transfer function defined as σ (x) = 0 for x < 0, σ (x) = x for 0 x 1, and σ (x) = 1 for x > 1.In such systems, the largest eigenvalue λ of the adjacency matrix A determines whether the system is in a regime of low activity (λ < 1) or high activity (λ > 1) or at criticality, the boundary between high and low activity regimes (λ = 1).In past work, this model (or slight variations of it) has been used to explain how the rich population dynamics of cortical neural networks depends on changes in excitatory and inhibitory interactions and several aspects of how such networks process and adapt to changing external stimuli [3,5,[8][9][10].
We now introduce control to the dynamics and analyze how its effectiveness depends on the network's topology and principal eigenvalue λ.To this end, we consider a modified evolution equation, and s t+1 n = 0 otherwise.The term μ n ( Ŝ − N m=1 b m s t m ) represents an external proportional control signal which attempts to bring the global variable S = N m=1 b m s t m to a target value Ŝ, where b n represents the weight given to node n in the global activity S. For example, this control term in the model may be interpreted as optogenetic control mechanisms in recent experiments [16].The dependence of the control term μ n on n allows one to include the effects of natural heterogeneities, experimental constraints, or targeted application of the control signal.Similarly, the dependence of the weights b n on the node n allows for the possibility that the measured nodal activities are given different weights or that one controls the activity of a subset of nodes.
There has been a significant amount of previous work on studying the control of brain dynamics from the point of view of classical control theory [21].For example, by assuming linear dynamics on functional brain networks, one can study the question of whether or not these networks are globally controllable, i.e., whether the dynamical state of the network can be steered to any desired state in finite time by a control signal [18].In contrast, here we consider only the problem of keeping a macroscopic variable close to a desired target value.Purposefully, this task is simpler than those considered in those previous works; this simplicity allows us to focus on the effects of the nonlinearity, the network's degree distribution, and the excitation-inhibition balance.

A. Linear regime
As a first step, we will ignore the nonlinearities of the function σ by assuming that its argument is always between 0 and 1, the regime in which σ (x) = x.Later, we will discuss the conditions under which this assumption is valid and consider the effects of the nonlinearity.
In this linear regime, each node's expected activity p t n ≡ E [s t n ], taken over the stochastic dynamics, is, at steady state, where we have introduced vectors m, b, and p with entries μ n , b n , and p n , respectively.This can be solved for p using the Sherman-Morrison inversion formula to get Equation ( 3) predicts the steady-state averages of nodal activity p n = E [s t n ] in terms of c and therefore in terms of the network A and the control vector m.
Our goal is to determine how control error depends on the network's degree distribution and proximity to the critical point λ = 1.We begin with analytical predictions and then test these predictions with numerical simulations.
We first assume that the in-and out-degrees, k out n = N m=1 A mn and k in n = N m=1 A nm , and the largest eigenvalue λ are specified.Then we consider the "annealed network" approximation [25][26][27] where we replace A by an averaged adjacency matrix Ā which preserves the eigenvalue λ and the degree distribution in the sense that which has the eigenvector k in with eigenvalue λ.Using this matrix, we effectively average over the ensemble of networks with largest eigenvalue λ and the desired degree distributions.This approximation assumes that there is no assortative mixing, community structure, or other additional structure in the network.While these features are always present to some extent in real-world networks, we focus here on the simplest case in which they are absent.
Replacing A with Ā in the definition of c and using the Sherman-Morrison inversion formula again, we get We quantify control error as the expected relative error Substituting Eq. ( 3) for p, we find R = −1/(1 + b c).Then, inserting Eq. ( 4) for c yields This equation makes two key predictions for the control error R based on the network's degree distribution (k in , k out ), the control vector m, the nodal weights b, and the principal eigenvalue λ.First, Eq. ( 6) predicts that control error is minimized to R = 0 whenever λ = 1, independent of the target Ŝ.Second, Eq. ( 6) predicts that when λ = 1, the magnitude of the nonzero control error depends on correlations between the in-and out-degrees of the nodes, the weights b, and the nodal control strengths m.We confirm these predictions via simulation by varying both the principal eigenvalue of the network adjacency matrix and the correlation between in-and out-degree sequences.Figure 1(a) illustrates the effects of λ on dynamics with and without control, with sample time series of network activity S t = N n=1 s t n /N for simulated dynamics with λ = 0.98, 1.0, and 1.02 on directed Erdös-Rényi random networks (simulation details are given in the caption).6) (solid lines) compared to simulations (squares and triangles) in which we systematically varied λ between 0.9 and 1.1 for two classes of networks which were designed to maximize and minimize the effects of degree correlations, respectively.To construct such networks, we considered a sequence of target degrees k1 > k2 > .We refer to the second (first) case as the (anti)correlated case since there is a positive (negative) correlation between kin n and kout n .We generated adjacency matrices A by choosing the entries A nm to be 1 with probability k in n k out m / n k out n [28] and then rescaled the resulting matrices so that they had the desired eigenvalue λ.We emphasize that while the simulations in Fig. 1(b) were done with specific network realizations, the predictions of Eq. ( 6) were derived by considering an averaged network.Finally, we note that analogous expressions for the control energy in terms of network control vectors can be derived for the related problem of reaching a given dynamical state (see, e.g., [29]).

B. Nonlinear regime
In the two cases described above, the control error is well described by Eq. ( 6), but this is not the case if the network's degree distribution is heterogeneous (even mildly).For example, if the degrees kn are sampled from a power-law distribution, P(k) ∝ k −γ , with a minimum degree of 50 and exponent γ = 4 [circles in Fig. 1(b)], the control error is no longer minimized at λ = 1.To understand this discrepancy, we assume for simplicity that each node is weighted equally, i.e., b n = 1/N, and consider the evolution of the network activity S t = N n=1 s t n /N.In the absence of control input and nonlinearity, the expected network activity for a homogeneous network with largest eigenvalue λ evolves, given S t , as E [S t+1 ] = λS t ; that is, the largest eigenvalue λ acts as a branching parameter for the fraction of excited nodes.Thus, our linear results above can be understood as follows: when λ = 1, S t evolves stochastically with no bias to increase or decrease (i.e., E [S t+1 ] = S t ), and therefore, control can effectively pin it to any value of the target Ŝ.When λ is larger (less) than 1, S t has a tendency to increase (decrease).To establish a (statistical) steady state, this tendency needs to be balanced with a nonzero control term, which can be accomplished only if Ŝ − S = 0, i.e., if there is control error.When the nonlinearity of the transfer function is taken into account, the relationship E [S t+1 ] = λS t is no longer valid, but one can still define an effective branching ratio, the branching function , that given S t , satisfies E [S t+1 ] = S t : The branching function (S) can be interpreted as an effective eigenvalue or branching ratio that applies when S t = S: when (S) > 1 [ (S) < 1], activity tends to increase (decrease) on average.In Ref. [13] it was shown that for directed FIG. 2. Experiments with three networks with N = 5000 show the validity of Eqs. ( 12) and (8).Network degree distributions were either uniform in [100,200] (solid blue line, λ = 1) or power laws with γ = 4 and a minimum degree of 50 (dot-dashed green line, λ = 1; dashed red line, λ = 1.05).(a) Branching functions calculated using Eq. ( 12) for the three networks.(b) Control error measured numerically (symbols) and theoretically using S obtained from Eq. ( 8) (lines) as a function of target Ŝ for the same three networks.Note that error is small whenever the branching function is close to 1.
Erdös-Rényi networks, depending on the fraction of inhibitory neurons and the average network degree, the branching function is a decreasing function of S with (S) λ for small S, (S) ≈ λ for intermediate S, and (S) < λ for S ≈ 1.The region where (S) ≈ 1 starts approximately at S ∼ 1/ k , where k = i k i /N is the mean degree of the network [13].If the mean degree is small, this region could vanish.If there is no inhibition, (0) = 1.Branching functions for three different networks are shown in Fig. 2(a).
Because the statistical steady state of S t is still established by a balance between its tendency to increase [if (S) > 1] or decrease [if (S) < 1] and the control term, we conjecture that Eq. ( 6) is a good approximation for the control error in the nonlinear case if one replaces λ by (S), where S is found self-consistently from the implicit equation that defines the error, Below we provide numerical evidence that shows the usefulness of this heuristic approach to predict how the controllable range is affected by heterogeneity in the network's degree distribution.
Our goal now is to estimate the branching function in terms of the network's degree distribution and the transfer function σ .For simplicity we will consider the case where k in n and k out n are uncorrelated, so that k out k in = N k 2 .From the definition of the branching function and Eq. ( 1), we obtain Replacing A by the averaged matrix Ā and assuming that s t m is uncorrelated with k out m , we get We note that, since we are replacing A by Ā, we do not expect this estimation to be quantitatively precise.Nevertheless, as we will see, it gives useful predictions on the controllable range and control stability.Expressing the sum as an integral over the in-degree distribution P(k), we obtain In particular, for σ (x) = x for x 1 and σ (x) = 1 for x > 1, we obtain an expression for the branching function (S): Note that if k max k /(λS), then Eq. ( 12) collapses to (S) = λ, showing how the nonlinear case reduces to the linear one, but if k max > k /(λS), then (S) < λ.This means that the more heterogeneous the in-degree distribution is, the more (S) differs from λ.For the case λ = 1, (S) ≈ 1 as long as S < k /k max .Therefore, the range of target values Ŝ that can be controlled with zero error is reduced for heterogeneous networks.
We illustrate how heterogeneity in the degree distribution modifies the controllable range of targets by considering networks of size N = 5000 generated as before but using b = 1 N 1, m = 1 2 1, and { kn } sampled from either a power-law distribution with exponent γ = 4 with minimum degree k min = 50 or a uniform distribution in [100,200], and we choose {k in n }, {k out n } to be two independent permutations of the { kn } sequence.(In light of recent discussion [30], note that we choose power-law degree distributions just as a convenient example of heavy-tailed distributions.) Figure 2(a) shows branching functions (S) calculated using Eq. ( 12) for three combinations of λ and degree distributions (see the caption).Of the two networks with λ = 1, the most homogeneous network has a branching function which is approximately 1 for a wide range of values of S. For the network with λ = 1.05, the branching function crosses 1 only at a value of S approximately equal to 0.53, effectively reducing the controllable range to a single point.All this is reflected in the error R, plotted as a function of the target Ŝ in Fig. 2(b).In all cases, the region where the error is approximately zero corresponds to the region where the branching function is approximately 1.
In addition, the error computed theoretically from the activity S found by solving Eq. ( 8) (symbols) agrees qualitatively with the one obtained numerically from simulations of the full system (solid line), with quantitative agreement except for large values of S. Networks with a homogeneous degree distribution and λ = 1 have the largest controllable range.

IV. STABILITY
Here we discuss how heterogeneity in the degree distribution affects the stability of the controlled fixed point.Considering the case m = μ1 and b = 1 N 1 for simplicity and assuming that the effect of the nonlinear function σ can be neglected [i.e., when (S) ≈ λ], we can study the stability of the controlled system by studying the stability of the fixed point of p t+1 = Āp t + μ1( Ŝ − b p t ).Its stability is determined by the eigenvalues of the Jacobian J = λk in k out /k out k in − μ11 /N.Since the range of J is spanned by k in and 1, to solve for the eigenvalue α in Ju = αu we set u = ak in + b1, with a and b not both zero.We get Grouping coefficients of k in and 1 gives the system which is satisfied for nonzero values of a, b for where h ≡ k 2 / k out k in .For undirected networks, h = k 2 / k 2 1 can be thought of as a measure of the network's degree heterogeneity, with more heterogeneous distributions giving smaller h.In principle, the controlled state becomes unstable when α = ±1.The eigenvalue α = 1 yields instability if μ < (λ − 1)/[1 + λ(h − 1)].This instability corresponds to applying insufficient control, as can be seen in the limit μ → 0, where one obtains α → λ.In this case, network activity S t = b p t will increase (λ > 1) or decrease (λ < 1) until it saturates due to the nonlinearity.A more dramatic instability occurs when too much control is applied and one alternatively overshoots and undershoots the target state, which corresponds to α = −1 and gives the condition From (18), one can see that heterogeneous (smaller h) networks become unstable at smaller control strengths.To illustrate this, in Fig. 3 we plot the rms relative error R 2 t (where • t denotes a time average) obtained from numerical simulations with networks constructed as described above with k out n = k in n = kn and kn sampled from a power-law degree distribution with minimum degree k min = 50 and exponent γ = 3 (blue circles), 4 (red squares), and 8 (green diamonds).We used λ = 1 and λ = 1.2 in Figs.3(a FIG. 3. The rms relative error R 2 t obtained from numerical simulations with undirected networks with a power-law degree distribution P(k) ∝ k −γ with minimum degree k min = 50 and γ = 3 (blue circles), 4 (red squares), and 8 (green diamonds).We used Ŝ = 0.3 and λ = 1 and λ = 1.2 in (a) and (b), respectively.The vertical dashed lines indicate the value of μ predicted from Eq. ( 18), at which the controlled state becomes unstable.state becomes unstable.As μ increases past the value given by Eq. ( 18), which is smaller for the more heterogeneous degree distributions, the error increases sharply.For small values of μ the error also increases as the control is insufficient, and this effect is more pronounced for λ = 1.2, for which there is a strong tendency for network activity to increase.From the cases shown here, the case with the smallest error for the largest range of control values μ is the network with λ = 1 and γ = 8, i.e., the most homogeneous network at criticality.

V. DISCUSSION
In summary, we studied how two key network properties (excitability and degree distribution heterogeneity) affect the ability to control the macroscopic activity of a network of excitable elements with inhibition.We found that for networks poised at the critical point of a phase transition (λ = 1), control was optimal (control error was minimized) for the widest range of control targets.Our results suggest that the cerebral cortex may take advantage of maximum controllability by tuning itself to operate near criticality.One limitation of our work is that we do not answer the question of what mechanisms a neural system may use to tune itself to criticality.Nonetheless, previous experiments [3,5] and theory [31][32][33] suggest that such mechanisms may exist, and we leave this question for further exploration in future work.
A second finding in our work is that heterogeneity in the network's degree distribution reduces the range of target values and the range of control strengths that yield stable control.While heterogeneity can be beneficial for robustness to random node failures [34], our results suggest that a more homogeneous degree distribution might be preferable for situations where control of a large range of macroscopic network activity levels is important.
A common critique of the hypothesis that the cerebral cortex may operate near criticality is that critical dynamics are too noisy, as reflected in the large fluctuations in Fig. 1(a).For many aspects of brain function it is easy to imagine that these large fluctuations would cause trouble.However, our primary result here is that the noisy dynamics of criticality are, in fact, easy to control.This suggests that a brain might be able to take advantage of the other functional benefits of criticality while controlling its own noise to remain at a manageable level.

FIG. 1 .
FIG. 1.(a) Network activity S t as a function of t for three values of λ as indicated, with and without control.Control is turned on at t = 15, 000, with target Ŝ = 0.5 (dashed line), using an Erdös-Rényi network of size N = 5000 and mean degree of 200 with m = 1 2 1 and b = 1 N 1.(b) Absolute control error obtained numerically from direct evolution of Eq. (1) (symbols) and theoretically from Eq. (6) (solid lines) for networks with different correlations (see text).
Figure 1(b) shows the accuracy of Eq. ( . .> kN , with N = 2000 sampled from a distribution uniform in [50, 250].To maximize the term b k in k out m/k out k in , we chose k out n = kn , k in n = kN−n , b ∝ k in , and m ∝ k out ; to minimize that term, we chose k out n = k in n = kn and μ n , b n ∝ kN−n . ) and 3(b), respectively, and Ŝ = 0.3.The vertical dashed lines indicate the value of μ predicted from Eq. (18), at which the controlled