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\times N$ dimensional adjacency matrix $\bm{A}$ which is scaled by a global coupling strength $\sigma$. We allow for an instantaneous and additive control input to the network, which is described by $N\times d$ independent control variables $\bm{u}=(\bm{u_{1}},\ldots,\bm{u_{N}})$ with $\bm{u_{i}}=(u_{i1},\dots,u_{id})$. The equations describing the controlled network dynamics thus read:

where $\otimes$ denotes the Kronecker product. The state of the network is described by a vector $\bm{x}=(\bm{x_{1}},\ldots,\bm{x_{N}})$, where $\bm{x_{i}}=(x_{i1},\dots,x_{id})$ characterize the individual states of the nodes. The local node dynamics is given by $\bm{h}(\bm{x})=(\bm{h}(\bm{x_{1}}),\ldots,\bm{h}(\bm{x_{N}}))$ with $\bm{h}(\bm{x_{i}})=(h_{1}(\bm{x_{i}}),\dots,h_{d}(\bm{x_{i}}))$. All nodes additionally receive Gaussian white noise of similar intensity $\eta$, which may be correlated within but is uncorrelated across the nodes of the network. The stochastic variables $\bm{\xi}=(\bm{\xi_{1}},\dots,\bm{\xi_{N}})$ with $\bm{\xi_{i}}=(\xi_{i1},\dots,\xi_{id})$ are independently drawn from a Gaussian distribution, $\xi_{ij}\in\mathcal{N}(0,1)$, with zero mean and unit variance. Within node correlations are quantified by the $d\times d$ dimensional local noise scheme $\bm{D}$, while the across-node statistical independence is assured by a N-dimensional identity matrix $\bm{I_{N}}$. The internode coupling term consists of the Kronecker product of the adjacency matrix $\bm{A}$ and the $d\times d$ dimensional local coupling scheme $\bm{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\times N$ matrix $\bm{B}$ in the control term allows for the control of multiple nodes with different strength. The $d\times d$ dimensional matrix $\bm{K}$ implements the local control scheme, which is similar for every node in the network. Initial conditions are denoted by $\bm{x}(t=0)=\bm{x_{0}}$.

The control $\bm{u}$ is considered to be optimal ($\bm{u}=\overline{\bm{u}}$), if it minimizes an appropriate cost functional. To this end, we construct a cost functional $F(\bm{x}(\bm{u}),\bm{u})$ for a state switching task, where the control input drives the network model from one stable state to another (see Section IV.1), and for a node synchronization task, where the control input increases the degree of synchronization among its nodes (see Section IV.2). For finite noise, i.e. for finite values of $\eta$ in Eq. (1), the overall cost functional $F$ is defined as a mean over noise realizations [30],

where $F_{n}$ denotes the cost functional for one noise realization $n$ and the angle brackets $\langle\cdot\rangle$ denote the mean. In the noise-free case no averaging is performed. The state dependent term $F_{n}^{x}(\bm{x})$ only implicitly depends on the control and penalizes the deviation from the desired output. It will be different for the switching, $F_{n,\textrm{{sw}}}^{x}$, and synchronization, $F_{n,\textrm{{syn}}}^{x}$, tasks (see below). The control dependent term $F^{u}(\bm{u})$ accounts for the cost of the control itself.

For the state switching task, we consider a noise-free system ($\eta=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 $\bm{x}_{T}(t)$:

where we consider the control being active within the time period $0\leq t\leq 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.1).

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\in[0,\dots,N]$ and where $\langle.\rangle_{t}$ and $\sigma_{\langle.\rangle_{t}}$ denotes the temporal mean and the standard deviation. The mean $R^{m}$ over all values $R_{ij}^{m}$

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}(\bm{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 L2-regulation 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\leq t\leq 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 $\overline{\bm{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. [31, 29, 50], we analyze the gradient $g=\bm{\nabla_{u}}F$ of the cost functional, which has to vanish when evaluated at the optimal control $\overline{\bm{u}}$. By applying the method of Lagrange multipliers, we obtain an expression for the optimality condition that depends on the adjoint states $\bm{\phi}(\bm{x},{\bm{u}},t)$ corresponding to the Lagrange multipliers. The control $\overline{\bm{u}}$ is optimal, if

for $0\leq t\leq T$, where $\overline{\lambda}(t)$ is the derivative of the sparsity term, Eq. (7), with respect to the control inputs $\bm{u}$. The adjoint states are governed by a set of linear differential equations:

where $D_{\bm{x}}(\bm{h})$ is the Jacobian matrix of the state equations of the dynamical system, and $f^{x}_{n}(\bm{x})$ is the integrand of the cost functional, $F^{x}_{n}(\bm{x})=\int_{0}^{T}f^{x}_{n}(\bm{x})dt$. The adjoint state satisfies the boundary condition $\bm{\phi}(T)=\bm{0}$, and the differential equation is solved backwards in time.

Following Ref. [28], $\overline{\lambda}(t)$ is given by

Here, $k\in[1,N]$, and $E_{k}$ is the nodewise control energy of the resulting optimal control $\overline{\bm{u}}$, which is defined as

resulting in the total control energy $E=\sum_{k}E_{k}$. A detailed derivation of Eqs. (9) and (10) is provided in Appendix A.

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 $\bm{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 $\bm{u}_{0}$ and drawing $20$ noise realizations $\bm{\xi}_{n}(t)$. The corresponding states $\bm{x}_{0}(t)$, Eq. (1), and adjoint states $\bm{\phi}_{0}(t)$, Eq. (10), are calculated for every noise realization. Then $\bm{g}(t)$ as given in Eq. (9) is evaluated. The descent direction is initialized with $\bm{d}_{0}(t)=-\bm{g}(t)$. We then iterate until convergence:

• while $\|\bm{g}(t)\|_{\infty}>\epsilon$:

  1. 1.

     compute the step size $s_{l}$ using bisection

  2. 2.

     set $\bm{u}_{l+1}(t)=\bm{u}_{l}(t)+s_{l}\bm{d}_{l}(t)$

  3. 3.

     calculate the states $\bm{x}_{l+1}(t)$ and the adjoint states $\bm{\phi}_{l+1}(t)$ for every noise realization and compute the mean $\langle\bm{\phi}_{l+1}(t)\rangle$

  4. 4.

     evaluate the gradient $\bm{g}_{l+1}(t)$

  5. 5.

     compute $\beta_{l}$ using the Polak-Ribiere method: $\beta_{l}=\frac{\bm{g}_{l+1}(\bm{g}_{l+1}-\bm{g}_{l})}{\|\bm{g}_{l+1}\|^{2}}$

  6. 6.

     set the direction $\bm{d}_{l+1}(t)=-\bm{g}_{l+1}(t)+\beta_{l}\bm{d}_{l}(t)$

  7. 7.

     if $\bm{d}_{l+1}(t)$ is not a descent direction, set $\bm{d}_{l+1}(t)=-\bm{g}_{l+1}(t)$

  8. 8.

     set $l=l+1$
     

Gradient-based optimization is only guaranteed to converge to local optima of the cost function. In general, however, multiple initial conditions $\bm{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 $\bm{u}$ coexists with the solution $\bm{u}(t)=0$ (see Appendix B for details).

We consider a single FitzHugh-Nagumo (FHN) oscillator, with the activity variable $x_{1}$ and the linear recovery variable $x_{2}$:

where $\mu$ is a node-independent, constant background input and $R(x)=-\alpha x^{3}+\beta x^{2}-\gamma x$. The parameters in $R$ are chosen to obtain the bifurcation diagram shown in Fig. 1. This bifurcation diagram, depending on the background input $\mu$, shows three distinct states. In the down 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.

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 $\bm{A}$ shown in Fig. 2. Since DTI data carries no directional information, $\bm{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}=\sum_{i=1}^{N}\bm{A}_{ik}$.

As motivated in Section III.1, 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\in[1,N]$. This implies setting the control matrix $\bm{B}=\bm{I_{N}}$, which means that all individual nodes can potentially receive independent control inputs. The local coupling, control, and noise schemes are set to $\bm{G}=\bm{K}=\bm{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 V 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.

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\,\,\forall\,k,t$].

Initially we consider the noise free case ($\eta=0$). Here, only two free model parameters remain: the global coupling strength $\sigma$ and the time independent background input $\mu$ 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 $\bm{x_{0}}$ are drawn randomly from a uniform distribution in the interval $\mathopen{[}0,1\mathclose{)}$ or taken as the state vector $\bm{x^{\prime}}(t_{end})$ of the last time step of the previous simulation with same $\sigma$ and slightly smaller $\mu$ (continuation). In order to avoid analyzing transient effects, we show and evaluate the network states only after a sufficiently long transient time $\bar{t}$.

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}=\left|(\mathcal{F}x_{i1})(f)\right|^{2},$ with Fourier transform $\mathcal{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 $\mu$ and $\sigma$ the activity variable $x_{k1}$ has a low value for all nodes [down state, cf. Fig. 3c(1)]; for large values of $\mu$ and $\sigma$ the activity value is high [up state, cf. Fig. 3c(8)]. 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 $\mu$ or $\sigma$. 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 $\mu\in[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 $\sigma$ 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 effects 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 lead – even 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 $\bm{x_{0}}$ and comparing the resulting time series by calculating their node-wise correlation in the activity variable. We ignore the first $\bar{t}=5000$ arb. units of each time series to avoid analyzing transient effects and then compare the interval $t\in[5000,6000]$ (equiv. to 10–35 periods) of one initial condition with $t\in[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 noise-free 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].

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-\tau)\leq t\leq T$ and $I_{p}=0$ else. The convergence criterion of the numerical solution of the minimization problem, as described in Section II.3, is set to $\epsilon=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 $\epsilon$. 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}(\overline{u}_{1})=35.0$] corresponds to the frequency of the fixed point’s focus [$f_{dom}(x^{FP}_{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). When switching from the laLC to the mixLC state (Figs. 6a, c), the frequency of the control input [$f_{dom}(\overline{u}_{1})=30.0$] is again adapted to the frequency of the initial state [$f_{dom}(x^{laLC}_{1})=30.0$], utilizing resonance effects. While the nodes in the laLC state all oscillate with the same frequency ($f_{k}=f_{dom}=30.0\ \forall\ k$), their phase speed is not the same at all times [reflected by lower synchrony $R(x^{laLC}_{1})=0.861$, best observed in Video 3 in [78]]. This results in a less synchronous oscillation of the control signals in Fig. 6c ($R(\overline{u}_{1})=0.714$) compared to the control inputs at the low bifurcation [Fig. 5c, $R(\overline{u}_{1})=0.999$]. In the other direction (Figs. 6b, d, Video 4 in [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\in[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 $\overline{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 driver nodes for attractor switching. At the high bifurcation (Figs. 7b, d), on the other hand, nodes with low weighted degree receive the strongest control inputs.

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.

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: $\langle f\rangle_{N}=32.2,\,\sigma_{\langle f\rangle_{N}}=1.1,$ and at S2: $\langle f\rangle_{N}=27.0,\,\sigma_{\langle f\rangle_{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 $\mu$ 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\geq 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\in[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 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 $\eta$, we observe a gradual transition of the optimal control strategy. Instead of balancing out 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 $\bm{x_{0}}$ (cf. Fig. S1 in [78]), as it depends on details of the exact network dynamics.