Front. Comput. Neurosci.Frontiers in Computational NeuroscienceFront. Comput. Neurosci.1662-5188Frontiers Media S.A.10.3389/fncom.2024.1483100NeuroscienceOriginal ResearchA framework for optimal control of oscillations and synchrony applied to non-linear models of neural population dynamicsSalfenmoserLena1*ObermayerKlaus121Institute of Software Engineering and Theoretical Computer Science, Technische Universitaet Berlin, Berlin, Germany2Bernstein Center for Computational Neuroscience Berlin, Berlin, Germany
Edited by: Qingdu Li, University of Shanghai for Science and Technology, China
Reviewed by: Qingyun Wang, Beihang University, China
Jr-Shin Li, Washington University in St. Louis, United States
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
We adapt non-linear optimal control theory (OCT) to control oscillations and network synchrony and apply it to models of neural population dynamics. OCT is a mathematical framework to compute an efficient stimulation for dynamical systems. In its standard formulation, it requires a well-defined reference trajectory as target state. This requirement, however, may be overly restrictive for oscillatory targets, where the exact trajectory shape might not be relevant. To overcome this limitation, we introduce three alternative cost functionals to target oscillations and synchrony without specification of a reference trajectory. We successfully apply these cost functionals to single-node and network models of neural populations, in which each node is described by either the Wilson-Cowan model or a biophysically realistic high-dimensional mean-field model of exponential integrate-and-fire neurons. We compute efficient control strategies for four different control tasks. First, we drive oscillations from a stable stationary state at a particular frequency. Second, we switch between stationary and oscillatory stable states and find a translational invariance of the state-switching control signals. Third, we switch between in-phase and out-of-phase oscillations in a two-node network, where all cost functionals lead to identical OC signals in the minimum-energy limit. Finally, we (de-) synchronize an (a-) synchronously oscillating six-node network. In this setup, for the desynchronization task, we find very different control strategies for the three cost functionals. The suggested methods represent a toolbox that enables to include oscillatory phenomena into the framework of non-linear OCT without specification of an exact reference trajectory. However, task-specific adjustments of the optimization parameters have to be performed to obtain informative results.
nonlinear optimal controlcontrol of oscillationscontrol of synchronycontrol of neural dynamicsneural population modelsOB 102/31Deutsche Forschungsgemeinschaft10.13039/5011000016591 Introduction
OCT is a mathematical framework that offers methods to compute efficient stimulation for linear or non-linear dynamical systems (Kirk, 2004) with numerous applications in science, engineering, and operations research. With OCT, a stimulating signal can be designed such that the dynamical system under consideration behaves in a desired manner. Conventionally, a reference trajectory is defined that the state variable of the dynamical system should mimic. A cost functional trades the closeness between the actual and a reference state (typically, the squared difference between the reference and target states integrated over time) against the input strength of the control (typically, the squared control input integrated over time). The optimal control (OC) minimizes the cost functional.
Control of neural system has been studied intensively both theoretically and experimentally (Kao and Hennequin, 2019; Suppa et al., 2016; Liu et al., 2018). Closed-loop control and machine learning approaches have been applied, for example, to modulate brain activity (Grosenick et al., 2015; Tafazoli et al., 2020; Park et al., 2019) and optimize stimulation protocols for the treatment of neurological disorders (Yu et al., 2020; Chandrabhatla et al., 2023). Open-loop optimal control has been computed for non-linear models of neural population dynamics (Salfenmoser and Obermayer, 2022, 2023) and networks thereof Chouzouris et al. (2021). Beyond that, network control theory has widely been applied to linear models of brain networks to understand how network structure determines the role of brain circuits to control network function (McGowan et al., 2022), or how altered structural connectivity is linked to mental disorders (Zöller et al., 2021).
In neural systems, oscillatory phenomena commonly emerge and spread over various spatial and temporal scales. To control neural oscillations, OCT has mostly been applied to phase-reduced models of oscillatory systems (e.g., Dasanayake and Li, 2011; Zlotnik and Li, 2012; Pyragas et al., 2020; Wilson and Moehlis, 2014). Phase reduction techniques enable to drastically simplify complex dynamical system evolving along a stable limit cycle by parametrizing their dynamics via their phase only (Pietras and Daffertshofer, 2019). As the oscillation phase becomes the (only) state variable, synchronization directly translates to phase alignment, and optimal synchronization and desynchronization can be studied straightforwardly by conventional OCT methods. The low computational complexity of phase-reduced models enables OC studies of large-scale networks (e.g., Bomela et al., 2023). However, phase-reduced models can by construction only describe trajectories along the stable limit cycles of coupled oscillators. Control strategies are limited to phase shifts along these limit cycles and cannot capture larger deviations from it. The phase-reduced approximation breaks down at a bifurcation, and the model dynamics may not be described accurately close to a bifurcation. Beyond phase model studies, delayed feedback control (Popovych et al., 2006; Hövel et al., 2010) and model-free machine learning methods (Krylov et al., 2020; Vu et al., 2024) have been applied to discuss optimal synchronization or desynchronization in multi-dimensional models of neural activity.
There are different classes of methods to compute OC for non-linear systems (Rao, 2010). Model-based approaches use the mathematical description of a dynamical system: While indirect methods reformulate the optimization problem such that optimality conditions are obtained, from which the optimal control can be computed, direct methods discretize the continuous control problem such that it becomes an optimization problem in finite space. Model-free approaches learn control strategies from a sufficient amount of data of how the system behaves under stimulation, without using any information on the analytical model description (Koryakovskiy et al., 2017).
In this study, we compute OC using the indirect adjoint method, which enables to formulate an expression for the gradient of the cost functional with respect to the control. It can be applied to any linear or non-linear dynamical system and is, in principle, analytic. However, the arising equations cannot be solved analytically in most non-linear systems and require numerical solutions. Up to the uncertainties of the numerical computation, the method remains exact. Control shapes are entirely determined by the system dynamics at any point in the state space and for any dynamical regime. The method is applicable to arbitrary networks; however, in practice, computational complexity will limit network sizes.
Standard mathematical methods borrowed from OCT are limited in their applicability to oscillatory phenomena in (unreduced) complex dynamical systems: The specification of a reference trajectory is overly restrictive since the exact shape of the oscillatory trajectory, including properties such as amplitude and phase, may often be irrelevant. We therefore adapt the OCT framework and suggest to replace the cost functional that measures closeness between actual and reference state by one of three alternative cost functionals:
The Fourier cost evaluates the Fourier component corresponding to a target oscillation frequency.
The cross-correlation cost evaluates the pairwise cross-correlation of the dynamical state of all nodes in a network (Chouzouris et al., 2021).
The variance cost evaluates the variance of the activity throughout a network.
The Fourier cost enables new applications of OCT to induce oscillations when the system is in a fixed point without specifying a reference trajectory. All three cost functionals are evaluated and compared for switching between stationary states and for the control of network (de-) synchronization. The variance cost is inspired by studies on the control of frequency synchronization in phase oscillator models (e.g., Taher et al. (2019)), and a similar cost functional has recently been studied using a data-driven approach (see Vu et al., 2024). The cross-correlation cost has been studied before in the context of network synchronization (cf. Chouzouris et al., 2021) and is included for comparison.
We demonstrate how to include these cost functionals into the adjoint method, such that all benefits from this method are inherited. This enables to study efficient control strategies in complex systems for a repertoire of control tasks that cannot be captured by phase-reduced oscillator models. Potential applications go beyond (de-) synchronization of oscillating networks and may combine oscillatory and non-oscillatory states, e.g., in multistable systems, or across bifurcation boundaries.
In this study, we use two non-linear models of neural populations dynamics as example systems. We study the Wilson-Cowan (WC) model (Wilson and Cowan, 1972, 1973), a simple, phenomenological neural mass model, and the mean-field EIF model, a biophysically motivated mean-field model of excitatory and inhibitory exponential integrate-and-fire (EIF) neurons (Augustin et al., 2017; Cakan and Obermayer, 2020). The mean-field EIF model exhibits more dynamical variety and hence enables to study control tasks that are not accessible with the WC model only. In these models, the OC cannot be computed analytically. Instead, the cost gradient is computed numerically, and the OC is approached by gradient descent.
This article is structured as follows: Section 2.1 provides an introduction to OCT, and Section 2.2 presents our adjustments to target the problem of controlling oscillations and synchrony. In Section 2.3, we introduce the two models of neural population dynamics that serve as example systems in this study. Section 2.4 presents the implementation of our OC methods, and we provide details on the numerical simulations and on the convergence properties of the gradient descent procedure. In Section 3, we present our results. We first study the efficacy of the Fourier cost to enforce oscillations from stationary states in a single-node system (see Sections 3.1 and 3.2.1). Then, we apply the Fourier cost, the cross-correlation cost, and the variance cost to switch between oscillatory states (see Section 3.2.2). Finally, we synchronize or desynchronize oscillating networks with the three proposed cost functionals (see Section 3.3). In Section 4, we summarize our findings, discuss the limitations of our method, and give an outlook on potential future areas of application.
2 Methods2.1 Optimal control of non-linear dynamical systems
We study controlled non-linear dynamical systems that may be defined by ordinary differential equations (ODEs) or differential-algebraic equations (DAEs) and may include delays. We denote such systems by
with the state vector x→, its time derivative x→˙, and Nd potentially different delays di. In this study, we consider network systems with N network nodes and a finite simulation time. The state vector x→ contains the stacked state vectors x→n,n∈{1,2,...N}, of all nodes. Similarly, h→ contains the stacked system dynamics of all network nodes. We denote the dimensionality of the state vector of a single node as Nx, such that {x→,h→}∈ℝN·Nx. u→(t) is the control input, which may effect one or more components of x→n for each node in a network. The OC formalism provides methods to compute a control u→(t) that is most efficient for a particular purpose. A cost functional F measures the efficiency of a control. It is conventionally defined as
The first term FP measures the closeness of the trajectory of the state vector x→(t) to a given target trajectory of state vectors x~(t)∈ℝN·Nx and is referred to as the precision cost. The second term FE measures the strength of the control signal and is referred to as the energy cost. T denotes the considered time interval. Inaccuracy is penalized for t0 ≤ t ≤ T, while control strength is penalized for the complete time interval. When controlling a switch between two stable states, we set t0>0 to account for a finite transition time. The OC u→* minimizes the cost,
u→*=argminu→F.
The weight wP can be tuned to trade precision against the required control strength.
In the following Section 2.2, we will consider new cost functionals FX(x→) which go beyond the precision cost FP(x→) defined in Equation 2, i.e., FP(x→)→FX(x→). To be able to apply the adjoint method, we require that the derivative of the cost functional with respect to the control can be written in the form
dFX(x→)du→=∫0Tg→(x→)dx→du→dt,
where g→(x→) is some function of the state variable x→ (see Supplementary material).
The adjoint method (Kirk, 2004) enables to compute the derivative ddu→F of the cost functional with respect to the control. The adjoint state λ→∈ℝN·Nx is defined by the differential equation
with the final condition λ→(T)=0. χ[ta, tb] denotes the indicator function of the time interval [ta, tb]. Above expression contains vectors {λ→,λ→˙,g→(x→)}∈ℝN·Nx, and matrices {∂h→∂x→˙,∂h→∂x→}∈ℝ(N·Nx)×(N·Nx), with (∂h→∂x→˙)ij=∂hi∂ẋj (indices i and j denote the ith and jth vector components). Note that for the dynamical systems studied in this work, ∂h→∂x→˙ is a constant such that the last term in Equation 5 vanishes. The adjoint state λ→ is computed by backward integration. It enables to compute the cost derivative
ddu→F=∫0T∂f∂u→+λ→T∂h→∂u→dt.
The integrand ∂f∂u→+λ→T∂h→∂u→ is the cost gradient as a function of time, which is required for the gradient descent. A derivation of the adjoint method is provided in the Supplementary Section 1.
When studying oscillatory phenomena and the control thereof, a well-defined target trajectory x~n(t) for each network node implies that the exact shape of the oscillation, including amplitude and phase, is set. This might be overly restrictive in certain scenarios. In the next section, we will, therefore, introduce alternative cost functionals to replace the precision cost FP in Equation 2.
2.2 Optimal control of oscillations and synchrony
For the cost functionals that measure oscillation and synchronization, we consider only one relevant state variable for each network node. The formalism we present is, however, not generally limited to this simplification and can be extended straightforwardly. We drop the vector arrow and denote the observable component of the state variable by xn(t) for network node n.
To enforce phase and frequency synchronization in a network at frequency f~, we replace FP in Equation 2 by the synchronization Fourier cost,
FFsync=−1N2(T−t0)2|∫t0TΣx(t)·e−iωtdt|2
with Σx(t)=∑n=1Nxn(t) and ω=2πf~. FFsync is the squared Fourier component corresponding to the frequency f~ of the sum of observables over all N network nodes. The computation of the adjoint state using Equation 5 requires
(gFsync(x→))k=-2N2(T-t0)2∫t0TΣx(τ)·cos(ω(t-τ))dτ.
The derivation of Equation 8 is provided in the Supplementary Section 2.1.
To enforce frequency synchronization in a network at frequency f~ irrespective of the relative phases of its nodes, we replace FP in Equation 2 by the oscillation Fourier cost
FFosc=−1N(T−t0)2∑n=1N|∫t0Txn(t)·e−iωtdt|2.
FFosc is the sum over all N network nodes of the node-wise squared Fourier component corresponding to the frequency f~. The computation of the adjoint state using Equation 5 requires
(gFosc(x→))k=-2N(T-t0)2∫t0Txk(τ)·cos(ω(t-τ))dτ.
The derivation of Equation 10 is provided in the Supplementary Section 2.2. For N = 1, synchronization and oscillation Fourier costs are identical and enforce oscillations at frequency f~. We will denote the cost by FF in this case. The oscillation Fourier cost FFosc improves, if the power corresponding to the target frequency f~ increases in any network node, irrespective of the relative phases. The synchronization Fourier cost FFsync improves, either if oscillatory network activity at frequency f~ synchronizes, or if the power in the network nodes corresponding to f~ increases1.
To enforce phase synchronization in a network irrespective of the oscillation frequency, we replace FP in Equation 2 by either the cross-correlation cost or the variance cost. The cross-correlation cost (Chouzouris et al., 2021) is defined as
where x¯=1T-t0∫t0Tx(t)dt is the temporal mean of x(t), and σ2(x)=1T-t0∫t0T(x(t)-x¯)2dt is its variance. The computation of the adjoint state using Equation 5 requires
with Ikl=∫t0T(xk(t)-x¯k)(xl(t)-x¯l)dt. The derivation of Equation 12 is provided in the Supplementary Section 2.3.
The variance cost is defined as
Fvar=1N(T-t0)∫t0T∑n=1N(xn(t)-x¯(t))2dt,
where x¯(t)=1N∑n=1Nxn(t) is the network mean of x(t). The computation of the adjoint state using Equation 5 requires
(gvar(x→))k=2(xk(t)-x¯(t))N(T-t0).
The derivation of Equation 14 is provided in the Supplementary Section 2.4.
Both Fcc and Fvar are defined for networks with N≥2 and cannot be applied to a single-node system.
2.2.1 Control weights
When replacing FP in Equation 2 by one of the abovementioned cost functionals, we also introduce the respective weights, wFsync, wFosc, wcc, or wvar replacing wP in Equation 2. By using a negative weight, we can suppress a specific oscillation mode (oscillation Fourier cost) and enforce desynchronization instead of synchronization (synchronization Fourier cost), enforce a small cross-correlation, or enforce a large variance (e.g., in the case of a desynchronization task).
2.3 Models of neural population dynamics
We apply our methods to two models of neural population dynamics, the WC model (see Section 2.3.1) and the mean-field EIF model (see Section 2.3.2). The former is a simple, two-dimensional system of ordinary differential equations (ODEs), while the latter is a high-dimensional system of delay differential-algebraic equations (DDAEs).
2.3.1 The Wilson-Cowan model
The WC model describes the activity of coupled excitatory and inhibitory neural populations (Wilson and Cowan, 1972, 1973). The dynamics are defined by
The activity variables En(t) and In(t) denote the fraction of excitatory and inhibitory neurons that are active at time t in node n in a network of N nodes. τE and τI are the decay time constants of the excitatory and inhibitory activity, respectively. cβα denotes the coupling from population α to population β with α, β∈{E, I}. Eext and Iext are static external inputs that serve as order parameters for the network. un(t) is the time-dependent control, cgl denotes the global coupling strength, and C and D are the N×N coupling and delay matrices. The transfer function S(x) = (1+exp−γ(x−μ))−1 computes a synaptic input to En or In from the sum of all inputs. Numerical values for the parameters are given in Table 1. Figure 1 sketches the dynamical interactions between the excitatory and inhibitory populations of a single WC node in a network. Excitatory and inhibitory population are recurrently coupled. Only the excitatory population receives network and control inputs.
Parameters of the WC model.
Parameter
Explanation
Value
τE
decay constant of E
2.5
τI
decay constant of I
3.75
γ
gain
1.5
μ
firing threshold
3.0
cEE
E to E coupling
16
cEI
I to E coupling
12
cIE
E to I coupling
15
cII
I to I coupling
3
Cartoon of a single node of the WC network model. Static external inputs are indicated by Eext and Iext, network inputs by Enw, and control inputs by u(t).
Figure 2A shows a slice of the state space of the one-node WC model with parameters as given in Table 1. Depending on the external inputs Eext and Iext, the system exhibits a state of constant low activity (“down state”), a state of constant high activity (“up state”), an oscillatory state, or a bistable state, where stable states of constant low and high activity coexist. Points (A)–(E) mark locations in state space, for which control tasks are studied (see Section 3).
(A) Slice of state space of one node of the WC model. The horizontal (vertical) axis corresponds to the external input Eext (Iext). The four dynamical regimes are separated by black lines. Red markers indicate the points in state space at which control tasks are studied. Their coordinates are (1, 1) for point (A); (3, 1) for point (B); (1.8, 0.8) for point (C); (1.6, 0.4) for point (D); (1, 0.4) for point (E). (B) Slice of state space of one node of the mean-field EIF model. The horizontal (vertical) axis corresponds to the external input μEext (μIext). The five dynamical regimes are separated by black and blue lines. In the inset, we zoom into the regime bounded by the blue line, where stable oscillations coexist with a stable up state. Red markers indicate the points in state space at which control tasks are studied. Their coordinates are (0.2 nA, 0.2 nA) for point (A); (0.4 nA, 0.04 nA) for point (B); (0.38 nA, 0.3 nA) for point (C).
2.3.2 The mean-field EIF model
The mean-field EIF model (Augustin et al., 2017; Cakan and Obermayer, 2020) captures the dynamics of a network of randomly and sparsely connected exponential excitatory and inhibitory integrate-and-fire neurons (Brette and Gerstner, 2005) in the limit of infinitely many neurons. In this limit, a dimensionality reduction can be performed, and the mean firing rates rE and rI of the excitatory and inhibitory populations can be obtained as the activity state variables. rE and rI are functions of the mean membrane currents and their variances. The model dynamics are described by a 16-dimensional system of delay differential-algebraic equations (DDAEs), which is provided in the Supplementary Section 3. Schematically, the dynamical interactions of the mean-field EIF model can be represented similarly as shown in Figure 1. Note that external inputs are denoted differently as μEext and μIext.
Figure 2B shows a slice of state space of the one-node mean-field EIF model. Depending on the external inputs μEext and μIext, the system exhibits a state of constant low activity (“down state”), a state of constant high activity (“up state”), an oscillatory state, or a bistable state, where stable states of constant low and high activity coexist. Beyond that, there is a small region between the oscillatory and bistable regimes where we find bistability between an oscillatory and a down state (marked in blue). Points (A)–(C) mark locations in state space, for which control tasks are studied (see Section 3).
2.3.3 Parameters for network control
We study two control tasks for a WC network with nodes coupled via the excitatory population. In a preliminary state-space exploration for various combinations of network and coupling parameters, we hand-picked adequate parameters such that the desired control tasks can be studied.
For the two-node network studied in Section 3.2.2, coupling and delay matrices are given by
C=(0110),D=(09.59.50).
The coupling scheme is sketched in Figure 3A. We set cgl = 1.8 (cf. Equation 15) and find a bistable state at point (C) (see Figure 2A), where a stable in-phase oscillation (IP) with period fIP-1=13.89 coexists with a stable out-of-phase oscillation (OOP) with period fOOP-1=22.72.
(A) Cartoon of the neural mass model with two WC nodes studied in Section 3.2.2. Static external inputs are indicated by Eext and Iext and control inputs by u0,1 (t). The nodes are coupled via the excitatory population with a global coupling strength cgl. (B) Six-node system studied in Section 3.3.
For the six-node network studied in Section 3.3, coupling and delay matrices are given by
The coupling scheme is sketched in Figure 3B. We set cgl = 0.8 (cf. Equation 15) and find a state of asynchronous oscillation at point (D) and a state of synchronous oscillation at point (E) (see Figure 2A).
Note that the state-space diagram in Figure 2A only shows the bifurcations boundaries of the one-node system. Bifurcation boundaries change depending on the choice of network and coupling parameters.
Parameters were chosen such that relevant control problems could be defined for evaluating the proposed cost functionals. Note, however, that the presented method can be applied to any choice of coupling and delay matrices.
2.4 Implementation and numerical simulations2.4.1 Open-source implementation
Numerical computations are based on the neurolib simulation framework (Cakan et al., 2023). neurolib is an open-source library that contains several models of neural population dynamics and enables to combine them to network models of arbitrary size and structure. Beyond that, methods to compute OC in its standard formulation (see Section 2.1) are included. In this study, we expand neurolib by an extension that enables to compute OC signals for oscillation and synchronization tasks. This extension is openly available on GitHub at https://github.com/lenasal/neurolib/tree/OC_osc_sync.
2.4.2 Simulation accuracy
All simulations in this study have a duration of several hundred time units (dimensionless in the WC model, ms in the mean-field EIF model). We chose an integration step size dt = 0.1 (WC) or dt = 0.1ms (mean-field EIF) and validated that there are no qualitative differences if a smaller integration time step is used.
2.4.3 Fourier cost in limited-time simulations
When evaluating the Fourier spectrum of a trajectory obtained in a limited-time simulation, the spectral peaks are the sharper, the longer the simulation duration. If, however, the simulation time is relatively short (as in the control tasks considered in this study), the Fourier spectrum of any oscillation shows broad peaks (see Supplementary Figure S1). For tasks, in which we want to switch to a specific oscillatory state (Section 3.2), this leads to a tolerance against variations in the frequency f~ compared to the natural frequency of the target state.
2.4.4 Initialization
As an initial control guess, we use u0 = 0 in most tasks investigated in Section 3, with two exceptions.
First, in Section 3.1, we control the mean-field EIF model in its stationary down state. In this state, we encounter numerical problems when the gradient of the cost functional ∂f∂u→ is computed. This is due to the fact that the activity, the excitatory firing rate rE(t), is not given in terms of an analytical function, but in terms of a transfer function, for which no analytical closed-form expression exists (for details on the mathematical model, see also Supplementary Section 3). In simulations, rE(t) is, therefore, interpolated using a pre-computed table. In the down state, gradients almost vanish, preventing an effective OC computation if the control is not initialized reasonably. Hence, we initialize the OC algorithm with a sinusoidal initial control u0 which oscillates at the target frequency f~.
Second, in Section 3.2.2, we study a control task where all network nodes exhibit the exact same time evolution in the uncontrolled, initial state. When applying the variance cost, both cost and cost gradient vanish (see Equations 13 and 14), and the gradient descent cannot be performed. To circumvent this problem, we initialize the OC algorithm with an initial control that randomly fluctuates around zero with a small amplitude.
2.4.5 Choice of weights
The numerical values of the weights wF, wcc, or wvar determine how accuracy and control strength are traded against each other. On the one hand, when choosing the weight below a certain threshold value, which is individual for each cost functional and task, the OC equals zero, u→*(t)=0. Any finite input would increase the total cost via FE more than it might improve FF, Fcc, or Fvar. This threshold value is a lower limit to the choice of weights. On the other hand, when choosing a very large weight, the algorithm will return very strong control signals, and the controlled activity will reach its upper limits determined by the system dynamics. Physically, such results might be difficult to interpret. Computationally, the numerical integration might fail when the activity variable is continuously pushed against its upper bound. Hence, too large weight values must be avoided.
We study two different types of control tasks: First, there are control tasks that force a system to behave in a way that cannot be maintained naturally (driving oscillations from a stationary state, Sections 3.1, and synchronization (desynchronization) of desynchronized (synchronized) oscillations, Section 3.3). For these tasks, the weights are chosen after preliminary investigations such that control signals are reasonably strong. Second, there are control tasks in which we initiate a switch between two stable states (Section 3.2). For these tasks, we want to find the minimum-energy transition. Therefore, we dynamically adjust the weights, starting with a relatively large value of wF, wcc, or wvar, and decrease the respective weight every few (hundred) iterations until convergence, such that low-energy transitions are enforced. For these cases, we only provide the initial and the final numerical weight value.
2.4.6 Measurement interval and control interval
In our simulations, we evaluate FFsync, FFosc, Fcc, or Fvar in the time interval [t0, T] (“measurement interval”). Similarly, we enable control only in a limited time interval, which we denote by [t0C,t1C] (“control interval”), i.e., u(t)≠0 only for t∈[t0C,t1C]. Naturally, we chose t0≥t0C.
For control tasks that force a system to behave in way that cannot be maintained without control input (driving oscillations from a stationary state, Section 3.1, and synchronization (desynchronization) of desynchronized (synchronized) oscillations, Section 3.3), we set t0=t0C and T=t1C. For such tasks, no sustainable, long-lasting effect can be achieved, and we only consider the interval in which control is active.
For control tasks that initiate a switch between two stable states (Section 3.2), the algorithm might fail if there is an overlap of measurement interval and control interval as we might encounter a (local) optimum at which the control reduces the cost for t∈[t0,t1C] without initiating the intended state switch. On the other hand, the gradient computation with the adjoint method becomes less precise if there is no overlap between measurement interval and control interval. Preliminary simulations help to determine well-suited values for t0C,t1C,t0, and T. We chose a setup with no overlap between measurement interval and control interval for single-node systems (Section 3.2.1), and a setup with overlap for multi-node systems (Section 3.2.2).
2.4.7 Local optima
We use gradient descent to reach a minimum of the cost functional. Hence, the algorithm is only assured to find a local minimum in the cost landscape. We find evidence for multiple local minima (for an example, see Supplementary Section 5), but the results presented in the following are minimum-cost solutions that were repeatedly obtained with different initializations and gradient descent parameters. We will denote the control solutions therefore as the “OC”. However, we want to emphasize that there is no way to guarantee that the global minimum was found.
2.4.8 Translational invariance
If a control signal initiates a transition between two stable states, we expect the success of the transition to be invariant under shifts in time, since earlier or later control signals would lead to the same state transition. If the initial state is a stationary state (see Section 3.2.1, up-to-oscillation task), time shifts can be continuous. If the initial state is an oscillatory state (see Section 3.2.1, oscillation-to-up-task, and Section 3.2.2), time shifts can only be a multiple of the oscillation period. By shifting control signals back in time, we might find solutions with smaller FF, Fcc, or Fvar, since an earlier transition implies a longer transition time before inaccuracy is penalized, while the control cost FE remains unchanged. The results presented for the state switch in the mean-field EIF model (Section 3.2.2) are obtained by shifting the originally obtained control signal back in time. A repeated optimization does not change the results.
2.4.9 Computational complexity
Table 2 summarizes the computational complexity for computing the cost functional and its derivative computation:
The complexity for computing the synchrony and oscillation Fourier cost computation scales linearly with the number NT = T/dt of total integration time steps and with the number N of network nodes. The computation of its derivative scales quadratically with NT and linearly with N.
The computational complexity for computing the cross-correlation cost scales linearly with NT and quadratically with N. The same holds true for its derivative.
The computational complexity for computing the variance cost scales linearly with both NT and N. The same holds true for its derivative.
Computational complexity of cost functional and cost derivative computation.
FFsync~NT
FFsync~N
dFFsyncdu~NT2
dFFsyncdu~N
FFosc~NT
FFosc~N
dFFoscdu~NT2
dFFoscdu~N
Fcc~NT
Fcc~N2
dFccdu~NT
dFccdu~N2
Fvar~NT
Fvar~N
dFvardu~NT
dFvardu~N
For simulations of long time series, the Fourier cost is outperformed in terms of computational costs by the other two cost functionals. Similarly, for large networks, the cross-correlation cost is outperformed by the other two cost functionals.
3 Results3.1 Induction of oscillations
We first apply the Fourier cost defined in Equation 7 to induce oscillations in a single-node system from a constant stationary (up or down) state.
Figures 4, 5 show the computed OC and the time series of the activity variables, E and I or rE and rI, respectively, of the corresponding controlled system at the two points (A) and (B) (cf. Figures 2A, B) for the WC and the mean-field EIF models. The algorithm successfully computes a control input that creates oscillations with the target frequency in all cases. We observe that the shapes of the control signals are comparable for points (A) and (B) within a model but do not transfer across models: For the WC, the periodic OC signals for both tasks have almost vertical slopes and for each period, a broad, secondary peaks follow a sharp, initial peak but in opposite directions, u→down→osc≈-u→up→osc. In comparison, the mean-field EIF OC signals resemble a more or less distorted sine curve. These observations hold true also for other state-space locations within the up- and down-state regimes and for other target frequencies (results not shown).
OC using the Fourier cost to induce oscillations at f~=30 Hz from a stationary state in the WC model. The bottom row shows the computed OC, and the top row shows the resulting time series of the activity variables E and I. The activity of and the control inputs to the excitatory population are plotted in red, and the activity of the inhibitory population is plotted in blue. The gray-shaded regions (t < 50 and t > 350) indicate the time intervals during which the control is not active and inaccuracy is not penalized. The energy cost FE is provided in the top-left corner of the control plots. The plot on the left-hand side relates to point (A), and the plot on the right-hand side relates to point (B) (cf. Figure 2A). The weight is wF=8·104 for both points.
OC using the Fourier cost to induce oscillations at f~=30 Hz from a stationary state in the mean-field EIF model. The bottom row shows the computed OC, and the top row shows the resulting time series of the activity variables rE and rI. The activity of and the control inputs to the excitatory population are plotted in red, and the activity of the inhibitory population is plotted in blue. The gray-shaded regions (t < 50 ms and t > 350 ms) indicate the time intervals during which the control is not active and inaccuracy is not penalized. The energy cost FE is provided in the top-left corner of the control plots. The plot on the left-hand side relates to point (A), and the plot on the right-hand side relates to point (B) (cf. Figure 2B). The weight is wF = 0.555 for both points. We initialized the OC computation with an oscillatory control signal for point (A).
3.2 Switch between stable states
In this Section, we study tasks, in which the control initiates a state switch between coexisting stable states, comparing the effects of the cost functionals FF, Fcc, and Fvar.
3.2.1 Switch between stationary and oscillatory states
The mean-field EIF model exhibits a regime in state space, where a stable state of constant high activity and a stable oscillatory state coexist (see Figure 2B). We apply the Fourier cost to initiate a state switch between these stable states at point (C) in the mean-field EIF state space. For the up-to-oscillations task, we first measure the frequency fC of the oscillation at point (C) to then enforce an oscillation at this frequency. For the oscillations-to-up task, we suppress oscillations at frequency fC by setting wF < 0, thus penalizing this Fourier mode. We compute the minimum-energy control, i.e., the control with the smallest possible FE by varying wF (see Section 2.4).
Figure 6 shows the OC and the corresponding activities of the excitatory and inhibitory populations for the two control tasks. We observe that the state switch is successful for both tasks and that the control signal exhibits a short pulse close to the end of the control interval. The control input is approximately equally strong for both tasks (see Figure 6 for numerical values). For the oscillations-to-up task, FE differs by less than 2% when we vary the oscillation phase at which control is activated (results not shown).
OC to switch between the up state and the oscillatory state in the mean-field EIF model at point (C) (cf. Figure 2B). The bottom row shows the computed OC, and the top row shows the resulting activity variables rE and rI, when this control is applied. The activity of and control inputs to the excitatory population are plotted in red, and the activity of the inhibitory population is plotted in blue. The control signals are shifted back in time by two periods (see Section 2.4), and the original control signal and the corresponding original activity are plotted as a dotted line. The shift improves the Fourier cost by 0.67% (left) and 99.36% (right). The Fourier cost FF corresponding to the control shifted back in time (solid line) is provided in the top-left corner of the activity plots. The gray-shaded regions indicate the time intervals during which the control is not active (t < 100 ms and t > 500 ms) and inaccuracy is not penalized (t < 600 ms). The energy cost FE is provided in the top-right corner of the control plots. The switch from the up state to the oscillatory state is shown on the left-hand side. Here, we initialize the system in the stationary up state and set the target frequency to fC. The weight is initialized at wF = 1, 000 and decreased to wF = 0.0429 (cf. Section 2.4). The switch from the oscillatory state to the up state is shown on the right-hand side. Here, we initialize the system in the oscillatory state and suppress the frequency FC with a negative weight initialized with wF = −1 and increased to wF = −0.0502 (cf. Section 2.4). For this task, we initially enable control for 400 < t < 500 and set the control interval as shown after a few iterations. Otherwise, we only find the local OC shown in the Supplementary Section 5.
3.2.2 Switch between in-phase and out-of-phase oscillation
We consider a symmetric WC two-node network (coupling scheme sketched in Figure 3A) at point (C) in state space (see Figure 2A), where a stable IP oscillation coexists with a stable OOP oscillation. We first apply the OC algorithm using FFsync, Fcc, and Fvar to compute control signals to switch from OOP and IP. We set the target frequency to f~=fIP for the synchronization Fourier cost. We then use FFosc, Fcc, and Fvar to switch from IP to OOP. We set the target frequency to f~=fOOP for the oscillation Fourier cost, and chose negative weights wcc, wvar < 0 for cross-correlation and variance cost to penalize synchrony. We compute the minimum-energy control, i.e., the control with the smallest possible FE, by varying wF (see Section 2.4). We find that all cost functionals produce the same OC in this limit.
Figure 7 shows the results. For the switch from IP to OOP, the total required energy input FE is less than half as strong as for the switch from OOP to IP (see Figure 7 for numerical values). Similarly as in Section 3.2.1, control costs differ only marginally when varying the oscillation phase at which control is activated (results not shown).
OC switches between in-phase (IP) and out-of-phase (OOP) oscillations in the two-node WC network at point (C) (cf. Figure 2A). We use FFsync, Fcc, and Fvar for OOP → IP and FFosc, Fcc, and Fvar for IP → OOP. For any of the cost functionals, the same OC is obtained in the minimum-energy limit. The bottom row shows the computed OC, and the top row shows the resulting excitatory activity. Excitatory activity of and control inputs to node 0 (node 1) are plotted in solid red (dotted black). The control signals are shifted back in time by two periods (see Section 2.4). Control is enabled between 50 and 350 (disabled in the gray-shaded regions in the control plots), Fourier, cross-correlation, or variance costs are evaluated between 200 and 600 (ignored in the gray-shaded region in the activity plots). The switch from OOP to IP is shown on the left-hand side, and the switch from IP to OOP is shown on the right-hand side. The energy cost FE is provided in the top-left corner of the control plots. For the OOP-to-IP switch, the weights are initialized as wFsync=4000, wcc = 250, and wvar = 30000 and are dynamically reduced to wFsync=17.01, wcc = 0.8361, and wvar = 12.47. For the IP-to-OOP switch, the weights are initialized as wFosc=2000, wcc = −500, and wvar = −1000 and are dynamically adjusted to wFosc=7.284, wcc = −2.260, and wvar = −3.979.
3.3 (De-) Synchronization of oscillating networks
To study how networks can be synchronized or desynchronized, we consider a six-node WC network (coupling scheme sketched in Figure 3B), for which we find a state of asynchronous oscillation at point (D) and a state of synchronous oscillation at point (E) in state space (see Figure 2A). We apply FFsync, Fcc, and Fvar with wFsync,wcc,wvar>0 to compute the OC to synchronize (at point (D)). To desynchronize (at point (E)) the activity of the network, we set wFsync,wcc,wvar<0. We evaluate the Fourier spectrum of the network activity and use its peak frequency as f~ for the Fourier cost.
The results are shown in Figure 8. We observe that, while all three cost functionals succeed to synchronize the network at point (D), neither the synchronization Fourier cost nor the variance cost can drive the system into an asynchronous state, even though numerically, the cost contributions FFsync and Fvar lead to considerably smaller total costs. For the Fourier-controlled scenario, we analyze the spectrum of the controlled activity and find that the control slightly increases the frequency. Hence, the Fourier component of the original target frequency drops almost to zero. For the variance-controlled scenario, the control shifts the phases slightly, such that the activity of each node deviates from the network mean, while increasing the oscillation amplitude, such that the difference between node activity and mean increases. This results in a considerable increase in Fvar.
OC to enforce synchrony (left, point (D)) or asynchrony (right, point (E)) in a six-node WC network. The six colors represent one network node each. The top row shows the uncontrolled excitatory activity. Row three (five, seven) shows the computed OC, when the Fourier cost (cross-correlation cost, variance cost) is applied, and row two (four, six) the resulting excitatory activity, if this control is applied. The gray-shaded regions indicate the time intervals during which control is not active and inaccuracy is not penalized (t < 100 and t > 600). The energy costs FE are provided in the top-left corners of the control plots. FFsync, Fcc, and Fvar are provided in the top-left corner of the activity plots. The weights are wFsync=18600, wcc = 4711, and wvar = 45000 (left) and wFsync=1000, wcc = 500, and wvar = 2000 (right).
We compute the temporal mean of the Kuramoto order parameter (Acebrón et al., 2005) of the network activity to quantitatively compare the performance of the cost functionals (equations are provided in the Supplementary Section 6). The Kuramoto order parameter ranges from 0 (no synchronization) to 1 (full synchronization). The numerical values are provided in Table 3. All cost functionals succeed to improve the Kuramoto order parameter compared to the uncontrolled activity. In agreement with our observations from the activity plots (see Figure 8), the Fourier cost performs worst, and the cross-correlation cost performs best in terms of increasing (synchronization task) or decreasing (desynchronization task) the Kuramoto order parameter.
Temporal mean Kuramoto order parameter for the uncontrolled and controlled network activity for the control tasks shown in Figure 8.
Cost functional
Synchronization task
Desynchronization task
uncontrolled
0.32
0.81
FFsync
0.63
0.73
Fcc
0.99
0.33
Fvar
0.94
0.47
The values were computed for the period during which the control was active.
4 Discussion
We introduce novel cost functionals that enable to use OCT to induce oscillations or synchrony in non-linear dynamical systems without specifying a reference trajectory. We apply the cost functionals FF, Fcc, and Fvar, to two models of neural population dynamics to study different control tasks. We first enforce oscillations from stationary states using FF (see Section 3.1) and observe that OC signals drive oscillations at the given target frequency for all investigated paradigms but that the shapes of the control signals are not transferable across models. Next, we study how one can switch between stable states using FF, Fcc, and Fvar (see Section 3.2). First, we show that applying FF, one obtains OC signals that switch between a stable up state and a stable oscillatory state in a single-node mean-field EIF system (see Section 3.2.1). Second, we show that all three cost functionals can produce control signals that switch between stable states of IP and OOP oscillations in a two-node WC network. We observe that all cost functionals lead to the same OC in the minimum-energy limit. Moreover, we use the cost functionals FF, Fcc, and Fvar, to (de-) synchronize larger WC networks (see Section 3.3). We study an asynchronously oscillating six-node WC network and try to synchronize its activity using FF, Fcc, and Fvar and obtain effective control signals for all cost functionals. Finally, we study a synchronously oscillating six-node WC network and try to desynchronize its activity. Here, only the cross-correlation cost functional creates a control signal that effectively desynchronizes the network activity. Both Fourier and variance costs improve numerically, but do so by either shifting the oscillation frequency or by increasing the oscillation amplitude. In summary, we find that the Fourier cost is well suited to drive oscillations at a certain frequency. For the state-switching tasks, all three cost functionals produce satisfactory results. The same holds true for the synchronization tasks in the six-node network. For the desynchronization task, only the cross-correlation cost performs well.
Our results prove that the suggested methods enable to include oscillatory phenomena into the framework of non-linear OCT beyond (de-) synchronization of coupled phase oscillators. However, the list of cost functionals provided in this study is not exhaustive. In particular, we have evaluated further cost functionals within the framework of the adjoint method and observed inferior performances for at least one control task. A list of these discarded cost functionals is given in the Supplementary Section 7. In addition, we do not present a one-size-fits-all solution but a toolbox of methods. Depending on the control task, one needs to chose adequately from this toolbox, and one might have to make educated guesses for initializations, and on the optimization schedule (i.e., weight choices and changes, changes of the control interval or measurement interval, or shifts of the control signals). We hope that the presented examples help to pick appropriately from the toolbox and tailor the choice for the respective task. Future research may evaluate the applicability of our suggested cost functionals in other optimization schemes.
Physical neural systems are often noisy. This property is not covered here; however, the adjoint method of OCT can be extended to enable the computation of OC in noisy systems (see, e.g., Chouzouris et al., 2021). In this case, the cost functional is defined as the expected value,
F→E(F)=limM→∞1M∑m=1MF.
In numerical computations, the expected value is replaced by the mean value over a reasonably large number of realizations M,
E(F)≈1M∑m=1MF,M∈ℕ,M<∞.
This affects the computation of the adjoint state and the gradient. Preliminary studies suggest that reliable results can be obtained for the presented cost functionals also in noisy systems. The OC framework available within neurolib (Cakan et al., 2023) includes modules for OC computations in a noisy environment.
One apparent limitation of the presented Fourier cost method is that it requires to define a target oscillation frequency. Otherwise, the Fourier cost can neither drive oscillations nor synchronize a system. However, as peaks in Fourier spectra are the broader, the shorter the simulation duration (see Section 2.4), the requirement for a very precise determination of a target frequency becomes less strict. To drive synchronization, this limitation can be overcome by using other cost functionals.
Furthermore, the presented cross-correlation and variance costs do only enable phase synchrony but not frequency synchrony with two or more network nodes oscillating with fixed phase shifts. The oscillation Fourier cost FFosc is partially able to capture such cases as it can induce phase-locked oscillations in a network.
The presented methods could improve our understanding of the internal communication of neural circuits and offer new approaches to design brain stimulation protocols. Natural selection led to energy efficiency in neural communication (Quintela-López et al., 2022), and optimality principles enforced by OCT result in biologically plausible communication strategies. Hence, theoretical studies could enable conclusions on, for example, the role of single populations in controlling oscillatory patterns in the brain. Our framework might also enter research on neurological diseases in which synchrony and asynchrony play crucial roles (see, e.g., Jiruska et al., 2012; Sobayo et al., 2020; Uhlhaas and Singer, 2010; Nimmrich et al., 2015). OCT might improve therapeutic brain stimulation, as specific oscillatory brain dynamics can be targeted optimally, simultaneously reducing unintended side effects.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://github.com/lenasal/neurolib/tree/OC_osc_sync.
The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This work was supported by the DFG (German Research Foundation) via the Individual Research Grant OB 102/31.
We would like to thank our colleagues from the Neural Information Processing Group for the valuable exchange and fruitful discussions during this project.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher's note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncom.2024.1483100/full#supplementary-material
1If network nodes oscillate entirely anti-phase, the Fourier components cancel out, and an increase in power will not improve the synchronization Fourier cost.
ReferencesAcebrónJ. A.BonillaL. L.Pérez VicenteC. J.RitortF.SpiglerR. (2005). The Kuramoto model: a simple paradigm for synchronization phenomena. Rev. Mod. Phys. 77:137–185. 10.1103/RevModPhys.77.13732872818AugustinM.LadenbauerJ.BaumannF.ObermayerK. (2017). Low-dimensional spike rate models derived from networks of adaptive integrate-and-fire neurons: comparison and implementation. PLoS Comput. Biol. 13, 1–46. 10.1371/journal.pcbi.100554528644841BomelaW.SinghalB.LiJ. (2023). Engineering spatiotemporal patterns: information encoding, processing, and controllability in oscillator ensembles. Biomed. Phys. Eng. Express9:045033. 10.1088/2057-1976/ace0c937348467BretteR.GerstnerW. (2005). Adaptive exponential integrate-and-fire model as an effective description of neuronal activity. J. Neurophysiol. 94, 3637–3642. 10.1152/jn.00686.200516014787CakanC.JajcayN.ObermayerK. (2023). neurolib: a simulation framework for whole-brain neural mass modeling. Cogn Comput15:9. 10.1007/s12559-021-09931-9CakanC.ObermayerK. (2020). Biophysically grounded mean-field models of neural populations under electrical stimulation. PLoS Comput. Biol. 16, 1–30. 10.1371/journal.pcbi.100782232324734ChandrabhatlaA. S.PomeraniecI. J.HorganT. M.WatE. K.KsendzovskyA. (2023). Landscape and future directions of machine learning applications in closed-loop brain stimulation. NPJ Digit Med. 6:79. 10.1038/s41746-023-00779-x37106034ChouzourisT.RothN.CakanC.ObermayerK. (2021). Applications of optimal nonlinear control to a whole-brain network of FitzHugh-Nagumo oscillators. Phys. Rev. E. 104:024213. 10.1103/PhysRevE.104.02421334525550DasanayakeI.LiJ.-S. (2011). Optimal design of minimum-power stimuli for phase models of neuron oscillators. Phys. Rev. E. 83:061916. 10.1103/PhysRevE.83.06191621797412GrosenickL.MarshelJ. H.DeisserothK. (2015). Closed-loop and activity-guided optogenetic control. Neuron86:106–139. 10.1016/j.neuron.2015.03.03425856490HövelP.DahlemM.SchöllE. (2010). Control of synchronization in coupled neural systems by time-delayed feedback. Int. J. Bifurcat. Chaos 20. 10.1142/S0218127410026101JiruskaP.de CurtisM.JefferysJ.SchevonC.SchiffS.SchindlerK. (2012). Synchronization and desynchronization in epilepsy: controversies and hypotheses. J. Physiol. 591. 10.1113/jphysiol.2012.23959023184516KaoT.-C.HennequinG. (2019). Neuroscience out of control: control-theoretic perspectives on neural circuit dynamics. Curr. Opin. Neurobiol. 58, 122–129. 10.1016/j.conb.2019.09.00131563084KirkD. (2004). Optimal Control Theory: An Introduction. Dover Books on Electrical Engineering. New York: Dover Publications.KoryakovskiyI.KudrussM.BabukaR.CaarlsW.KirchesC.MombaurK.. (2017). Benchmarking model-free and model-based optimal control. Rob. Auton. Syst. 92, 81–90. 10.1016/j.robot.2017.02.00634398780KrylovD.DylovD. V.RosenblumM. (2020). Reinforcement learning for suppression of collective activity in oscillatory ensembles. Chaos30:033126. 10.1063/1.512890932237778LiuA.VöröslakosM.KronbergG.HeninS.KrauseM.HuangY.. (2018). Immediate neurophysiological effects of transcranial electrical stimulation. Nat. Commun. 9:5092. 10.1038/s41467-018-07233-730504921McGowanA. L.ParkesL.HeX.StanoiO.KangY.LomaxS.. (2022). Controllability of structural brain networks and the waxing and waning of negative affect in daily life. Biol. Psychiatry2, 432–439. 10.1016/j.bpsgos.2021.11.00836324655NimmrichV.DraguhnA.AxmacherN. (2015). Neuronal network oscillations in neurodegenerative diseases. Neuromolecular Med. 17, 270–284. 10.1007/s12017-015-8355-925920466ParkS.-E.LaxpatiN. G.GutekunstC.-A.ConnollyM. J.TungJ.BerglundK.. (2019). A machine learning approach to characterize the modulation of the hippocampal rhythms via optogenetic stimulation of the medial septum. Int. J. Neural Syst. 29:1950020. 10.1142/S012906571950020531505977PietrasB.DaffertshoferA. (2019). Network dynamics of coupled oscillators and phase reduction techniques. Phys. Reports819, 1–105. 10.1016/j.physrep.2019.06.001PopovychO. V.HauptmannC.TassP. A. (2006). Control of neuronal synchrony by nonlinear delayed feedback. Biol. Cybern. 95, 69–85. 10.1007/s00422-006-0066-816614837PyragasK.FedaravičiusA. P.PyragienėT.TassP. A. (2020). Entrainment of a network of interacting neurons with minimum stimulating charge. Phys. Rev. E102:012221. 10.1103/PhysRevE.102.01222132795011Quintela-LópezT.ShiinaH.AttwellD. (2022). Neuronal energy use and brain evolution. Curr. Biol. 32, R650–R655. 10.1016/j.cub.2022.02.00535728547RaoA. (2010). “A survey of numerical methods for optimal control,” in Astrodynamics 2009: Proceedings of the AAS/AIAA Astrodynamics Specialist Conference Held August 9–13 2009, Pittsburgh, Pennsylvania; Advances in the Astronautical Sciences 135 (San Diego, CA: Published for the American Astronautical Society by Univelt).SalfenmoserL.ObermayerK. (2022). Nonlinear optimal control of a mean-field model of neural population dynamics. Front. Comput. Neurosci. 16:931121. 10.3389/fncom.2022.93112135990368SalfenmoserL.ObermayerK. (2023). Optimal control of a Wilson Cowan model of neural population dynamics. Chaos33:043135. 10.1063/5.014468237097962SobayoT.FarahmandS.MogulD. J. (2020). Determining the Role of Synchrony Dynamics in Epileptic Brain Networks. Singapore: Springer Singapore, 1–28.SuppaA.HuangY.-Z.FunkeK.RiddingM.CheeranB.Di LazzaroV.. (2016). Ten years of theta burst stimulation in humans: established knowledge, unknowns and prospects. Brain Stimul. 9, 323–335. 10.1016/j.brs.2016.01.00626947241TafazoliS.MacDowellC.CheZ.LetaiK.SteinhardtC.BuschmanT. (2020). Learning to control the brain through adaptive closed-loop patterned stimulation. J. Neural Eng. 17:056007. 10.1088/1741-2552/abb86032927437TaherH.OlmiS.SchöllE. (2019). Enhancing power grid synchronization and stability through time-delayed feedback control. Phys. Rev. E. 100:062306. 10.1103/PhysRevE.100.06230631962463UhlhaasP.SingerW. (2010). Abnormal neural oscillations and synchrony in schizophrenia. Nat. Rev. Neurosci. 11:100–113. 10.1038/nrn277420087360VuM.SinghalB.ZengS.LiJ.-S. (2024). Data-driven control of oscillator networks with population-level measurement. Chaos34:033138. 10.1063/5.019185138526979WilsonD.MoehlisJ. (2014). Optimal chaotic desynchronization for neural populations. SIAM J. Appl. Dyn. Syst. 13, 276–305. 10.1137/120901702WilsonH. R.CowanJ. D. (1972). Excitatory and inhibitory interactions in localized populations of model neurons. Biophys. J. 12, 1–24. 10.1016/S0006-3495(72)86068-54332108WilsonH. R.CowanJ. D. (1973). A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue. Kybernetik13, 55–80. 10.1007/BF002887864767470YuY.XiaominW.WangQ.WangQ. (2020). A review of computational modeling and deep brain stimulation: applications to Parkinson's disease. Appl. Math Mech. 41, 1747–1768. 10.1007/s10483-020-2689-933223591ZlotnikA.LiJ.-S. (2012). Optimal entrainment of neural oscillator ensembles. J. Neural Eng. 9:046015. 10.1088/1741-2560/9/4/04601522791697ZöllerD.SandiniC.SchaerM.EliezS.BassettD. S.de VilleD. V. (2021). Structural control energy of resting-state functional brain states reveals less cost effective brain dynamics in psychosis vulnerability. Hum. Brain Mapp. 42, 2181–2200. 10.1002/hbm.2535833566395