A framework for optimal control of oscillations and synchrony applied to non-linear models of neural population dynamics

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.

1 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 Methods

2.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

$\begin{array}{ll}\overrightarrow{h}(\overrightarrow{x}(t),\dot{\overrightarrow{x}}(t),\overrightarrow{x}(t-d_1),\overrightarrow{x}(t-d_2),...,\overrightarrow{x}(t-d_{N_d}),\overrightarrow{u}(t))=0,& (1)\end{array}$

with the state vector $\overrightarrow{x}$, its time derivative $\dot{\overrightarrow{x}}$, and N_d potentially different delays d_i. In this study, we consider network systems with N network nodes and a finite simulation time. The state vector $\overrightarrow{x}$ contains the stacked state vectors ${\overrightarrow{x}}_n,n\in \left\{1,2,...N\right\}$, of all nodes. Similarly, $\overrightarrow{h}$ contains the stacked system dynamics of all network nodes. We denote the dimensionality of the state vector of a single node as N_x, such that $\left\{\overrightarrow{x},\overrightarrow{h}\right\}\in {ℝ}^{N·N_x}$. $\overrightarrow{u}(t)$ is the control input, which may effect one or more components of ${\overrightarrow{x}}_n$ for each node in a network. The OC formalism provides methods to compute a control $\overrightarrow{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

$\begin{array}{r}F=w_{\text{P}}{\displaystyle \underset{=F_{\text{P}}}{{\displaystyle \underbrace{\frac{1}{2(T-t_0)}{\displaystyle {\int }_{t_0}^T}{(\overrightarrow{x}(t)-\tilde{x}(t))}^{2}dt}}}}+{\displaystyle \underset{=F_{\text{E}}}{{\displaystyle \underbrace{\frac{1}{2}{\displaystyle {\int }_0^T}{\overrightarrow{u}}^2(t)dt}}}}\\ ={\displaystyle {\int }_0^T}f(\overrightarrow{x},\overrightarrow{u})dt.& (2)\end{array}$

The first term F_P measures the closeness of the trajectory of the state vector $\overrightarrow{x}(t)$ to a given target trajectory of state vectors $\tilde{x}(t)\in {ℝ}^{N·N_x}$ and is referred to as the precision cost. The second term F_E 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 t₀ ≤ t ≤ T, while control strength is penalized for the complete time interval. When controlling a switch between two stable states, we set t₀>0 to account for a finite transition time. The OC ${\overrightarrow{u}}^{*}$ minimizes the cost,

$\begin{array}{ll}{\overrightarrow{u}}^{*}=arg{\displaystyle \underset{\overrightarrow{u}}{min}}F.& (3)\end{array}$

The weight w_P can be tuned to trade precision against the required control strength.

In the following Section 2.2, we will consider new cost functionals $F_{\mathrm{\text{X}}}(\overrightarrow{x})$ which go beyond the precision cost $F_{\mathrm{\text{P}}}(\overrightarrow{x})$ defined in Equation 2, i.e., $F_{\mathrm{\text{P}}}(\overrightarrow{x})\to F_{\mathrm{\text{X}}}(\overrightarrow{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

$\begin{array}{ll}\frac{d{F}_{\text{X}}(\overrightarrow{x})}{d\overrightarrow{u}}={\displaystyle {\int }_0^T}\overrightarrow{g}(\overrightarrow{x})\frac{d\overrightarrow{x}}{d\overrightarrow{u}}dt,& (4)\end{array}$

where $\overrightarrow{g}(\overrightarrow{x})$ is some function of the state variable $\overrightarrow{x}$ (see Supplementary material).

The adjoint method (Kirk, 2004) enables to compute the derivative $\frac{d}{d\overrightarrow{u}}F$ of the cost functional with respect to the control. The adjoint state $\overrightarrow{\lambda }\in {ℝ}^{N·N_x}$ is defined by the differential equation

$\begin{array}{r}\dot{\overrightarrow{\lambda }}{(t)}^T\frac{\partial \overrightarrow{h}(t)}{\partial \dot{\overrightarrow{x}}}=\overrightarrow{g}(\overrightarrow{x})+{\overrightarrow{\lambda }}^T(t)\frac{\partial \overrightarrow{h}(t)}{\partial \overrightarrow{x}}\\ +{\displaystyle \sum _{i=1}^{N_d}}{\chi }_{\left[0,T-d_i\right]}{\overrightarrow{\lambda }}^T(t+d_i)\frac{\partial \overrightarrow{h}(t+d_i)}{\partial \overrightarrow{x}}-{\overrightarrow{\lambda }}^T\frac{d}{dt}\frac{\partial \overrightarrow{h}}{\partial \dot{\overrightarrow{x}}}& (5)\end{array}$

with the final condition $\overrightarrow{\lambda }(T)=0$. χ_{[ta, tb]} denotes the indicator function of the time interval [t_a, t_b]. Above expression contains vectors $\left\{\overrightarrow{\lambda },\dot{\overrightarrow{\lambda }},\overrightarrow{g}(\overrightarrow{x})\right\}\in {ℝ}^{N·N_x}$, and matrices $\left\{\frac{\partial \overrightarrow{h}}{\partial \dot{\overrightarrow{x}}},\frac{\partial \overrightarrow{h}}{\partial \overrightarrow{x}}\right\}\in {ℝ}^{(N·N_x)\times (N·N_x)}$, with ${(\frac{\partial \overrightarrow{h}}{\partial \dot{\overrightarrow{x}}})}_{ij}=\frac{\partial h_i}{\partial {ẋ}_j}$ (indices i and j denote the ith and jth vector components). Note that for the dynamical systems studied in this work, $\frac{\partial \overrightarrow{h}}{\partial \dot{\overrightarrow{x}}}$ is a constant such that the last term in Equation 5 vanishes. The adjoint state $\overrightarrow{\lambda }$ is computed by backward integration. It enables to compute the cost derivative

$\begin{array}{ll}\frac{d}{d\overrightarrow{u}}F={\displaystyle {\int }_0^T}\frac{\partial f}{\partial \overrightarrow{u}}+{\overrightarrow{\lambda }}^T\frac{\partial \overrightarrow{h}}{\partial \overrightarrow{u}}dt.& (6)\end{array}$

The integrand $\frac{\partial f}{\partial \overrightarrow{u}}+{\overrightarrow{\lambda }}^T\frac{\partial \overrightarrow{h}}{\partial \overrightarrow{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 ${\tilde{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 F_P 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 x_n(t) for network node n.

To enforce phase and frequency synchronization in a network at frequency $\tilde{f}$, we replace F_P in Equation 2 by the synchronization Fourier cost,

$\begin{array}{ll}{F}_{\text{F}}^{\text{sync}}=-\frac{1}{N^2{(T-t_0)}^2}|{\displaystyle {\int }_{t_0}^T{\Sigma }_x}(t)·e^{-i\omega t} dt{|}^2& (7)\end{array}$

with ${\Sigma }_x(t)=\sum _{n=1}^N{x}_n(t)$ and $\omega =2\pi \tilde{f}$. $F_{\mathrm{\text{F}}}^{\mathrm{\text{sync}}}$ is the squared Fourier component corresponding to the frequency $\tilde{f}$ of the sum of observables over all N network nodes. The computation of the adjoint state using Equation 5 requires

$\begin{array}{ll}{\left(g_{\text{F}}^{\text{sync}}(\overrightarrow{x})\right)}_k=-\frac{2}{N^2{(T-t_0)}^2}{\displaystyle {\int }_{t_0}^T}{\Sigma }_x(\tau )·cos(\omega (t-\tau ))d\tau .& (8)\end{array}$

The derivation of Equation 8 is provided in the Supplementary Section 2.1.

To enforce frequency synchronization in a network at frequency $\tilde{f}$ irrespective of the relative phases of its nodes, we replace F_P in Equation 2 by the oscillation Fourier cost

$\begin{array}{ll}{F}_{\text{F}}^{\text{osc}}=-\frac{1}{N{(T-t_0)}^2}{\displaystyle \sum _{n=1}^N|}{\displaystyle {\int }_{t_0}^T{x}_n}(t)·e^{-i\omega t} dt{|}^2.& (9)\end{array}$

$F_{\mathrm{\text{F}}}^{\mathrm{\text{osc}}}$ is the sum over all N network nodes of the node-wise squared Fourier component corresponding to the frequency $\tilde{f}$. The computation of the adjoint state using Equation 5 requires

$\begin{array}{ll}{\left(g_{\text{F}}^{\text{osc}}(\overrightarrow{x})\right)}_k=-\frac{2}{N{(T-t_0)}^2}{\displaystyle {\int }_{t_0}^T}{x}_k(\tau )·cos(\omega (t-\tau ))d\tau .& (10)\end{array}$

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 $\tilde{f}$. We will denote the cost by F_F in this case. The oscillation Fourier cost $F_{\mathrm{\text{F}}}^{\mathrm{\text{osc}}}$ improves, if the power corresponding to the target frequency $\tilde{f}$ increases in any network node, irrespective of the relative phases. The synchronization Fourier cost $F_{\mathrm{\text{F}}}^{\mathrm{\text{sync}}}$ improves, either if oscillatory network activity at frequency $\tilde{f}$ synchronizes, or if the power in the network nodes corresponding to $\tilde{f}$ increases¹.

To enforce phase synchronization in a network irrespective of the oscillation frequency, we replace F_P in Equation 2 by either the cross-correlation cost or the variance cost. The cross-correlation cost (Chouzouris et al., 2021) is defined as

$\begin{array}{ll}{F}_{\text{cc}}=-\frac{2}{N(N-1)(T-t_0)}{\displaystyle {\int }_{t_0}^T}{\displaystyle \sum _{\begin{array}{c}n=1\\ l=n+1\end{array}}^N}\frac{(x_n(t)-{\overline{x}}_n)(x_l(t)-{\overline{x}}_l)}{\sigma (x_n)\sigma (x_l)}dt,& (11)\end{array}$

where $\overline{x}=\frac{1}{T-t_0}{\int }_{t_0}^{T}x(t)dt$ is the temporal mean of x(t), and ${\sigma }^2(x)=\frac{1}{T-t_0}{\int }_{t_0}^T{(x(t)-\overline{x})}^{2}dt$ is its variance. The computation of the adjoint state using Equation 5 requires

$\begin{array}{ll}\begin{array}{r}\hfill {\left(g_{\text{cc}}(\overrightarrow{x})\right)}_k=\frac{2}{N(N-1)(T-t_0)}{\displaystyle \sum _{\begin{array}{c}l=1\\ l\ne k\end{array}}^N\frac{(x_l(t)-{\overline{x}}_l)}{\sigma (x_k)\sigma (x_l)}}\\ \hfill -\frac{(x_k(t)-{\overline{x}}_k){ℐ}_{kl}}{(T-t_0)\sigma {(x_k)}^3\sigma (x_l)},\end{array}& (12)\end{array}$

with $I_{kl}={\int }_{t_0}^T(x_k(t)-{\overline{x}}_k)(x_l(t)-{\overline{x}}_l)dt$. The derivation of Equation 12 is provided in the Supplementary Section 2.3.

The variance cost is defined as

$\begin{array}{ll}{F}_{\text{var}}=\frac{1}{N(T-t_0)}{\displaystyle {\int }_{t_0}^T}{\displaystyle \sum _{n=1}^N}{(x_n(t)-\overline{x}(t))}^{2}dt,& (13)\end{array}$

where $\overline{x}(t)=\frac{1}{N}\sum _{n=1}^N{x}_n(t)$ is the network mean of x(t). The computation of the adjoint state using Equation 5 requires

$\begin{array}{ll}{\left(g_{\text{var}}(\overrightarrow{x})\right)}_k=\frac{2(x_k(t)-\overline{x}(t))}{N(T-t_0)}.& (14)\end{array}$

The derivation of Equation 14 is provided in the Supplementary Section 2.4.

Both F_cc and F_var are defined for networks with N≥2 and cannot be applied to a single-node system.

2.2.1 Control weights

When replacing F_P in Equation 2 by one of the abovementioned cost functionals, we also introduce the respective weights, $w_{\mathrm{\text{F}}}^{\mathrm{\text{sync}}}$, $w_{\mathrm{\text{F}}}^{\mathrm{\text{osc}}}$, w_cc, or w_var replacing w_P 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

$\begin{array}{ll}\begin{array}{l}{\tau }_E{\dot{E}}_n(t)=-E_n(t)+(1-E_n(t))·S(c_{EE}{E}_n(t)\hfill \\ -c_{EI}{I}_n(t)+E^{\text{ext}}+u_n(t)+c_{\text{gl}}{\displaystyle \sum _{m=1}^N{C}_{nm}}{E}_m(t-D_{nm}))\hfill \\ {\tau }_I{\dot{I}}_n(t)=-I_n(t)+(1-I_n(t))·S(c_{IE}{E}_n(t)-c_{II}{I}_n(t)+I^{\text{ext}}).\hfill \end{array}& (15)\end{array}$

The activity variables E_n(t) and I_n(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}. E^ext and I^ext are static external inputs that serve as order parameters for the network. u_n(t) is the time-dependent control, c_gl 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−μ))⁻¹ computes a synaptic input to E_n or I_n 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.

Table 1

Table 1. Parameters of the WC model.

Figure 1

Figure 1. Cartoon of a single node of the WC network model. Static external inputs are indicated by E^ext and I^ext, network inputs by E^nw, 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 E^ext and I^ext, 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).

Figure 2

Figure 2. (A) Slice of state space of one node of the WC model. The horizontal (vertical) axis corresponds to the external input E^ext (I^ext). 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 ${\mu }_E^{\mathrm{\text{ext}}}$ (${\mu }_I^{\mathrm{\text{ext}}}$). 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 r_E and r_I of the excitatory and inhibitory populations can be obtained as the activity state variables. r_E and r_I 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 ${\mu }_E^{\mathrm{\text{ext}}}$ and ${\mu }_I^{\mathrm{\text{ext}}}$.

Figure 2B shows a slice of state space of the one-node mean-field EIF model. Depending on the external inputs ${\mu }_E^{\mathrm{\text{ext}}}$ and ${\mu }_I^{\mathrm{\text{ext}}}$, 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

$\begin{array}{ll}C=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),D=\left(\begin{array}{cc}0& 9.5\\ 9.5& 0\end{array}\right).& (16)\end{array}$

The coupling scheme is sketched in Figure 3A. We set c_gl = 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 $f_{\mathrm{\text{IP}}}^{-1}=13.89$ coexists with a stable out-of-phase oscillation (OOP) with period $f_{\mathrm{\text{OOP}}}^{-1}=22.72$.

Figure 3

Figure 3. (A) Cartoon of the neural mass model with two WC nodes studied in Section 3.2.2. Static external inputs are indicated by E^ext and I^ext and control inputs by u_0,1 (t). The nodes are coupled via the excitatory population with a global coupling strength c_gl. (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

$\begin{array}{ll}C=\left(\begin{array}{cccccc}0& 1& 0& 0& 0& 1\\ 1& 0& 1& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\\ 0& 0& 1& 0& 0& 1\\ 1& 1& 0& 1& 0& 1\\ 0& 1& 1& 0& 0& 0\end{array}\right),D=\left(\begin{array}{cccccc}0& 12& 0& 0& 0& 8\\ 8& 0& 13& 0& 1& 0\\ 0& 0& 0& 0& 0& 9\\ 0& 0& 4& 0& 0& 11\\ 5& 17& 0& 14& 0& 18\\ 0& 0& 3& 0& 0& 0\end{array}\right).& (17)\end{array}$

The coupling scheme is sketched in Figure 3B. We set c_gl = 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 simulations

2.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 $\tilde{f}$ compared to the natural frequency of the target state.

2.4.4 Initialization

As an initial control guess, we use u₀ = 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 $\frac{\partial f}{\partial \overrightarrow{u}}$ is computed. This is due to the fact that the activity, the excitatory firing rate r_E(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, r_E(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 u₀ which oscillates at the target frequency $\tilde{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 w_F, w_cc, or w_var 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, ${\overrightarrow{u}}^{*}(t)=0$. Any finite input would increase the total cost via F_E more than it might improve F_F, F_cc, or F_var. 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 w_F, w_cc, or w_var, 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 $F_{\mathrm{\text{F}}}^{\mathrm{\text{sync}}}$, $F_{\mathrm{\text{F}}}^{\mathrm{\text{osc}}}$, F_cc, or F_var in the time interval [t₀, T] (“measurement interval”). Similarly, we enable control only in a limited time interval, which we denote by $\left[t_0^{\mathrm{\text{C}}},t_1^{\mathrm{\text{C}}}\right]$ (“control interval”), i.e., u(t)≠0 only for $t\in \left[t_0^{\mathrm{\text{C}}},t_1^{\mathrm{\text{C}}}\right]$. Naturally, we chose $t_0\ge t_0^{\mathrm{\text{C}}}$.

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 $t_0=t_0^{\mathrm{\text{C}}}$ and $T=t_1^{\mathrm{\text{C}}}$. 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\in \left[t_0,t_1^{\mathrm{\text{C}}}\right]$ 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 $t_0^{\mathrm{\text{C}}},t_1^{\mathrm{\text{C}}},t_0$, 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 F_F, F_cc, or F_var, since an earlier transition implies a longer transition time before inaccuracy is penalized, while the control cost F_E 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 N_T = T/dt of total integration time steps and with the number N of network nodes. The computation of its derivative scales quadratically with N_T and linearly with N.

• The computational complexity for computing the cross-correlation cost scales linearly with N_T and quadratically with N. The same holds true for its derivative.

• The computational complexity for computing the variance cost scales linearly with both N_T and N. The same holds true for its derivative.

Table 2

Table 2. Computational complexity of cost functional and cost derivative computation.

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 Results

3.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 r_E and r_I, 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, ${\overrightarrow{u}}^{\mathrm{\text{down}}\to \mathrm{\text{osc}}}\approx -{\overrightarrow{u}}^{\mathrm{\text{up}}\to \mathrm{\text{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).

Figure 4

Figure 4. OC using the Fourier cost to induce oscillations at $\tilde{f}=30\text{ 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 F_E 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 $w_{\mathrm{\text{F}}}=8·10^4$ for both points.

Figure 5

Figure 5. OC using the Fourier cost to induce oscillations at $\tilde{f}=30\text{ 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 r_E and r_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 ms and t > 350 ms) indicate the time intervals during which the control is not active and inaccuracy is not penalized. The energy cost F_E 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 w_F = 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 F_F, F_cc, and F_var.

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 f_C 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 f_C by setting w_F < 0, thus penalizing this Fourier mode. We compute the minimum-energy control, i.e., the control with the smallest possible F_E by varying w_F (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, F_E differs by less than 2% when we vary the oscillation phase at which control is activated (results not shown).

Figure 6

Figure 6. 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 r_E and r_I, 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 F_F 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 F_E 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 f_C. The weight is initialized at w_F = 1, 000 and decreased to w_F = 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 F_C with a negative weight initialized with w_F = −1 and increased to w_F = −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 $F_{\mathrm{\text{F}}}^{\mathrm{\text{sync}}}$, F_cc, and F_var to compute control signals to switch from OOP and IP. We set the target frequency to $\tilde{f}=f_{\mathrm{\text{IP}}}$ for the synchronization Fourier cost. We then use $F_{\mathrm{\text{F}}}^{\mathrm{\text{osc}}}$, F_cc, and F_var to switch from IP to OOP. We set the target frequency to $\tilde{f}=f_{\mathrm{\text{OOP}}}$ for the oscillation Fourier cost, and chose negative weights w_cc, w_var < 0 for cross-correlation and variance cost to penalize synchrony. We compute the minimum-energy control, i.e., the control with the smallest possible F_E, by varying w_F (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 F_E 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).

Figure 7

Figure 7. 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 $F_{\mathrm{\text{F}}}^{\mathrm{\text{sync}}}$, F_cc, and F_var for OOP → IP and $F_{\mathrm{\text{F}}}^{\mathrm{\text{osc}}}$, F_cc, and F_var 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 F_E is provided in the top-left corner of the control plots. For the OOP-to-IP switch, the weights are initialized as $w_{\mathrm{\text{F}}}^{\mathrm{\text{sync}}}=4000$, w_cc = 250, and w_var = 30000 and are dynamically reduced to $w_{\mathrm{\text{F}}}^{\mathrm{\text{sync}}}=17.01$, w_cc = 0.8361, and w_var = 12.47. For the IP-to-OOP switch, the weights are initialized as $w_{\mathrm{\text{F}}}^{\mathrm{\text{osc}}}=2000$, w_cc = −500, and w_var = −1000 and are dynamically adjusted to $w_{\mathrm{\text{F}}}^{\mathrm{\text{osc}}}=7.284$, w_cc = −2.260, and w_var = −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 $F_{\mathrm{\text{F}}}^{\mathrm{\text{sync}}}$, F_cc, and F_var with $w_{\mathrm{\text{F}}}^{\mathrm{\text{sync}}},w_{\mathrm{\text{cc}}},w_{\mathrm{\text{var}}}>0$ to compute the OC to synchronize (at point (D)). To desynchronize (at point (E)) the activity of the network, we set $w_{\mathrm{\text{F}}}^{\mathrm{\text{sync}}},w_{\mathrm{\text{cc}}},w_{\mathrm{\text{var}}}<0$. We evaluate the Fourier spectrum of the network activity and use its peak frequency as $\tilde{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 $F_{\mathrm{\text{F}}}^{\mathrm{\text{sync}}}$ and F_var 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 F_var.

Figure 8

Figure 8. 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 F_E are provided in the top-left corners of the control plots. $F_{\mathrm{\text{F}}}^{\mathrm{\text{sync}}}$, F_cc, and F_var are provided in the top-left corner of the activity plots. The weights are $w_{\mathrm{\text{F}}}^{\mathrm{\text{sync}}}=18600$, w_cc = 4711, and w_var = 45000 (left) and $w_{\mathrm{\text{F}}}^{\mathrm{\text{sync}}}=1000$, w_cc = 500, and w_var = 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.

Table 3

Table 3. Temporal mean Kuramoto order parameter for the uncontrolled and controlled network activity for the control tasks shown in Figure 8.

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 F_F, F_cc, and F_var, to two models of neural population dynamics to study different control tasks. We first enforce oscillations from stationary states using F_F (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 F_F, F_cc, and F_var (see Section 3.2). First, we show that applying F_F, 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 F_F, F_cc, and F_var, 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 F_F, F_cc, and F_var 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,

$\begin{array}{ll}F\to E(F)={\displaystyle \underset{M\to \infty }{lim}}\frac{1}{M}{\displaystyle \sum _{m=1}^M}F.& (18)\end{array}$

In numerical computations, the expected value is replaced by the mean value over a reasonably large number of realizations M,

$\begin{array}{ll}E(F)\approx \frac{1}{M}{\displaystyle \sum _{m=1}^M}F,M\in \mathbb{N},M<\infty .& (19)\end{array}$

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 $F_{\mathrm{\text{F}}}^{\mathrm{\text{osc}}}$ 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.

Author contributions

LS: Validation, Conceptualization, Writing – review & editing, Writing – original draft, Software, Methodology, Investigation, Data curation. KO: Writing – review & editing, Supervision, Funding acquisition, Conceptualization.

Funding

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.

Acknowledgments

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

Footnotes

1. ^If network nodes oscillate entirely anti-phase, the Fourier components cancel out, and an increase in power will not improve the synchronization Fourier cost.

References

Acebrón, J. A., Bonilla, L. L., Pérez Vicente, C. J., Ritort, F., and Spigler, R. (2005). The Kuramoto model: a simple paradigm for synchronization phenomena. Rev. Mod. Phys. 77:137–185. doi: 10.1103/RevModPhys.77.137

PubMed Abstract | Crossref Full Text | Google Scholar

Augustin, M., Ladenbauer, J., Baumann, F., and Obermayer, K. (2017). Low-dimensional spike rate models derived from networks of adaptive integrate-and-fire neurons: comparison and implementation. PLoS Comput. Biol. 13, 1–46. doi: 10.1371/journal.pcbi.1005545

PubMed Abstract | Crossref Full Text | Google Scholar

Bomela, W., Singhal, B., and Li, J. (2023). Engineering spatiotemporal patterns: information encoding, processing, and controllability in oscillator ensembles. Biomed. Phys. Eng. Express 9:045033. doi: 10.1088/2057-1976/ace0c9

PubMed Abstract | Crossref Full Text | Google Scholar

Brette, R., and Gerstner, W. (2005). Adaptive exponential integrate-and-fire model as an effective description of neuronal activity. J. Neurophysiol. 94, 3637–3642. doi: 10.1152/jn.00686.2005

PubMed Abstract | Crossref Full Text | Google Scholar

Cakan, C., Jajcay, N., and Obermayer, K. (2023). neurolib: a simulation framework for whole-brain neural mass modeling. Cogn Comput 15:9. doi: 10.1007/s12559-021-09931-9

Crossref Full Text | Google Scholar

Cakan, C., and Obermayer, K. (2020). Biophysically grounded mean-field models of neural populations under electrical stimulation. PLoS Comput. Biol. 16, 1–30. doi: 10.1371/journal.pcbi.1007822

PubMed Abstract | Crossref Full Text | Google Scholar

Chandrabhatla, A. S., Pomeraniec, I. J., Horgan, T. M., Wat, E. K., and Ksendzovsky, A. (2023). Landscape and future directions of machine learning applications in closed-loop brain stimulation. NPJ Digit Med. 6:79. doi: 10.1038/s41746-023-00779-x

PubMed Abstract | Crossref Full Text | Google Scholar

Chouzouris, T., Roth, N., Cakan, C., and Obermayer, K. (2021). Applications of optimal nonlinear control to a whole-brain network of FitzHugh-Nagumo oscillators. Phys. Rev. E. 104:024213. doi: 10.1103/PhysRevE.104.024213

PubMed Abstract | Crossref Full Text | Google Scholar

Dasanayake, I., and Li, J.-S. (2011). Optimal design of minimum-power stimuli for phase models of neuron oscillators. Phys. Rev. E. 83:061916. doi: 10.1103/PhysRevE.83.061916

PubMed Abstract | Crossref Full Text | Google Scholar

Grosenick, L., Marshel, J. H., and Deisseroth, K. (2015). Closed-loop and activity-guided optogenetic control. Neuron 86:106–139. doi: 10.1016/j.neuron.2015.03.034

PubMed Abstract | Crossref Full Text | Google Scholar

Hövel, P., Dahlem, M., and Schöll, E. (2010). Control of synchronization in coupled neural systems by time-delayed feedback. Int. J. Bifurcat. Chaos 20. doi: 10.1142/S0218127410026101

Crossref Full Text | Google Scholar

Jiruska, P., de Curtis, M., Jefferys, J., Schevon, C., Schiff, S., and Schindler, K. (2012). Synchronization and desynchronization in epilepsy: controversies and hypotheses. J. Physiol. 591. doi: 10.1113/jphysiol.2012.239590

PubMed Abstract | Crossref Full Text | Google Scholar

Kao, T.-C., and Hennequin, G. (2019). Neuroscience out of control: control-theoretic perspectives on neural circuit dynamics. Curr. Opin. Neurobiol. 58, 122–129. doi: 10.1016/j.conb.2019.09.001

PubMed Abstract | Crossref Full Text | Google Scholar

Kirk, D. (2004). Optimal Control Theory: An Introduction. Dover Books on Electrical Engineering. New York: Dover Publications.

Google Scholar

Koryakovskiy, I., Kudruss, M., Babuka, R., Caarls, W., Kirches, C., Mombaur, K., et al. (2017). Benchmarking model-free and model-based optimal control. Rob. Auton. Syst. 92, 81–90. doi: 10.1016/j.robot.2017.02.006

PubMed Abstract | Crossref Full Text | Google Scholar

Krylov, D., Dylov, D. V., and Rosenblum, M. (2020). Reinforcement learning for suppression of collective activity in oscillatory ensembles. Chaos 30:033126. doi: 10.1063/1.5128909

PubMed Abstract | Crossref Full Text | Google Scholar

Liu, A., Vöröslakos, M., Kronberg, G., Henin, S., Krause, M., Huang, Y., et al. (2018). Immediate neurophysiological effects of transcranial electrical stimulation. Nat. Commun. 9:5092. doi: 10.1038/s41467-018-07233-7

PubMed Abstract | Crossref Full Text | Google Scholar

McGowan, A. L., Parkes, L., He, X., Stanoi, O., Kang, Y., Lomax, S., et al. (2022). Controllability of structural brain networks and the waxing and waning of negative affect in daily life. Biol. Psychiatry 2, 432–439. doi: 10.1016/j.bpsgos.2021.11.008

PubMed Abstract | Crossref Full Text | Google Scholar

Nimmrich, V., Draguhn, A., and Axmacher, N. (2015). Neuronal network oscillations in neurodegenerative diseases. Neuromolecular Med. 17, 270–284. doi: 10.1007/s12017-015-8355-9

PubMed Abstract | Crossref Full Text | Google Scholar

Park, S.-E., Laxpati, N. G., Gutekunst, C.-A., Connolly, M. J., Tung, J., Berglund, K., et al. (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. doi: 10.1142/S0129065719500205

PubMed Abstract | Crossref Full Text | Google Scholar

Pietras, B., and Daffertshofer, A. (2019). Network dynamics of coupled oscillators and phase reduction techniques. Phys. Reports 819, 1–105. doi: 10.1016/j.physrep.2019.06.001

Crossref Full Text | Google Scholar

Popovych, O. V., Hauptmann, C., and Tass, P. A. (2006). Control of neuronal synchrony by nonlinear delayed feedback. Biol. Cybern. 95, 69–85. doi: 10.1007/s00422-006-0066-8

PubMed Abstract | Crossref Full Text | Google Scholar

Pyragas, K., Fedaravičius, A. P., Pyragienė, T., and Tass, P. A. (2020). Entrainment of a network of interacting neurons with minimum stimulating charge. Phys. Rev. E 102:012221. doi: 10.1103/PhysRevE.102.012221

PubMed Abstract | Crossref Full Text | Google Scholar

Quintela-López, T., Shiina, H., and Attwell, D. (2022). Neuronal energy use and brain evolution. Curr. Biol. 32, R650–R655. doi: 10.1016/j.cub.2022.02.005

PubMed Abstract | Crossref Full Text | Google Scholar

Rao, A. (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).

Google Scholar

Salfenmoser, L., and Obermayer, K. (2022). Nonlinear optimal control of a mean-field model of neural population dynamics. Front. Comput. Neurosci. 16:931121. doi: 10.3389/fncom.2022.931121

PubMed Abstract | Crossref Full Text | Google Scholar

Salfenmoser, L., and Obermayer, K. (2023). Optimal control of a Wilson Cowan model of neural population dynamics. Chaos 33:043135. doi: 10.1063/5.0144682

PubMed Abstract | Crossref Full Text | Google Scholar

Sobayo, T., Farahmand, S., and Mogul, D. J. (2020). Determining the Role of Synchrony Dynamics in Epileptic Brain Networks. Singapore: Springer Singapore, 1–28.

Google Scholar

Suppa, A., Huang, Y.-Z., Funke, K., Ridding, M., Cheeran, B., Di Lazzaro, V., et al. (2016). Ten years of theta burst stimulation in humans: established knowledge, unknowns and prospects. Brain Stimul. 9, 323–335. doi: 10.1016/j.brs.2016.01.006

PubMed Abstract | Crossref Full Text | Google Scholar

Tafazoli, S., MacDowell, C., Che, Z., Letai, K., Steinhardt, C., and Buschman, T. (2020). Learning to control the brain through adaptive closed-loop patterned stimulation. J. Neural Eng. 17:056007. doi: 10.1088/1741-2552/abb860

PubMed Abstract | Crossref Full Text | Google Scholar

Taher, H., Olmi, S., and Schöll, E. (2019). Enhancing power grid synchronization and stability through time-delayed feedback control. Phys. Rev. E. 100:062306. doi: 10.1103/PhysRevE.100.062306

PubMed Abstract | Crossref Full Text | Google Scholar

Uhlhaas, P., and Singer, W. (2010). Abnormal neural oscillations and synchrony in schizophrenia. Nat. Rev. Neurosci. 11:100–113. doi: 10.1038/nrn2774

PubMed Abstract | Crossref Full Text | Google Scholar

Vu, M., Singhal, B., Zeng, S., and Li, J.-S. (2024). Data-driven control of oscillator networks with population-level measurement. Chaos 34:033138. doi: 10.1063/5.0191851

PubMed Abstract | Crossref Full Text | Google Scholar

Wilson, D., and Moehlis, J. (2014). Optimal chaotic desynchronization for neural populations. SIAM J. Appl. Dyn. Syst. 13, 276–305. doi: 10.1137/120901702

Crossref Full Text | Google Scholar

Wilson, H. R., and Cowan, J. D. (1972). Excitatory and inhibitory interactions in localized populations of model neurons. Biophys. J. 12, 1–24. doi: 10.1016/S0006-3495(72)86068-5

PubMed Abstract | Crossref Full Text | Google Scholar

Wilson, H. R., and Cowan, J. D. (1973). A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue. Kybernetik 13, 55–80. doi: 10.1007/BF00288786

PubMed Abstract | Crossref Full Text | Google Scholar

Yu, Y., Xiaomin, W., Wang, Q., and Wang, Q. (2020). A review of computational modeling and deep brain stimulation: applications to Parkinson's disease. Appl. Math Mech. 41, 1747–1768. doi: 10.1007/s10483-020-2689-9

PubMed Abstract | Crossref Full Text | Google Scholar

Zlotnik, A., and Li, J.-S. (2012). Optimal entrainment of neural oscillator ensembles. J. Neural Eng. 9:046015. doi: 10.1088/1741-2560/9/4/046015

PubMed Abstract | Crossref Full Text | Google Scholar

Zöller, D., Sandini, C., Schaer, M., Eliez, S., Bassett, D. S., and de Ville, D. 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. doi: 10.1002/hbm.25358

PubMed Abstract | Crossref Full Text | Google Scholar