{'fraction_non_alphanumeric': 0.03696902627962691, 'fraction_numerical': 0.013158054748209313, 'mean_word_length': 4.543107516009437, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Applications of optimal nonlinear control to a whole-brain network of FitzHugh-Nagumo oscillators

We apply the framework of optimal nonlinear control to steer the dynamics of a whole-brain network of FitzHugh-Nagumo oscillators. Its nodes correspond to the cortical areas of an atlas-based segmentation of the human cerebral cortex, and the inter-node coupling strengths are derived from Diffusion Tensor Imaging data of the connectome of the human brain. Nodes are coupled using an additive scheme without delays and are driven by background inputs with fixed mean and additive Gaussian noise. Optimal control inputs to nodes are determined by minimizing a cost functional that penalizes the deviations from a desired network dynamic, the control energy, and spatially non-sparse control inputs. Using the strength of the background input and the overall coupling strength as order parameters, the network's state-space decomposes into regions of low and high activity fixed points separated by a high amplitude limit cycle all of which qualitatively correspond to the states of an isolated network node. Along the borders, however, additional limit cycles, asynchronous states and multistability can be observed. Optimal control is applied to several state-switching and network synchronization tasks, and the results are compared to controllability measures from linear control theory for the same connectome. We find that intuitions from the latter about the roles of nodes in steering the network dynamics, which are solely based on connectome features, do not generally carry over to nonlinear systems, as had been previously implied. Instead, the role of nodes under optimal nonlinear control critically depends on the specified task and the system's location in state space. Our results shed new light on the controllability of brain network states and may serve as an inspiration for the design of new paradigms for non-invasive brain stimulation.


I. INTRODUCTION
The widespread use of noninvasive electrical brain stimulation in clinical applications has sparked ongoing interest in studying the effects of external inputs on brain activity. Stimulation with electric fields in the range of 1 − 2 V/m can already modulate oscillatory brain activity [1,2], affect cross-regional synchronization [3,4], and modulate cognitive performance [5]. Clinical studies have successfully demonstrated the efficacy of targeted transcranial stimulation in the treatment of neurological and psychiatric disorders and diseases such as epilepsy [6,7], schizophrenia [8,9], Alzheimer's disease [10] and depression [11].
Brain network models offer a way to simulate and understand the human brain as a nonlinear dynamical system, in which each brain region is represented by a node, and the node dynamics is defined by a model of the average neural activity in that region [12]. Nodes interact with each other according to empirically measured human structural neural connectivity, which quantifies how neural activity in one brain region is coupled to the activity of connected regions.
Parcellation of human brains has yielded various brain atlases [13], which provide information on spatial and functional segregation, dividing the brain into distinct areas [14,15]. The relative connectivity strength between these areas is defined by the number of axonal fibers, which can be estimated using structural neuroimaging scans of individual subjects. This results in a network model of the human brain where the edges reflect the relative strengths of axonal fiber bundles and the nodes represent individual brain regions. These nodes are then equipped with a dynamical system modeling the average neuronal activity in that region [16]. It has been repeatedly shown that when the parameters of brain network models are fitted to functional brain data, optimal operating points were close to the bifurcation lines of these models. This ensures that the model is in a state in which noise fluctuations can be amplified and produce realistic spatial correlation structures which are similar to empirical functional connectivity measurements [17][18][19][20].
Previous theoretical investigations into the impact of external perturbations mostly relied on the assumption of a linear node dynamics, allowing the application of methods from linear control theory [21]. Thus, one can draw conclusions on the effects of external inputs to the system, based on the network topology and independent of its dynamical state [22,23]. On the other hand, linear node dynamics cannot reproduce the dynamics of neural processes close to bifurcations [17,24]. In order to describe the transitions from one dynamical regime to another, Muldoon et al. [25] consider nonlinear node dynamics. They conclude that the effects of stimulationbased control can be predicted by diagnostics from linear control theory based only on the structural connectivity of the network.
In this work we go beyond linear control theory and explore the framework of optimal nonlinear control [26] for the assessing the impact of perturbations on networks of coupled nonlinear systems. Optimal control is an optimization method that derives control policies based on the minimization of a cost functional which depends on the state and control variables. Here, we define a cost functional that penalizes the deviation to the desired network dynamics, the control energy, and the spatial sparsity. The latter allows us to find optimal control signals that apply to only a few control sites, as introduced by [27] and further studied by [28]. We used methods presented by Tröltzsch in [29] for their calculation, while the handling of the stochastic term is based on the work of Stannat et al. [30]. The resulting mathematical formulation is analogous to the formulation presented by Casas et al. [31] for partial differential equations.
We evaluate the framework of optimal control for a brain network of FitzHugh-Nagumo (FHN) [32,33] oscillators with additive coupling and white Gaussian noise, where each oscillator represents a brain region and where the connections between them were chosen according to the human connectome derived from diffusion tensor imaging (DTI) measurements. FHN oscillators are well studied models for neural dynamics and detailed analyses of the dynamical states exist for single oscillators [34] and various network configurations, like two coupled units [35,36], lattices [37], rings and hierarchical architectures [38], as well as random and small-world topologies [39,40] and multiplex networks [41][42][43]. Similar empirical DTI-measured brain connectivities were used with FHN dynamics to model highly synchronized epilepticseizure-like states [44,45], unihemispheric sleep [46], and the functional organization of the resting brain [47][48][49].
We consider two different classes of control problems, targeted attractor switching between multistable network states and increasing network synchronization. The solutions obtained for given energy, precision, and sparseness constraints are well interpretable and result in intuitively sensible optimal control inputs for all network nodes over time.
We then correlate the control energies to controllability measures from linear control theory. We confirm the findings of Muldoon et al. [25] for the investigated state transition (from the low fixed-point state to the oscillatory regime). For other control tasks, however, we show, that diagnostics from linear control theory do not correlate with results from optimal nonlinear control. Conclusions drawn from the structural connectivity alone lead to contradictions, which can only be resolved if the dynamical state of the network and the nonlinear interactions between nodes are taken into account. Applications of nonlinear control theory to whole-brain models enables us to investigate control mechanisms also close to bifurcations. It thus may help in the search for more effective paradigms for realistic transcranial brain stimulation protocols.

A. Network model and control inputs
We consider networks of N equivalent d-dimensional and "noisy" nodes with additive and zero-delay internode coupling. Internode coupling strength is described by a N × N dimensional adjacency matrix A which is scaled by a global coupling strength σ. We allow for an instantaneous and additive control input to the network, which is described by N ×d independent control variables u = (u 1 , . . . , u N ) with u i = (u i1 , . . . , u id ). The equations describing the controlled network dynamics thus read: where ⊗ denotes the Kronecker product. The state of the network is described by a vector x = (x 1 , . . . , x N ), where x i = (x i1 , . . . , x id ) characterize the individual states of the nodes. The local node dynamics is given by h(x) = (h(x 1 ), . . . , h(x N )) with h(x i ) = (h 1 (x i ), . . . , h d (x i )). All nodes additionally receive Gaussian white noise of similar intensity η, which may be correlated within but is uncorrelated across the nodes of the network. The stochastic variables ξ = (ξ 1 , . . . , ξ N ) with ξ i = (ξ i1 , . . . , ξ id ) are independently drawn from a Gaussian distribution, ξ ij ∈ N (0, 1), with zero mean and unit variance. Within node correlations are quantified by the d×d dimensional local noise scheme D, while the across-node statistical independence is assured by a N-dimensional identity matrix I N . The internode coupling term consists of the Kronecker product of the adjacency matrix A and the d × d dimensional local coupling scheme G, where the purpose of the former is to describe the relative interaction strength between nodes while the purpose of the latter is to distribute the between-node interactions across the d variables which describe the local node dynamics. The N × N matrix B in the control term allows for the control of multiple nodes with different strength.
The d × d dimensional matrix K implements the local control scheme, which is similar for every node in the network. Initial conditions are denoted by x(t = 0) = x 0 .

B. The cost functional and the optimality condition
The control u is considered to be optimal (u = u), if it minimizes an appropriate cost functional. To this end, we construct a cost functional F (x(u), u) for a state switching task, where the control input drives the network model from one stable state to another (see Section IV A), and for a node synchronization task, where the control input increases the degree of synchronization among its nodes (see Section IV B). For finite noise, i.e. for finite values of η in Eq. (1), the overall cost functional F is defined as a mean over noise realizations [30], (2) where F n denotes the cost functional for one noise realization n and the angle brackets · denote the mean. In the noise-free case no averaging is performed. The state dependent term F x n (x) only implicitly depends on the control and penalizes the deviation from the desired output. It will be different for the switching, F x n,sw , and synchronization, F x n,syn , tasks (see below). The control dependent term F u (u) accounts for the cost of the control itself.
For the state switching task, we consider a noise-free system (η = 0), and the state dependent cost functional F x n,sw is defined in terms of the deviation of the controlled state to a predefined target state x T (t): where we consider the control being active within the time period 0 ≤ t ≤ T . In order to penalize the precision only towards the end of the controlled time interval rather than during the transition between the initial and target states, the precision-penalizing variable I p can be time-dependent (see Section IV A). The state dependent cost functional for the synchronization task F x n,syn is defined in terms of the deviations of the normalized pairwise cross-correlations R ij for all nodes i and j from R T , the desired mean cross-correlation in the synchronized target state: The cross-correlation in component m is defined as with i, j ∈ [0, . . . , N ] and where . t and σ . t denotes the temporal mean and the standard deviation. The mean is the component-wise network cross-correlation. Later we specialize to m = 1 and suppress the index m.
The input dependent cost functional F u (u) penalizes the energy of the control signal and enforces its directional sparsity [27,28]. It is given by where I e and I s are the weights for the energy and sparsity terms. The first term corresponds to the L2regulation and the second term to the L1-regulation of the cost-functional. Typically, increasing the sparsity I s reduces the number of controlled nodes, i.e. the nodes for which the control input is non-zero for at least part of the time period 0 ≤ t ≤ T , while a higher penalty on the energy term typically leads to an overall reduction of control strengths. Our goal is to find the optimal control u that minimizes the cost functional for chosen I p , I e , and I s , leading to the minimization problem Similarly to Troeltzsch et al. and Casas et al. [29,31,50], we analyze the gradient g = ∇ u F of the cost functional, which has to vanish when evaluated at the optimal control u. By applying the method of Lagrange multipliers, we obtain an expression for the optimality condition that depends on the adjoint states φ(x, u, t) corresponding to the Lagrange multipliers. The control u is optimal, if is the derivative of the sparsity term, Eq. (7), with respect to the control inputs u. The adjoint states are governed by a set of linear differential equations: where D x (h) is the Jacobian matrix of the state equations of the dynamical system, and f x n (x) is the integrand of the cost functional, F x n (x) = T 0 f x n (x)dt. The adjoint state satisfies the boundary condition φ(T ) = 0, and the differential equation is solved backwards in time.
Following Ref. [28], λ(t) is given by Here, k ∈ [1, N ], and E k is the nodewise control energy of the resulting optimal control u, which is defined as resulting in the total control energy E = k E k . A detailed derivation of Eqs. (9) and (10) is provided in Appendix A.

C. Numerical solution of the minimization problem
The optimization problem is numerically solved using the conjugate gradient method. We integrate the equations for the network state and the adjoint given in Eqs. (1) and (10) (see Appendix B for details). The direction d l for each step of the conjugate gradient algorithm is defined by the Polak-Ribiere method [51], while its step size s l is derived using simple bisection [52]. We apply the Fletcher and Reeves algorithm [53] as presented in the following.
We initialize at iteration l = 0 by choosing an initial control u 0 and drawing 20 noise realizations ξ n (t). The corresponding states x 0 (t), Eq. (1), and adjoint states φ 0 (t), Eq. (10), are calculated for every noise realization. Then g(t) as given in Eq. (9) is evaluated. The descent direction is initialized with d 0 (t) = −g(t). We then iterate until convergence: while g(t) ∞ > : 1. compute the step size s l using bisection 3. calculate the states x l+1 (t) and the adjoint states φ l+1 (t) for every noise realization and compute the mean φ l+1 (t) 4. evaluate the gradient g l+1 (t) 5. compute β l using the Polak-Ribiere method: Gradient-based optimization is only guaranteed to converge to local optima of the cost function. In general, however, multiple initial conditions u 0 converged to the same optimum. The only exception occurred for the state-switching task for high values of the sparseness parameter I s , for which a solution with finite values of u coexists with the solution u(t) = 0 (see Appendix B for details).

A. The local node dynamics
We consider a single FitzHugh-Nagumo (FHN) oscillator, with the activity variable x 1 and the linear recovery variable x 2 : where µ is a node-independent, constant background input and R(x) = −αx 3 + βx 2 − γx. The parameters in R are chosen to obtain the bifurcation diagram shown in Fig. 1. This bifurcation diagram, depending on the background input µ, shows three distinct states. In the down FIG. 1: Bifurcation diagram of an uncoupled FitzHugh-Nagumo oscillator, as given in Eq. (13), with the parameters α = 3, β = 4, γ = 1.5, δ = 0.5, and τ = 20. The blue line shows the minimal and maximal values of the activity variable x 1 (identical in fixed points) for the respective background input µ to the node. The red line shows the frequency of the oscillation (time and therefore also frequency are measured in arbitrary units).
state, the node is in a stable fixed point and has a low constant value of the activity variable x 1 (in blue). In the oscillatory state, the activity variable x 1 oscillates at an input-dependent frequency (in red). In the up state, the node is again in a stable fixed point with a high constant value of the activity variable x 1 . The succession of these states in the single node dynamics -a supercritical Andronov-Hopf bifurcation from the down state to the oscillatory state and another supercritical Andronov-Hopf bifurcation from the oscillatory state to the up state [34,54] -closely resembles the states found in large random networks of excitatory and inhibitory spiking neuron models and their corresponding mean-field description [55]. The mean firing rate of these models changes as a function of background input in a way which is qualitatively similar to the value of the activity variable shown in Fig. 1. For this reason one can interpret the value of the activity variable x 1 as the difference between the output firing rate of a cortical node to a baseline value.

B. The brain network model
The topology of the whole-brain network model is derived from diffusion tensor imaging (DTI) data of 12 human subjects (all male and 26-30 years old) from the Human Connectome Project [56] (see Appendix D for subject IDs). The N = 94 nodes in the network correspond to the cortical and subcortical regions defined by the AAL2 atlas-based segmentation [14] (cerebellum excluded). We performed probabilistic fiber tracking (using the Functional Magnetic Resonance Imaging of the Brain Software Library (FSL) [57], details on the processing pipeline in Appendix D) to determine the relative connection strength (edges) between these brain regions (nodes). The pairwise structural connectivity is summarized in the weighted adjacency matrix A shown in Fig. 2. Since DTI data carries no directional information, A must be symmetric. For the following, we define the weighted degree d k of a node k as the sum over all afferent connection strengths, i.e. d k = N i=1 A ik . As motivated in Section III A, the dynamics of each node in the network is described by a FHN oscillator [Eq. (13)]. The general network dynamics, Eq. (1), then simplifies to for all nodes k ∈ [1, N ]. This implies setting the control matrix B = I N , which means that all individual nodes can potentially receive independent control inputs. The local coupling, control, and noise schemes are set to G = K = D = [ [1,0], [0, 0]], i.e. the indvidual FHN nodes receive node-external inputs only through their activity variables x k1 . We do not consider finite delays for reasons of simplicity. Finite delays induce additional dynamical states. Although of high interest in terms of whole-brain modeling, these would not add to the conclusions drawn in Sections IV and IV B about the impact of non-linear optimal control and its comparison with diagnostics [23] derived from linear control and based on connectome properties only. Please refer to Section VI for a discussion of the biological plausibility of this assumption.

C. State space exploration
In this section, we describe and characterize the different dynamical states that can emerge in a whole-brain network of coupled FHN oscillators for a large range of parameter configurations. Having such an overview of the dynamical landscape of the system will allow us to formulate well defined control tasks in Section IV. Throughout the state space exploration we do not apply any control input [u k (t) = 0 ∀ k, t].
Initially we consider the noise free case (η = 0). Here, only two free model parameters remain: the global coupling strength σ and the time independent background input µ which is the same for each FHN oscillator. The state space of the FHN network is explored by simulating the network dynamics for wide ranges of these parameters (see Appendix B for details). The initial conditions x 0 are drawn randomly from a uniform distribution in the interval [0, 1) or taken as the state vector x (t end ) of the last time step of the previous simulation with same σ and slightly smaller µ (continuation). In order to avoid analyzing transient effects, we show and evaluate the network states only after a sufficiently long transient timet.
If oscillations are present, the dynamical state is characterized by the strength of synchronization quantified by the total cross-correlation, Eq. (6), and by the dominant frequency f dom . The latter is given by the frequency of the highest peak of the combined power spectrum of all nodes, where S xx,i1 = |(Fx i1 )(f )| 2 , with Fourier transform F, denotes the power spectral density of an individual node. Figure 3 provides an overview of the state space of the FHN network. Similar to the case of a single FHN oscillator (cf. Fig. 1), the network shows two regions in parameter space, for which a stable fixed point (FP) exists. For small values of µ and σ the activity variable x k1 has a low value for all nodes [down state, cf.  Figure 3a shows that the dynamics can transition from an down state to an oscillatory state, as well as from an oscillatory state to an up state by either increasing µ or σ. For convenience, we call these transitions low and high bifurcation, respectively.
Within the oscillatory regime, we observe network states that are qualitatively different in their appearance, which we explain in the following by their different underlying mechanisms. A single, uncoupled FHN node oscillates in its limit cycle (LC) if it receives a background input in the range of µ ∈ [0.73, 1.33] (cf. Fig. 1). We therefore expect a node in the network to show sustained oscillation (indiv. LC) -and consequently affecting the dynamics of the network -if its combined input is in this interval. We thus call a node "to be in its individual LC regime" if for at least one point in time within the considered time interval. This criterion is further motivated by observations of the dynamics (cf. Fig. S1 in [78]) which show qualitative changes in behavior depending on whether it is met by all, some, or none of the nodes in the network. If all nodes fulfill Eq. (16), we observe the network to be in a synchronous, high amplitude limit cycle [haLC in Fig. 3a, cf. Fig. 3c(5)]. If, however, no node fulfills Eq. (16), the network either inherits the fixed point solution from the individual nodes, or a new state emerges due to network effects. We observe both cases: The first one is true in both down and up state of the network in Fig. 3a. The second case appears for low and intermediate coupling strength σ at the transition from down or up state to oscillatory state. In these states the FP of the network dynamics is unstable due to coupling ef-fects and a self-sustained low amplitude limit cycle [laLC in Fig. 3a, cf. Fig. 3c(2,7)] emerges. The geometric interpretation of this laLC is a rotation around the FP that would be transient for an isolated node, but is prevented from converging to a fixed-point by the network couplings.
It is also possible that only a fraction of the network nodes fulfill the criterion in Eq. (16), which can leadeven in the absence of noise and without network delays -to asynchronous and apparent aperiodic behavior [mixLC in Fig. 3a, cf. Fig. 3c (3,6), further indications for aperiodicity are provided in [78]]. On the network level, this can be seen as the result of the different frequencies of the two coexisting network limit cycles, laLC and haLC, interacting. Close to the bifurcation the frequency of the individual nodes (cf. Fig. 1), and their trajectories (cf. Fig. S3 in [78]), are particularly sensitive to their respective input. If the additive coupling input between the nodes is not sufficient to entrain a common frequency, this may result in asynchronous and potentially aperiodic network oscillations. States that are classified as mixLC are, however, not necessarily asynchronous. If the driving force of nodes that fulfill the criterion in Eq. (16) is high enough, it will lead to frequency entrainment. Therefore, the oscillations of the driven nodes [where Eq. (16) is not fulfilled] are similar to the ones of the driver nodes, resulting in dynamics that are similar to the haLC state [cf. Fig. 3c (4,5)].
The distinction between different dynamical states and types of network oscillations (laLC, mixLC, haLC) is particularly interesting since the transitions between these are the regions in state space, where we find multistable network states. Multistability is detected by simulating the network dynamics for 21 different initial conditions x 0 and comparing the resulting time series by calculating their node-wise correlation in the activity variable. We ignore the firstt = 5000 arb. units of each time series to avoid analyzing transient effects and then compare the interval t ∈ [5000, 6000] (equiv. to 10-35 periods) of one initial condition with t ∈ [5000, 6200] of all other traces (second interval is sufficiently longer to account for all phase shifts). If the auto-correlation of at least one initial condition is close to one, this indicates the existence of a stable (periodic) state, and if the cross-correlation between two (or more) initial conditions is different from one, this indicates bi-(or multi-)stability. Red pixels in Fig. 3a show states with more than one stable solution, which are observed along both, the low and high bifurcation.
Numerically we found regions of bistability between FP and mixLC (A1), laLC and mixLC (A2, B2), and mixLC and haLC (B1), as well as multistability in mixLC with different numbers of nodes fulfilling the criterion in Eq. (16) (A3, B3) (cf. Fig. 3a). The corresponding time series of the activity variables x k1 are shown in Fig. 4. The automatic detection of multistability might not capture all multistable states, especially because some of the states may be very similar (cf. Fig. 4b, bottom panel). More detailed analyses could thus also uncover bistabilities between up state and mixLC, as well as between mixLC and haLC at the low bifurcation. We inspected the detected multistable states visually to assure that their difference are not due to transient effects and confirm their assignment to the specified multistable regions in Fig. 3a. We do not observe state switches in this noisefree case even for very long simulation times (200, 000 arb. units simulated for multiple initial conditions in each multistable region, A1-3, B1-3, of Fig. 3a), meaning that the states are stable and no temporal intermittency is found.
Adding noise to the network most strongly affects the dynamics of the states close to the bifurcations. The influence of noise on synchronous, mono-stable oscillatory states is, as expected, comparably small (cf. Figs. S4d, e in [78]). If the network dynamics is in a FP far away from the bifurcation line, the additional Gaussian white noise results in uncorrelated fluctuations in the activity variable (cf. Fig. S4i in [78]). Parametrizations that without noise would lead to a stable FP close to the bifurcation, show oscillations when sufficient noise is introduced (cf. Figs. S4a, h in [78]). The clear distinction between the FP and laLC network states in the noise-free case (cf. Fig. 3a) is therefore blurred for noisy dynamics. For asynchronous mixLC states, the additional noise can have a synchronizing effect (cf. Fig. S4c in [78]). This can be explained by more nodes fulfilling the criterion in Eq. (16), resulting in a stronger drive for the remaining nodes and therefore more synchronous oscillation. These effects lead to a shift of the bifurcation lines compared to the noise-free case and consequently also affects the location of multistable network states.
In the case of multistability, we observe that this stochastic additional input can lead to sufficient perturbations that drive the dynamics from one attractor to the other. This results in noise-induced state switching, as shown in Fig. 4c, at the bifurcation lines of the noisy state space. Such noise-induced transitions are a known phenomenon in systems of coupled oscillators [58] and may be exploited for the control of the neural dynamics in practice [59].

A. Switching between bistable network states
In this section, we present optimal control inputs that induce a switch between previously identified multistable states. All control inputs optimize the cost functional given in Eq. (2) with the state dependent terms given in Eq. (3) and the energy and sparsity terms given in Eq. (7). Energy and sparsity terms are evaluated over the whole time interval during which the control is active. The cost functional itself, however, considers the deviation from the target state only the end of the control period. For this we set I p (t) = I * p for (T −τ ) ≤ t ≤ T and I p = 0 else. The convergence criterion of the numerical solution of the minimization problem, as described in Section II C, is set to = 10 −5 for all our applications. We ensured for all presented examples that the method actually converges and that results do not change for lower values for . We computed the optimal control for different phase shifts of the initial and target states with respect to the control onset and present the results for the phase shift which leads to the smallest total control energy E. Figure 5 shows the result of applying optimal control at point A1 in Fig. 3a, where a low activity fixed point coexists with an oscillatory mixed state. The weight I s for the sparseness term in Eq. (7) was set to zero. Since the target state in this task is always stable, the network remains in this state once it is reached also after the control is turned off. If the initial state is the down state FP, the obtained optimal control input oscillates synchronously for all nodes with increasing amplitude. The frequency of the control input [f dom (u 1 ) = 35.0] corresponds to the frequency of the fixed point's focus [f dom (x F P 1 ) = 35.0, measured by perturbing the FP], inducing resonance effects in the node activity. This strategy can be well observed in Video 1 (in [78]), showing the node oscillations and the respective optimal control inputs in state space.
When switching from the oscillatory to the down-state, as in Fig. 5d, the optimal control input consists of one short biphasic pulse applied to all of the nodes followed by minor corrections. It is a known result from control theory, that applying a biphasic control pulse around an extreme point is an efficient way to drive an oscillating system to a stop. Interestingly, the optimal control strategy in this example is not to apply the pulse at an extreme point of the activity variables x k1 , but rather in a way that the gradient of the control is highest when the phase velocity of the limit cycle oscillation is the lowest (best observed in Video 2 in [78]). Although the deviation from the target state is only penalized at the end of the control interval, the pulse is applied early in order to give the system enough time to come to rest. Figure 6 shows the result of applying optimal control at point B1 in Fig. 3a, where a low amplitude limit cycle (laLC) around the high activity up-state coexists with an oscillatory mixed state (mixLC, cf. Section III).  [78]), the optimal control finds a short, off-phase biphasic pulse analogous to Fig. 5d. As a result, the amplitude of the network oscillation decreases to almost zero (for t ∈ [200, 250] in Fig. 6b). Since there is no stable FP, the oscillation amplitude increases again over time and converges to the laLC target state approx. 550 time units after the control pulse (see also Fig. S5 in [78]).
By increasing the sparsity parameter I s of the cost functional, we can tune the number of controlled nodes. Figure 7 shows the control energy E k [Eq. (12)] of the optimal control input u k1 to each node k as a function of the sparsity parameter I s . For low values of I s , all nodes receive a finite control signal. When I s is increased, less nodes are controlled, until I s becomes so large, that the target state can no longer be reached. As expected, decreasing the number of controlled nodes needs to be compensated for with a higher control energy.
The results also show that at the low bifurcation (Figs. 7a, c), nodes with a high weighted degree receive, independent of the switching direction, a stronger control signal and remain controlled even for high values of I s . This shows that at the low bifurcation it is most important to control the network hubs since they act as This different behavior can be explained based on the additive coupling scheme between the nodes (cf. Section III), which causes nodes with higher degree to typically receive stronger inputs. When choosing parameters close to the low bifurcation, these high degree nodes therefore have higher oscillation amplitudes (cf. Fig. 1) and are consequently driving the oscillation of the remaining nodes in the network. Close to the high bifurcation the nodes with lower degree, who receive less inputs, are more likely to remain in the limit cycle regime. Consequently -and in contrast to the dynamics close to the low bifurcation -the low degree nodes drive the oscillation of the network hubs. Videos 1 and 3 (in [78]) illustrate how the control inputs on the respective driver nodes force the network dynamics to the high-amplitude oscillation target states.
The optimal control inputs to the control sites have well interpretable shapes. When the objective is to transition from a low-amplitude state to one with a higher amplitude (mixLC, Figs. 5a and 6a), the optimal control strategy utilizes the characteristics of the flow field in state space and synchronously drives the network with its resonant frequency from one attractor towards the other. Switching in the opposite direction, and therefore leaving this basin of attraction, is achieved most efficiently with one strong, biphasic pulse. This deflects the system from its initial mixLC trajectory and causes it to drop to the attractor with a lower amplitude.

B. Synchronizing the network dynamics
In this section, we apply the nonlinear control method to find the optimal inputs to synchronize the dynamics of the individual nodes. For this application we use the state dependent cost functionals given in Eq. (4) to penalize the deviation from the target cross correlation (R T = 1, fully synchronous state) and in Eq. (7) to penalize the energy of the control and enforce its sparsity in space. We parameterize the system to be in the two asynchronous regions marked by the labels S1 and S2 in Fig. 3a, close to the low and high bifurcation lines. The optimal control method is then applied to synchronize the dynamics for the noise-free (Fig. 8) and noisy (Fig. 10) case. Both figures show the controlled (synchronous) and uncontrolled (asynchronous) time series and the corresponding optimal control inputs. Since the dynamics in points S1 and S2 are monostable, the synchronous state is not a stable solution and the system returns to its original state as soon as the control is switched off. The network dynamics and the optimal control strategies are best seen in Videos 5-8 (in [78]). Figure 8 shows how the optimal control inputs synchronize the oscillation of all nodes in the network. In the uncontrolled state (in gray, Figs. 8a, b) the nodes at both bifurcations oscillate with similar frequencies (mean and standard deviation of oscillation frequencies across nodes at S1: f N = 32.2, σ f N = 1.1, and at S2: f N = 27.0, σ f N = 0.9). Close inspection of the control inputs in Figs. 8c and d shows that the optimal control acts in two phases. First, the control aligns the phases of all oscillators until all nodes are synchronized for the first time, which we defined as the critical time t c (please refer to [78] for details about the computation of t c ). Second, a periodic input maintains the synchronous oscillation during the period the control is active.
For times t < t c , Fig. 9a shows that at the low bifurcation the energies E k of the control inputs to the nodes are positively correlated with the weighted node degrees. Thus, it is beneficial to focus the control input on the network hubs for alignment. A possible interpretation for this strategy is that the optimal control utilizes the influence of the network hubs on the remaining nodes to force them on the synchronous limit cycle trajectory. However, we again observe a different behavior at the high bifurcation, as shown in Fig. 9b. Here, E k is negatively correlated with the weighted node degree for times t < t c . In this case, the network coupling is much smaller compared to the background input µ and the node's phase space oscillations are similar to the trajectory of the limit cycle of an uncoupled node (cf. Video 6 in [78]). The negative correlation in Fig. 9b suggests that if the dynamics of the nodes are similar in the first place, it is beneficial to focus the control on more weakly coupled nodes and align their phase space trajectory with the trajectory of the network hubs. Similar to the attractor switching at the high bifurcation (cf. Figs. 7b, d), the control of the low degree nodes drives the network oscillation towards the target state.
In the second phase, for t ≥ t c , the high degree nodes without control would have a higher (or lower) phase velocity than the nodes with intermediate degree, while the opposite is true for nodes with low degree. The optimal control strategy is thus to act on the low and high degree nodes in opposite directions, which can be observed as antiphase control inputs shown in Figs. 8c and d, where nodes with intermediate degree receive only small control inputs. (Videos 5 and 6 in [78] show how the control inputs keep the nodes together in phase space.) This is reflected in the arch-like form of the mean energies E k of the control inputs at both bifurcations, as seen in Figs. 9c and d, with high control energies for nodes with low or high degree. As a result, the optimal control achieves a mean network cross-correlation of R = 0.995 in the time interval t ∈ [t c , T ], compared to the uncontrolled scenario with R = 0.289 at point S1 (Fig. 8a), and R = 0.996 compared to R = 0.368 at point S2 (Fig. 8b). Figure 10 shows results for the application of optimal  Fig. 3a. All parameters are as in Fig. 8. control to the synchronization task for the case of finite noise. We obtain a mean network cross-correlation of R = 0.83 (analysis of the deviation from R T in [78], uncontrolled: R = 0.25) at point S1 (Fig. 10a) and R = 0.85 (uncontrolled: R = 0.35) at point S2 (Fig. 10b). We conclude that the optimal control input successfully synchronizes the network dynamics, just as in the noise-free case.
The control strategy that leads to the respective target state, however, changes substantially when noise is added to the system. The control input in the noise-free case is tailored individually for each node, slowing the phase space velocity of some nodes down and speeding others up at the same time (cf. blue inset in Fig. 8d, mean cross-correlation of the control inputs during this control period is 0.03). In contrast to that, the control input for the noisy network dynamics has a similar shape for all nodes (cf. blue inset in Fig. 10d, mean cross-correlation of the control inputs during control period is 0.88) with varying amplitude depending on the weighted node degree d k .
This change in control strategy is also reflected in Fig. 11, which shows the control energy E k per node for different noise levels. The control energy increases for all nodes with increasing noise strength, both at the low (Fig. 11a) and the high bifurcation (Fig. 11b). With increasing noise level η, we observe a gradual transition of the optimal control strategy. Instead of balancing out  12)] as a function of the weighted node degree d k for different noise levels η for the synchronization task. The error bars show the standard deviation over 10 different initial conditions of the network dynamics with 20 independent noise realizations of the optimization each. E k is, in contrast to Fig. 9, calculated over the whole interval during which the control is active. (a) Point S1 (low bifurcation) in Fig. 3a. (b) Point S2 (high bifurcation) in Fig. 3a. Parameters were I p = 0.1, I e = 1.0, I s = 0, all other parameters are as in Fig. 8. the phase velocity differences between nodes towards the phase of nodes with intermediate degree (cf. Videos 5 and 7 in [78]), the control in the noisy case forces the system on a limit cycle with a slightly higher oscillation amplitude (cf. Videos 6 and 8 in [78]). This strategy requires a higher control energy but is independent of the specific realization of the noise.
By imposing sparsity constraints on the system, we can investigate to what extent synchronization can be achieved with fewer control sites. Figure 12 shows the relation of the synchrony in the network, measured by the average cross-correlation R [Eq. (6)], to the sparsity parameter I s . The average cross-correlation R decreases with increasing sparsity parameter I s at both locations in state space and for all noise levels. This shows that, under the imposed constraints, it is not possible to fully synchronize the network with sparse control. Optimally synchronizing aperiodic states is a collective effort. Unlike state switching it cannot be achieved by controlling only a few sites, because the reduced number of control sites cannot sufficiently be compensated for by a higher control energy. Especially in the case of noisy network dynamics, where the control must drive all nodes, the network's synchrony quickly drops to the value of the uncontrolled network when I s is increased. For fixed I s the number of controlled nodes and the control sites optimal for synchronizing the network changes for each initial condition x 0 (cf. Fig. S1 in [78]), as it depends on details of the exact network dynamics.

V. COMPARISON WITH CONTROLLABILITY MEASURES DERIVED FROM LINEAR CONTROL THEORY
When relating functional properties of a neural system to properties of the underlying connectome, neural activity is often approximated by linear dynamical systems [23,60,61]. This has obvious benefits for analyzing the effects of perturbations. Calculating optimal control inputs for linear systems, as in Refs. [12,22], can be done analytically and with little computational effort, and conclusions about the effects of external inputs can be drawn from controllability measures, which depend on network topology only [21].
Two of these measures have previously been applied to quantify the impact of perturbations in a whole-brain network setting [22,23,25]. Modal controllability refers to the ability of a node to control each evolutionary mode of a dynamical network [62], and the average controllability is given by the average control input energy to the respective node over all possible target states [23,25] (mathematical description in Appendix C). Nodes with high average controllability require only low energy input to move a linear system into "easy-to-reach" states, while nodes with high modal controllability require a high control energy input to have an effect on the dynamics but are crucial when the target state of the system is "difficult-to-reach" [21,23,63]. It has also been shown that the average controllability of a node is strongly correlated with its weighted degree d k (r = 0.85, p < 10 −26 for our structural connectivity matrix) while the modal controllability is known to have a strong inverse correlation (r = −0.82, p < 10 −23 , see Fig. S7 in [78]).
We now apply the diagnostics from linear control theory to the brain network model, Eq. (14). Figure 13 shows the correlations of the energies E k of the optimal control inputs with the average and modal controllability for the attractor switching tasks. At the low bifurcation A1 (Figs. 13a, c) we find, that the energy of the control input for a node is positively correlated with its average controllability and negatively correlated with its modal controllability, both irrespective of the switching direction. Therefore, the most efficient strategy for switching between the attractors is to control the network hubs with high average controllability which then force the other nodes towards the target state. These results can be considered consistent with predictions from linear control theory [23] and previous results on the global impact of stimulation [25].
At the high bifurcation at location B2 (Figs. 13b, d), however, we observe the opposite trend, with the control energy being negatively correlated with the average and positively correlated with the modal controllability. Predictions based on linear control theory alone are indifferent to the actual dynamics and therefore not able to distinguish between the two cases. Yet, nodes with high modal controllability play a much more important role in affecting the global dynamics at this location in state space.
For the synchronization task, we observe a similar pattern. While there is a clear correlation between the optimal node-wise control energies E k with the average controllability when aligning the phases (i.e. t < t c ) at the low bifurcation (Fig. 14a), no correlation is observed at the high bifurcation (Fig. 14b). (The corresponding results involving modal controllability are shown in Fig. S8 in [78].) In the case of maintaining the synchronous network dynamics (t ≥ t c , Figs. 14c, d), we observe -at both bifurcations -a similar but less obvious arclike shape as in Figs. 9c, d. For noisy network dynamics (cf. Fig. S9 in [78]), however, E k is negatively correlated with the average controllability at both bifurcations (cf. in [78]). This shows that in our nonlinear setting, linear controllability measures do not provide additional insights compared to the weighted node degree, and that intuitions based on these measures can be misleading. Diagnostics from linear control theory have previously been claimed to be predictive for a node's role in driving brain state transitions [22,23] or for the global impact of local brain stimulation [25]. Our results however show, that the optimal control inputs and sites in nonlinear systems not only depend on the structural network connectivity, but also on the location in state space, the control task, and other factors like the amount of noise.

VI. DISCUSSION
In this contribution we apply techniques from the optimal control of nonlinear dynamical systems to the dynamics of brain network models. Nodes were equipped with FitzHugh-Nagumo oscillators, since they are simple and well studied nonlinear models for neural dynamics.
FIG. 14: Mean energies E k of the optimal control inputs [Eq. (12)] plotted against the average controllability for the synchronization task (cf. Fig. 9). (a) E k vs. average controllability, with control signal considered during the initial time interval t ∈ [0, t c ), at the low bifurcation (Point S1 in Fig. 3a). (b) Same as (a) but at the high bifurcation (Point S2 in Fig. 3a). Solid red lines denote the linear regression results with the following obtained coefficients: (a) r = 0.87, p < 10 −29 , (b) r = −0.16, p = 0.12 (not significant). (c) E k vs. average controllability, with control signal considered for time t ∈ [t c , T ], at S1. (d) Same as (c) but at S2. All parameters as in Fig. 9.
Changing the background input for the FHN nodes or the global coupling strength of the network can both lead to transitions between two different stable fixed points and two different limit cycle attractors (laLC and haLC). The interaction between nodes in different oscillation states (mixLC) can lead to an asynchronous network dynamics. At the bifurcations, we also find different coexisting stable states.
The general mathematical framework of the optimal control of partial differential equations [29] is adapted for noisy dynamical systems on graphs, where the local network dynamics, the noise level, the network connectivity, and the local coupling schemes [Eq. (1)] can be freely chosen. The state dependent part of the cost functional is averaged over multiple noise realizations and penalizes the deviation of the network dynamics from a task dependent control target [Eq. (3) for attractor switching, Eq. (4) for synchronization]. The part of the cost functional which depends on the control input [Eq. (7)] penalizes the control energy and non-sparse solutions. The presented method is applicable to any network that can be described in the form of Eq. (1), including models of power grids, social networks, or climate dynamics.
A common problem for gradient based methods are potential local optima of the cost functional which may prevent the convergence to a globally optimal solution. To alleviate this problem we performed the minimization as described in Section II C with different initial condi-tions u 0 for the control time-series and choose the result with minimal cost. The set of initial conditions included u 0 = 0 as well as valid control time-series taken from different parametrizations.
We use the optimal control method to cause targeted attractor switching between previously identified coexisting stable states. When no sparsity is enforced, we show that it is optimal to resonantly drive all nodes to transition to the high amplitude oscillatory state. When the task is to switch to a state with lower amplitude or no oscillation, the optimal strategy is to apply a precisely timed biphasic pulse. The nodes that receive the largest control energy are the same for both of these switching directions and are also the ones that are still controlled when we enforce sparsity in space. Depending on the location in state space, either nodes with high degree (at the low bifurcation) or low degree (at the high bifurcation) are the ones that most efficiently drive the network dynamics from the initial to the target state. When sparsity is enforced, controlling only a small number of these driving nodes with increased control energies E k is sufficient to switch from one attractor to another in an optimal way.
In the second application, we show that our method can also be used to control global properties of the network dynamics, such as the average cross-correlation between node activities in the oscillating regime. Immediately with control onset, the control signal acts on all nodes to align their phases. As soon as this is achieved, the control maintains the synchronous oscillation with periodic control signals. Which nodes receive larger control input for the initial alignment again depends on the location in the state space. Both at the low and high bifurcation, individually adapted control inputs lead to a successful synchronization of the network dynamics. While the average cross-correlation R is increased also with sparse control, the synchronous target state is only achieved in an optimal way, when all nodes in the network receive a finite control input. This suggests that synchronizing all nodes needs collective intervention, while attractor switching can be caused by controlling a few selected nodes only. The introduction of noise to the system makes the dynamics of each node less predictable, resulting in a loss in specificity and more similar control inputs to the nodes. Consequently, the optimal control of noisy network dynamics requires higher total control energy E (cf. Fig. 11) while resulting in a lower precision (cf. Fig. 12 for I s = 0).
The information on the different states and bifurcations, which we show to be a decisive factor for choosing the optimal control sites in our applications, is lost when techniques from linear control theory are applied. While predictions do qualitatively agree for certain dynamical systems, control tasks, and locations in state space (cf. [25]), this agreement does not hold in a general setting. It would, however, be worthwhile investigating for what classes of dynamical systems and control tasks "controllability measures" can be defined, which only depend on the properties of the connectome.
In this contribution, the techniques of nonlinear optimal control were applied to a simplified model of the global brain dynamics. The bifurcations in our FHN model (cf. Fig. 1) phenomenologically capture the state transitions found in more complex, biophysically motivated network models [55] as we have adapted the FHN parameters to resemble their dynamics. Consequently, the activity variables of the FHN nodes can be interpreted as relative output firing rates of cortical nodes, and their total input µ + σ N i=1 A ki x i1 (t) can -at least qualitatively -be interpreted as a proxy of the local field potential [64]. Given this interpretation, model down-(up-)states are related to down-(up-)states in cortical physiology, while an oscillatory regime of the wholebrain model corresponds to a brain state with oscillations present. Locations at the lower bifurcation of the network model were implicated as proper "operating points" for modeling the brain's resting state activity [17][18][19], while locations at the high bifurcation can be considered a cartoon model for the global brain dynamics during non-Rapid Eye Movement (non-REM) sleep (see [20] for a biophysically more detailed model).
In this work we consider a network with instantaneous couplings, although coupling delays in human white matter are finite. The interpretation of results regarding the model's state space should thus be limited to cortical states whose dynamics is slow compared to these delays. Given that delays are of the order of approx. 5-15 ms [65] brain oscillations with periods of approx. 1 s or less would qualify. These include the so-called slow oscillations, which are the prominent brain rhythm during non-REM sleep and which have already been subject to external perturbation experiments in human neuroscience and clinical settings [66,67].
In lieue of a biologically more detailed model of the effects of control inputs on whole-brain dynamics, delays must be included as soon as brain oscillations with frequencies above a few Hz are of interest. Relative delays can be computed from DTI-based estimates of the length of white matter fiber tracts and included in the coupling terms of Eqs. (1) and (14).
When relative delays are scaled by a global delay constant, the ratio of the coupling delays to the FHN oscillation period can be adapted to match the physiologically observed ratios for any given brain rhythm of interest. The formalism of nonlinear optimal control summarized in Section II must extended, though, to cover finite delays by modifying the coupling term in Eq. (1) to, e.g., σ(A ⊗ G)x(t − d) for the case of constant delays. The adapted coupling term reappears when computing the adjoint state with the differential equation [Eq. (10)], and thus the iterative algorithm for the numerical optimization (cf. Section II C) has to be adapted accordingly in the calculation of step 3.
Here we used a highly simplified model of the global brain dynamics to showcase the applicability of nonlinear optimal control and its added value beyond connectome-based diagnostics derived from linear control theory. When applied to biophysically more realistic models of whole-brain activity (cf. [20]), nonlinear optimal control may serve as a tool to evaluate the impact of external stimulation and to facilitate the design of new brain stimulation protocols. Non-invasive brain stimulation like transcranial current stimulation [68][69][70] is a highly promising technique to perturb the global brain dynamics with the goal to improve sensory [71], motor [72], and cognitive abilities [66] of human subjects. Using a cost functional which penalizes spatial sparseness and control energy may -for example -help reducing the required current applied to a subject's brain as well as focusing the electrical perturbation on the relevant brain areas only. While exact target trajectories are unlikely to be available in a practical setting, state-dependent cost functionals which refer to global quantities [cf. Eq. (4) for the synchronization task] may well be. Particularly, reformulating the control formalism in frequency space will be beneficial here. It would for example allow for computing optimal interventions for changing the power of certain brain rhythms, which can be monitored by electroencephalography and which are common control targets in clinical settings (cf. [67]). As shown in other computational and physiological studies [55,73] and emphasized again by our results, the timing of the control inputs may also be crucial for successful interventions.
Here we see a high potential of nonlinear optimal control in guiding electrical brain stimulation under electro-or magnetoencephalography [74,75], a setting which allows for such precisely timed control inputs.
For quantitative predictions, an application of nonlinear optimal control to biophysically grounded network models (cf. [20]) is desirable, where the coupling of neurons to externally applied electric fields is included in a biophysically realistic way. Still, results obtained with simplified brain-network models will strongly facilitate biologically detailed but computationally expensive insilico experiments and may already inspire physiological studies on brain-stimulation to explore new control paradigms for clinical stimulation protocols with potentially better control results, reduced energy expenditure, and higher spatial sparseness.

VII. ACKNOWLEDGMENTS
We thank Dr. Michael Scholz, Cristiana Dimulescu, and Lena Salfenmoser for scientific discussions, and Prof. Dr. Eckehard Schöll for his comments on the manuscript. In this section the mathematical foundations for controlling the network dynamics are laid. We use nonlinear sparse optimal control. All derivations below are similar to the ones in Refs. [29,31], where sparse optimal control was applied to a system of partial differential equations with diffusive coupling.
We consider a network of N nodes with d-dimensional dynamics each. Its dynamics is defined by a system of N stochastic differential equations, with the state vector x = (x 1 , ..., x N ) and x i = (x i1 , . . . , x id ). ⊗ denotes the Kronecker product. The local node dynamics is governed by The coupling term consists of the N × N -dimensional adjacency matrix A and the d×d-dimensional local coupling scheme G. u = (u 1 , ..., u N ) with u i = (u i1 , . . . , u id ) denotes the vector of control inputs, B is a diagonal N × N dimensional control matrix, and K the d×d-dimensional local control scheme. The noise term consists of the Kronecker product of an N-dimensional identity matrix A and the d × d-dimensional local noise scheme G. The independent white Gaussian stochastic process is given by ξ with standard deviation η. Boundary conditions are given by (A2) Since our system is stochastic while the control is deterministic, the optimal control must minimize a cost functional which is an expectation over many realizations of the noise. It is defined as where the angle brackets denote the expectation. Our goal is to find the optimal control u that minimizes the above mean cost functional. Hence, the variational inequality holds when regarding every dimension i separately. We linearize around the optimal solution and write δu i = Since above inequality (A5) has to hold for all δu i and −δu i , we find the stronger necessary condition for the optimal solution. Inserting the expression from the right side of Eq. (A3) yields Since Eq. (A7) must hold for all components i, we write where • denotes the element-wise multiplication (Schur product). Here we assume that the cost functional can be separated into a part F u (u), which depends on control inputs only and which is deterministic, and a part F x (x(u)) , which is a function of the network state x and only indirectly depends on the applied control, i.e.: Equivalently, one can write Applying the chain rule in Eq. (A8) yields where D u (x(u)) is the Jacobian matrix of the state x with respect to the control vector u, evaluated at the optimal control u.
To calculate the above derivative we must find an expression for the Jacobian matrix D u (x(u)) . However, the dependence of the network state on the control inputs cannot be provided in closed form. The Lagrange Formalism enables us to derive an alternative expression for the second term on the right-hand side of Eq. (A11), which can be evaluated to find the optimal control.

The method of Lagrange multipliers
We can use the formalism of Lagrange multipliers to formulate conditions for the optimization of the cost functional, Eq. (A3), that is subject to the constraint given in Eq. (A1) and the boundary conditions given in Eq. (A2). Our approach is similar to the method used in [68], and we will obtain a result that resembles what is found in [31] for a system reaction-diffusion equations.
We define the mean Lagrange function L(x(t), u(t)) as Equation (A1) for the dynamics of the network state x has to be satisfied for all times 0 ≤ t ≤ T , which is ensured by the time integration of the constraint. φ(t) = (φ 1 (t), ..., φ N (t)) with φ i (t) = (φ i1 (t), . . . , φ id (t)) is the d × N -dimensional vector of time-dependent Lagrange multipliers.
After inserting the Lagrange function, Eq. (A12), into Eq. (A14), we apply partial integration and the chain rule, and obtain The boundary term must hold for every realization. Since δx(T ) could be finite, the boundary condition φ(T ) for the vector of Lagrange multipliers is Equation (A15) is fulfilled if the integrand vanishes. We thus obtain a d × N -dimensional system of linear ordinary differential equations for the so-called adjoint states φ, where D x (h) denotes the Jacobian matrix of the local dynamics with respect to the network state x. The adjoint states are obtained for every realization by solving Eq. (A18) backwards in time obeying the boundary conditions (A17).
Inserting the Lagrange function [Eq. (A12)] into Eq. (A13), we obtain The first term in Eq. (A19) vanishes [see Eq. (A8)]. Furthermore, after applying partial integration and the chain rule, above Eq. (A19) simplifies to where the boundary term of the partial integration was set to zero because of the boundary conditions imposed on φ and x.
Inserting the differential equations (A18) that govern the dynamics of the adjoint states into Eq. (A20), we obtain The first term corresponds to the expression found in Eq. (A11). We can solve Eq. (A21) for this term and insert it into Eq. (A11) to obtain the optimality condition Equation (A22) is fulfilled if the stronger optimality condition When solving the minimization problem for the optimal control numerically (see Section II C) using gradient based optimization, we may converge to a local minimum of the cost functional. In order to alleviate this problem, we choose multiple initial conditions u 0 for the conjugate gradient algorithm. In the case of synchronizing the network dynamics, the control scheme was robust to changes in u 0 . The coexistence of multiple local minima of the cost functional was observed only in the case of switching between network states using sparse control (cf. Fig. 7). For larger values for I s [cf. Eq. (7)] we find one minimum at u(t) = 0 (smaller cost for the sparsity term) coexisting with another minimum corresponding to a control input u(t) affecting a small number of nodes (smaller cost for the precision term). Which of these local optima has minimal cost depends on the value of I s .
To ensure that the global optimum was found in the latter case for a given value of I s , we systematically varied the initial conditions. Starting with I s = 0, we first calculated the optimal control signal, then increased the value of I s , and used the previous optimal control as the new initial condition (continuation). Likewise, starting with a high enough value for I s , this process is repeated, but for decreasing values of I s . For a given value of I s the control input with the lower cost was chosen for further analysis.

Appendix C: Definitions of average and modal controllability
We consider a simplified linear network dynamics of the form with activity vector x, adjacency matrix A, control matrix B and control inputs u. In order to calculate the average controllability of node k in the network, we choose the input matrix B (k) to target a single node, i.e. B (k) ij = 1 for i, j = k and B (k) ij = 0, otherwise. The controllability Gramian of this system is defined as and the average controllability as its trace (cf. [23]) Network nodes with high average controllability have a large impact on the network dynamics. This makes them important control sites since by controlling these nodes a large number of possible activity vectors x can be reached with little control energy. They are consequently said to be able to move the system to many easy-to-reach states [23].
The modal controllability is defined based on on an eigenvalue decomposition of the network adjacency matrix A, where v kj are the elements of the eigenvector matrix V = [v kj ] and λ 1 . . . λ N are the corresponding eigenvalues [23]. Nodes with high modal controllability tend to be sparsely connected. They are capable of pushing the network dynamics towards activity vectors x which require substantial control energy to be reached. Thus, nodes with high modal controllability are said to be capable of steering the network into difficult-to-reach states [23]. For a detailed description and derivation of these regional controllability measures please refer to Refs. [23,25,63]. We use the matlab code provided with the publication of Ref. [23] to compute the controllability measures [79]. For their calculation, only the adjaciency matrix A (cf. Fig. 2) is required.

Data processing pipeline
Data acquisition and preprocessing is described in Ref. [56]. We used BEDPOSTX (Bayesian Estimation of Diffusion Parameters Obtained using Sampling Techniques) from the FSL toolbox [57] for building up distributions on diffusion parameters, which automatically determines the number of crossing fibres per voxel. Based on these diffusion parameters at each voxel, we then applied the PROBTRACKX (probabilistic tractography) algorithm from the same toolbox with 5000 samples per voxel. The resulting connectivity matrix A s for each subject s was normalized by dividing the connections between any two regions by the number of voxels in the source region multiplied by 5000. Self connections were deleted, A s ii ! = 0, ∀ i, and we averaged the structural connectivity matrices of the individual subjects to obtain the adjacency matrix A = 1 12 N =12 s A s . Probabilistic fiber tracking does not provide information about directionality; hence the resulting graph must be undirected. Following standard procedures (cf. [76]), the structural connectivity matrix was symmetrized A ← A+A T 2 , as any deviations from symmetry are artifacts of the processing method. Since application of the PROBTRACKX algorithm typically results in a certain number of false positive fibers and thus assigns non-zero connection strengths to all pairs of brain regions, we enforced a sparsity of 20% [77] by discarding all connections with a relative connection strength smaller than 0.00071 (0.15% of strongest relative connection strength).

Supplemental Material
1 Whole-brain network -state space exploration 1.1 Limit cycle criterion Figure S1 shows traces of the activity variables x k1 and the respective total inputs for each node over time for different background inputs µ (see bifurcation diagram shown in Fig. 3b of the main paper). The total input to each node consists of the sum of the background input µ, the input σ N i A ki x i1 (t) from other nodes in the network, and the input η ξ k (t) from the noise process. We expect qualitatively different dynamics if the total input of a node lies inside vs. outside the interval given by Eq. (17) of the main paper (dashed blue lines in Fig. S1), because total inputs inside this interval correspond to an individual FHN oscillator being in its limit cycle.
For small µ, as shown in Fig. S1a, the total input is below the lower threshold, and the network is in a stable fixed point. If µ is larger with the total input still below but closer to the lower threshold, a low amplitude limit cycle (laLC) emerges (Fig. S1b). These oscillations cannot be explained by the dynamics of the individual FHN oscillator but emerge as a result of network effects. If the total inputs of some of the nodes enter the described interval, as it happens for some (Fig. S1c) or all (Fig. S1d) of the maxima of activity, the network dynamics exhibits high amplitude oscillations. The elevated amplitudes for nodes with total input within the interval increase the coupling inputs to other nodes in the network and this creates oscillations in nodes that are not in their individual LC. We call this regime the mixed limit cycle (mixLC) state. If µ is high enough so that all nodes satisfy the condition in Eq. (17), as shown in Fig. S1e, the network is in the high amplitude limit cycle (haLC) state. For higher µ, the total input of some nodes is exceeds the upper threshold and we observe mixLC states, as in Fig. S1f, that show similar aperiodic behavior as mixLC states at the low bifurcation (cf. Fig. S1c). We also observe a laLC at the high bifurcation (Fig. S1g) and a stable fixed point (Fig. S1h) for even higher values for the background input µ.
The interval criterion for the node input alone is not sufficient to explain the asynchronous or aperiodic network dynamics observed in Figs. S1c and f, but in these traces we see that the network dynamics change qualitatively, depending on whether some of the node inputs fall into the considered interval. Consequently, small changes in the input to individual nodes can result in large effects on the network dynamics. Figure S2a shows the activity variables x k1 for nodes k ∈ {21, 22, 23} plotted against each other for the asynchronous mixed state at the low bifurcation in Fig. S1c (equivalent to Fig. 3c(3) of the main paper). We observe open trajectories, which differ between different oscillation periods. This is a strong indication that the analyzed network state is not only asynchronous but also aperiodic. Figure S2b shows the equivalent phenomenon for the trajectories of the investigated asynchronous states at the low bifurcation S1 (cf. Fig. 3a of the main paper). Figure   point S1 (cf. Fig. 3a of the main paper). The large differences in ranges of the trajectories are caused by total inputs [Eq. (S1)] of different strength due to their different weighted in-degree. Even small changes in the input to a node can have a large effect on the shape of the limit cycle, especially close to the bifurcations (see also Fig. S3b). The starting and end points of these trajectories (all plotted for t ∈ [5000, 5064]) show that the oscillations of these nodes are not in phase and that the phase difference between them does not remain constant.

Influence of noise on the network dynamics
Figures S4a-h show the activity variables x k1 and the total inputs [Eq. (S1)] for each node for the same parametrizations of the network dynamics as Figs. S1a-h, with additional Gaussian white noise (η = 0.024). If the network without noise is in the stable down state close the the bifurcation, as in Fig. S1a, small fluctuations in the node input can lead to large changes in the activity variable. Such oscillations of individual nodes are then sustained by the coupling and drive the network oscillation. This effect blurs the distinction between FPs and the laLC for increasing noise levels. The oscillation amplitudes of nodes in the laLC fluctuate due to the noise, which makes it even harder to distinguish the states based on the time-series of activity variables. Interestingly, the mean network cross-correlation of the noise-free asynchronous state in Fig. S1c is increased from R = 0.47 (η = 0) to R = 0.58 (η = 0.024) in Fig. S4c. A potential explanation for the higher synchrony with additional Gaussian noise is the increased number of nodes that fulfill the criterion in Eq. (17) of the main paper (36 for η = 0.024, compared to 15 for η = 0). If a random fluctuation drives a node into the respective individual LC, this increases the coupling input to the other nodes in the network. With more nodes oscillating in the individual LC, the driving input to the remaining nodes is increased, which can then be sufficient to entrain a synchronous network oscillation. Close to the high bifurcation in Fig. S4f, a similar effect counteracts the aperiodic collapsing of the network oscillation (cf. Fig. S1f). The periodic and synchronous states in Figs. S4d, e are robust against the influence of noise, showing no qualitative change in their dynamics compared to the noise-free case in Figs. S1d, e. Figure S4h   at the low bifurcation, a self-sustained noise induced network oscillation at the high bifurcation. Figure S4i shows that if the network dynamics are in a fixed point far away from a bifurcation, we observe uncorrelated fluctuations in the activity variable and the combined input, instead of the emerging oscillations in Figs. S4a and h.
When evaluating the results of the optimal control method for the synchronization task in the case of finite noise, it is unreasonable to expect that the optimal control algorithm will achieve the target mean network cross-correlation of R T = 1. In the synchronization task in Section IV B of the main paper, the optimal control method increases the average network cross-correlation from R = 0.25 in the uncontrolled case to R = 0.83 (parameters µ = 0.7, σ = 0.025, η = 0.024, location S1, Fig. 10a of the main paper). We compare the deviation of this result from the target R T = 1 with the effect of noise on R of a previously synchronous network state. Figure S4j shows the effect of noise on an uncontrolled, synchronous state with low global coupling strength as S1 (σ = 0.025, used for synchronization task of the main paper). The mean network cross-correlation of R = 0.99 in the noise-free (η = 0) case is substantially reduced to R = 0.75 when noise is introduced (η = 0.024, same as for the synchronization task of the main paper).
2 Optimal control of the brain network dynamics

Recovery of laLC target state
The network dynamics as shown in Fig. 6b of the main paper does not achieve the target state (cf. initial state in Fig. 6a of the main paper) within the time interval in which the control is active. The biphasic pulse of the optimal control method occurs around t = 150, after which the network oscillation almost comes to rest. Since, however, no stable fixed point exists in the absence of control input, the oscillation amplitude grows until it reaches the laLC target state. Figure S5 shows that the laLC target state is stable and reached about 550 time units after the control pulse. Although the control target is eventually reached, this trajectory results in a nonzero precision cost, which is evaluated in the last τ = 25 time units of the control interval. Due to energy constraints, however, this is still the optimal control solution.

Computation of critical times t c
The critical time t c is defined as the first point in time when all nodes are "in synchrony". The cross-correlation measure for synchrony [Eq. (7) of the main paper] is defined for time intervals and does not allow to determine synchrony for a given point in time. We therefore compute the Kuramoto order parameter r, for N nodes of the network with their respective phases θ k . We calculate the phases for node k by assigning θ k (t) = 0 to all timest that correspond to its oscillation maxima. We then linearly interpolate the phases θ k ∈ [0, 2π) between subsequent maxima. The critical time t c is then defined as the first point in time in the control period, for which r(t) ≥ 0.999 holds. Figure S6 provides additional information to Fig. 12 of the main paper by showing the fraction of controlled nodes as a function of the sparsity parameter I s and the resulting average crosscorrelation R. As expected, the fraction of controlled nodes decreases with increasing penalty on sparsity in space. Plotting the average cross-correlations R against the fraction of controlled nodes intuitively shows how many nodes have to receive a finite control input in order to achieve a certain synchrony for different noise levels. It is not possible to achieve a significantly higher synchrony compared to the uncontrolled dynamics by controlling less than 15% of the network nodes.

Videos of the node dynamics in phase space with optimal control
Videos can be viewed at https://github.com/rederoth/nonlinearControlFHN-videos. The states (x k1 , x k2 ) for all network nodes k are plotted as dots in the (x k1 , x k2 ) phase plane. The color of each node corresponds to its weighted degree (green to blue indicates high to low values). The optimal control inputsū k are indicated by black arrows whose direction denotes the sign and whose length denote absolute strength (see top right corner for scale). The time t and the duration of the control interval is indicated in the bottom right corner. If applicable, the limit cycle of a single node (no coupling, σ = 0) is indicated in light gray.
• Video 1: Attractor switching at point A1 (low bifurcation) from down-state to oscillatory state Figure S6: Synchronization with sparse optimal control for different noise levels η, equivalent to Fig. 12 of the main paper. (a) Fraction of nodes in the network that receive finite control input as a function of the sparsity parameter I s at point S1 (low bifurcation). (b) Same for S2 (high bifurcation). (c) Achieved synchonization measured by the average cross-correlation R as a function of the fraction of controlled nodes at point S1. The mean and standard deviation are shown for 5 different initial conditions of the network dynamics, each with 20 independent realizations of the noise. Horizontal lines denote the average cross-correlations R for the uncontrolled case. (d) Same for S2. Parameters as in Fig. 12 of the main paper. • Video 2: Attractor switching at point A1 (low bifurcation) from oscillatory state to downstate • Video 3: Attractor switching at point B2 (high bifurcation) from low-amplitude oscillation to high-amplitude oscillation • Video 4: Attractor switching at point B2 (high bifurcation) from high-amplitude oscillation to low-amplitude oscillation • Video 5: Synchronization at point S1 (low bifurcation) with η = 0 (noise-free network) • Video 6: Synchronization at point S2 (high bifurcation) with η = 0 (noise-free network) • Video 7: Synchronization at point S1 (low bifurcation) with η = 0.024 (noisy network) • Video 8: Synchronization at point S2 (high bifurcation) with η = 0.024 (noisy network) 3 Linear network control theory 3.1 Correlation of average and modal controllability with weighted degree The considered linear controllability measures are known to be closely related to the weighted node degree d k [23,25]. Figure S7 shows the average and modal controllability measures as function of d k for each of the N = 94 network nodes, where the average controllability shows a positive, and the modal controllability a negative correlation with d k . Nodes with high average controllability tend to have a small modal controllability.
3.2 Modal controllability evaluated for the synchronisazion task (noisefree case) Figure S8 shows the correlation between the optimal node-wise control energies E k of the control inputs with the modal controllability for the synchronization task discussed in Fig. 14 of the main paper. Since nodes with high average controllability show a low modal controllability, the positive correlation at the low bifurcation for t < t c (Fig. 14a of the main paper) is reversed, while there is again no significant correlation at the high bifurcation (cf. Fig. S8b). In the second phase of the control, with t ≥ t c , we again observe nodes with an intermediate modal controllability receiving smaller control energies E k (analogous to Figs. 9c, d and 14c, d of the main paper). Figure S8: Mean of the optimal node-wise control energy E k plotted against the modal controllability for the synchronization task (equivalent to Fig. 14   . All parameters as in Fig. 11 of the main paper.

Controllability measures for the synchronization task (noisy network dynamics)
For the synchronization task under finite noise, we observe a change in the control strategy with increasing noise levels (cf. Section II B of the main paper). Figure S9 shows that the influence of noise also affects the correlation of the node-wise control energy E k with the introduced measures from linear control theory. At the low bifurcation (Fig. S9a), E k is not correlated with the average controllability over the whole time interval of the control for η = 0 (combination of Figs. 14a and c of the main paper). With increasing noise levels, however, we also observe an increased negative correlation. At the high bifurcation (Fig. S9b) the correlation of E k with the average controllability is even inverted from positive for η = 0 to increasingly negative for η = 0.012 and η = 0.024. Similar effects are observed for the changes in the changes in the correlation between the node-wise control energy E k and the modal controllability (Fig. S9c, d) when noise is increased.