Control of coupled neural oscillations using near-periodic inputs

Deep brain stimulation is a well-established treatment for medication resistant Parkinson’s disease.^4,37 DBS has also shown promise as a treatment for other neurological disorders, such as epilepsy,²⁵ depression,³¹ tinnitus,⁵⁶ and Tourette syndrome.⁵³ While the fundamental mechanisms of DBS are not fully understood, one prevalent hypothesis postulates that DBS disrupts pathological synchronization among neurons in the basal ganglia-cortical loop associated with the motor symptoms of Parkinson’s disease.^5,17,20,24 This hypothesis has led to the development of a variety of electrical stimulation techniques designed to desynchronize populations of neurons. Coordinated reset^28,57 is one such strategy, whereby high frequency stimulation is coordinated at different sites to desynchronize a population of oscillators that are distributed in space. The coordinated reset strategy has been shown to produce long lasting therapeutic effects in both animal and human trials.^1,60 Other methods for achieving neural desynchronization have focused on driving neurons to a phaseless set so that they desynchronize during the relaxation to their nominal limit cycles,³⁶ identifying energy optimal stimuli that achieve positive Lyapunov exponents and produce chaotic desynchronization,^64,71 controlling neurons to a prespecified phase distribution,³⁴ and separating neurons into distinct clusters.^12,29,72

In a clinical setting, DBS is usually applied by injecting high frequency (130–180 Hz) pulse trains of electrical current into an appropriate brain region.^21,59 While DBS is administered in an open-loop fashion in most cases, there has been a growing interest in developing demand-controlled and adaptive methods for applying DBS both to increase therapeutic efficacy and reduce the emergence of side effects. For instance, measured local field potentials (LFPs) have been used as a control variable with the goal of keeping the LFP power spectrum close to that of the tremor-free condition.^32,44,46 While most adaptive DBS (aDBS) strategies are fairly simple from a control theoretic perspective, generally turning off or on depending on the value of the chosen biomarker relative to some threshold, other experimental^7,18,49 and computational studies^3,10,19 have shown that timing the application of DBS input relative to the phase of a patient’s tremor rhythm can also be an effective closed-loop stimulation strategy. Indeed, phase-amplitude reduction strategies can be used to predict how phase-triggered stimulation will influence the amplitude of pathological oscillations⁶³ that are thought to contribute to the symptoms of Parkinson’s disease; this information can subsequently be used to develop stimulation protocols for suppressing oscillation amplitude.^9,62 Additionally, slow modulation of the magnitude of the DBS input has been shown to be an effective aDBS strategy. For instance, increasing or decreasing the strength of the DBS pulses according to the beta band LFP power⁴⁶ or the severity of tremor²⁷ has been successfully demonstrated in human trials. Computational studies have also suggested that the magnitude of the DBS input could be controlled using a variety of delayed-feedback approaches;^41–43 related delayed-feedback approaches consider non-pulsatile inputs.^39,47

Motivated by these preliminary aDBS control strategies, this work focuses on understanding the dynamical behavior of a large population of coupled neurons in response to a nearly periodic input. By nearly periodic, we mean an input that is nominally periodic (for example, as is the case for conventional DBS) with the addition of a slowly varying amplitude and phase offset that can be used to coordinate the timing of an applied input with a specified phase of the population oscillation. As illustrated in the results to follow, these additional degrees of freedom allow for the development of a closed-loop control strategy that can be used to manipulate the phase cohesion (i.e., synchronization) of a large population of neural oscillators as measured by the Kuramoto order parameter. Additionally, we identify and validate a strategy for estimating the necessary terms of the proposed closed-loop control framework in real-time so that it can be implemented without knowledge of an underlying dynamical model.

The organization of this paper is as follows: Sec. II provides necessary background on the phase-amplitude reduction frameworks used for the theoretical analysis in this work. Section III derives a reduced order model to characterize the behavior of a large population of noisy, coupled, conductance-based neural oscillators. Here, we employ the recently developed adaptive phase-amplitude reduction framework^67,68 in conjunction with formal averaging techniques^16,51 to arrive at a low order, nonlinear model describing population-level behaviors of the neural oscillators. Section IV describes a closed-loop control strategy that can be used to manipulate the phase cohesion and disrupt synchronization in population-level oscillations. Additionally, we propose a data-based strategy to estimate the terms of the associated reduced order model in real-time. Theoretical results are validated by numerical simulations with results presented in Sec. V for a Fokker–Planck representation of the dynamics of a coupled oscillator population and for a full model of conductance-based neurons. Section VI provides concluding remarks.

Our primary focus in this work is on the analysis of coupled, noisy, conductance-based neural oscillators in response to near-periodic forcing. The general model is described below:

Above, $N$ sets the number of neurons in the population, $C$ is a constant membrane capacitance, $Vi$ is the transmembrane voltage of neuron $i$ (in mV), $ai∈RNa$ is a collection of auxiliary variables (gating variables, ionic concentrations, etc.) with dynamics governed by $fa$⁠, $Iion$ is a collection of ionic currents, $si$ is a synaptic variable where $c1,VT,σT$⁠, and $c2$ determine the dynamics of the synaptic variable, $Esyn=−100$ mV is the reversal potential of the neurotransmitter so that the coupling is inhibitory, $τ$ is a constant time delay, $g0>0$ is a constant conductance that sets the synaptic coupling strength, $Istim(t)$ is an externally applied current that is applied identically to each neuron, and $ηi(t)$ is a zero-mean white noise process associated with neuron $i$ for which $D=1μA2/cm4$ sets the intensity. Here, the collection $η1(t),…,ηN(t)$ is independent and identically distributed. The specific forms of $Iion$ and $fa$ are model dependent. All neurons in Eq. (1) are assumed to be identical. Because DBS is often applied to the subthalamic nucleus,¹⁷ for concreteness, we will use a model of thalamic neurons from⁵⁰ to specify the ionic dynamics. These model equations are given in the  Appendix. Terms governing the synaptic variables are taken to be $c1=3,VT=−20mV,σT=0.8mV,c2=1$⁠. Inhibitory synaptic coupling is chosen because it synchronizes the population in the absence of input. We emphasize that due to the lack of heterogeneity between neurons, the simple all-to-all coupling structure, and the lack of coupling between upstream and downstream neural populations, Eq. (1) is not meant to represent a specific target brain region for DBS. Rather, our goal is to address challenges associated with desynchronizing a population of pathologically synchronized neurons from both a dynamical systems and control theoretic perspective.

In this work, we will adopt the hypothesis that pathological synchronization among neurons in the basal ganglia-cortical loop of the brain contributes to the motor symptoms of Parkinson’s disease,^5,20,24 and that therapeutic DBS disrupts this synchronization. Consequently, we will consider the identification of control strategies using $Istim(t)$ to desynchronize a nominally synchronized population of neurons from Eq. (1). Our focus in this work is on the study of the overall population response of (1) to nearly periodic inputs $Istim(t)$⁠; the equations governing $Istim$ will take the general form

where $fs(t)$ is a $Ts$-periodic input with amplitude $μ$ and phase offset $ζ$⁠. Both $ζ$ and $μ$ are allowed to change slowly according to the relations $fζ$ and $fμ$⁠, respectively. The input $Istim(t)$ is nearly $Ts$-periodic, where the magnitude and phase shift are allowed to slowly vary in time. Further structure for the terms $fζ$ and $fμ$ will be considered in the analysis to follow.

Figure 1 shows the behavior of a population of $N=1000$ synaptically coupled neurons from Eq. (1) taking the synaptic coupling time delay to be $τ=4$ ms. Initially, the coupling strength $g0$ is set to zero allowing noise to fully desynchronize the population. At $t=100$ ms, the coupling strength is set to $g0=0.35mS/cm2$ yielding synchronized oscillations with a period of approximately 10.9 ms. During this simulation, panel (a) shows the transmembrane voltage of a representative sample of neurons (colored traces) along with the average voltage (black line). Panel (b) shows the average value of time-delayed synaptic variable,

The average synaptic variable will be used as a model observable in the development of the control strategies in Secs. III–V. In Fig. 1, the emergence of synchronization can clearly be seen from the individual neural traces in panel (a). In order to quantify the level of synchronization, we will consider the Kuramoto order parameter

where $θj$ is the asymptotic phase of neuron $j$ [as defined by the isochrons according to Eq. (8)]. Intuitively, $r$ takes values between 0 and 1 and gives a sense of the overall level of phase coherence. $r≈1$ corresponds states that are more synchronized, i.e., with a probability distribution that is well approximated by a delta function. A fully desynchronized splay state, having a uniform distribution, returns $r=0$⁠. While the Kuramoto order parameter is not a perfect measure of synchrony (for example, rotating block solutions can also have an order parameter of 0), this metric is sufficient for the purposes of this work to quantify the level of synchronization of a distribution of oscillators. The order parameter associated with the simulations shown in Fig. 1 is shown in panel (c). Indeed, before coupling is applied $r(t)$ starts near zero and grows to $r(t)≈0.9$ after the coupling is turned on indicating a high degree of phase cohesiveness.

Due to the complexity and high-dimensionality of Eq. (1), phase-amplitude-based model order reduction techniques will be used for further analysis. To this end, consider a general dynamical system of the form

where $x∈RN$ is the system’s state, $F$ sets the dynamics, $p0∈RM$ is a collection of nominal parameters, and $u(t)∈R$ is an external input. Suppose that Eq. (5) admits a stable, $T$-periodic limit cycle $xγ$ when $u(t)=0$⁠.

The timing of oscillations occurring in Eq. (5) can be analyzed in a reduced order setting using the well-established method of phase reduction.^11,38,74 When using this strategy, the phase $θ∈[0,2π)$ is defined on the limit cycle and scaled so that $dθ/dt=2π/T=ω$ when $u(t)=0$⁠. The phase can be extended to include the basin of attraction of the limit cycle using the notion of isochrons,^15,74 which are defined as follows: for a given $p0$ and taking $u(t)=0$⁠, for any initial condition $a(0)∈xγ$⁠, the isochron corresponding to $a(0)$ is defined to be the set of all $b(0)$ such that

where $b(t)$ gives the solution of (5) with initial condition $b(0)$ and $||⋅||$ is some vector norm. Intuitively, initial conditions on the same isochron have the same asymptotic behavior. Defining the phase $θ(x)$ using isochrons, one can arrive at a reduced order model by restricting attention to a close neighborhood of the periodic orbit,

where $z(θ)≡⟨∂θ∂x,∂F∂u⟩$ is known as the phase response curve, where each partial derivative is evaluated at $xγ(θ)$ and $u=0$⁠. Note that in the transformation to (7), Taylor expansion is used to obtain $x˙=F(x,p0,0)+∂F∂uu(t)+O(u2)$ as an intermediate step. Phase reduction is a tremendously powerful tool that can be used to analyze oscillations of the $N$-dimensional system (5) using a one-dimensional reduction. However, it is only valid in the weak perturbation limit, i.e., requiring the underlying assumption that $u(t)=O(ϵ)$ where $0<ϵ≪1$⁠.

Phase-only models of the form (7) typically fail when large magnitude inputs are required. Such inputs drive the state of the system far from the vicinity of the periodic orbit to states for which the local linearizations of the gradient of the phase are no longer accurate. To combat this issue, phase-amplitude based coordinate systems can be employed to capture the dynamics in directions transverse to the periodic orbit. A variety of such coordinate systems have been proposed and investigated in recent years.^{8,14,23,48,55,61,65} In this work, the notion of isostable coordinates will be used, which represent level sets of the slowest decaying eigenfunctions of the Koopman operator.³⁰ As described in Ref. 69 (cf., Ref. 73) one can define the slowest decaying isostable coordinates in the basin of attraction of the limit cycle by exploiting Floquet theory. By first defining $Δx=x−xγ(θ)$⁠, one can write $Δx˙=J(θ)Δx+O(||Δx||2)$⁠, where $J(θ)$ denotes the Jacobian of the vector field evaluated at $xγ(θ)$⁠. Letting $Υ$ be the fundamental matrix associated with this linear equation [i.e., that solves $Δx(T)=ΥΔx(0)$], define the left eigenvectors, right eigenvectors, and eigenvalues of $Υ$ as $wj$⁠, $vj$⁠, and $λj$⁠, respectively. Also, define $κj=log⁡(λj)/T$ as the Floquet exponent associated with $λj$⁠. These Floquet exponents will be sorted so that $κN=0$ (corresponding to the $λN=1$ eigenvalue of the fundamental matrix) and that $|Re(κj)|≤|Re(κj+1)|$ for all others. Provided $λ1$ is not a defective eigenvalue,⁶⁶ the $ψ1$ isostable coordinate can be defined in the basin of attraction of the periodic orbit as follows:

where $χ(t,x)$ represents the flow of (5) when $u(t)=0$⁠, $x0Γ$ is the intersection of the periodic orbit and the $θ=0$ isochron, and $tΓ$ gives the time of the $kth$ transversal of the $θ=0$ isochron. Intuitively, isostable coordinates give a sense of the temporal distance from the periodic orbit, with larger isostable coordinates encoding for states that will take longer to relax to the periodic orbit. Isostable coordinates associated with more rapidly decaying Floquet exponents can be defined implicitly as level sets of Koopman eigenmodes. The interested reader is directed to Refs. 35, 70, and 73 for a more detailed discussion of the isostable coordinate paradigm and related reduced order modeling strategies.

This work will employ the adaptive phase-amplitude reduction approach to consider the aggregate dynamical behavior coupled neural oscillators using an analytically tractable reduced order coordinate basis. Following the adaptive reduction strategy illustrated in Refs. 67 and 68, Eq. (5) can be considered in terms of the behavior of a shadow system

where

One can verify that Eq. (9) is mathematically equivalent to Eq. (5). Here, $F(x,p,0)$ is taken to be the nominal dynamical behavior and $Ue$ is an effective additive input. It will be assumed that for some nominal range of system parameters $p$⁠, periodic solutions $xpγ$ exists and are continuously differentiable with respect to both $θ$ and $p$⁠. In reference to each periodic orbit, a generalized phase, $θ(x,p)$ and a set of slowly decaying generalized isostable coordinates $ψ1(x,p),…,ψβ(x,p)$ can be considered. Assuming that $θ(x,p)$ and each $ψj(x,p)$ are continuously differentiable in a neighborhood of the limit cycle for all $x$ and $p$⁠, one can transform to phase and isostable coordinates via the chain rule

for $j=1,…,β$⁠. In the above transformation, the nominal parameter set is allowed to change as a function of time. A detailed description of strategies for the computation of the partial derivatives from (11) are presented in Ref. 68. Evaluation of these terms yields the adaptive reduction

Above, $ω(p)$ is the natural frequency of the unperturbed periodic orbit $xpγ$⁠, $Z(θ,p)$ is the gradient of the phase coordinate evaluated on $xpγ$⁠, $D(θ,p)=∂θ∂p$ and captures how changes to the adaptive parameter set influence the generalized phase, $κj(p)$ is the Floquet exponent associated with the $jth$ isostable coordinate, $Ij(θ,p)$ is the gradient of the isostable coordinate $ψj$ with respect to the state evaluated on $xpγ$⁠, $Qj(θ,p)=∂ψj∂p$ and captures how changes to the adaptive parameter set influence the isostable coordinates, and $Gp$ is a parameter update rule that will be discussed in further detail momentarily. In the transformation from (9) to (12), the most rapidly decaying isostable coordinates (i.e., with large magnitude, negative Floquet exponents) can be assumed to be zero and neglected from the model dynamics, thereby providing a reduction in dimensionality.

As mentioned earlier, standard phase-based reduction strategies implicitly require that the state remains close to the underlying periodic orbit. The critical advantage of the adaptive reduction is its ability to continuously update the adaptive parameter set $p$ in order keep the isostable (amplitude) coordinates small, thereby keeping the state close to the periodic orbit $xpγ$⁠. To accomplish this, $Gp$ must be chosen appropriately. Some general design heuristics for $Gp$ are proposed in Ref. 68. For instance, when only a single non-truncated isostable coordinate is considered and $p∈R1$⁠, it often works well in practice to choose

taking $ν>0$⁠. For higher dimensional reductions, $Gp$ must be designed on an ad hoc basis. Provided an appropriate $Gp$ can be found so that each $ψj$ remains an $O(ϵ)$ term at all times, one can write the state as

Substituting Eq. (14) into (12) yields a closed reduced order equation that no longer requires direct knowledge of the state $x$⁠.

Here, using the adaptive phase-amplitude reduction framework in conjunction with formal averaging techniques, we illustrate that the general model of form (1) with an externally applied current of the form (2) can be accurately represented according to a low-dimensional set of autonomous ordinary differential equations. Ultimately, this reduced order model will be used to investigate the influence of $fζ$ and $fμ$ from (2) on the efficacy of the control law.

In the absence of coupling, input, and noise, we assume that each neuron from (1) has a strongly attracting periodic orbit with unperturbed period $Tn$ so that the each neuron can be represented as a limit cycle oscillator. Using the phase reduction from (7), the equations describing (1) become

for $i=1,…,N$⁠. Here, $θjn$ is the phase of the $jth$ neuron. Because the periodic orbit is strongly attracting, functions $V(θin)$ and $s(θin)$ are used to represent the transmembrane voltage and synaptic variable associated with the $ith$ neuron as a function of the phase.²² Additionally, $ωn=2π/Tn$⁠. As part of the transformation to phase coordinates, it is assumed that $D$ is small enough so that terms proportional to the square of the noise intensity can be truncated²⁶ (cf., Ref. 2). In the formulation above, $g0$ and $Istim$ are assumed to be small enough so that $si(t−τ)$ is well approximated by $s(θi−2πτ/Tn)$ for all $i$⁠. Note that in this work, we use the deterministic definition of phase from (8), which is defined for each neural oscillator with respect to the stable limit cycle in the absence of noise and other inputs. We emphasize that alternative formulations of phase are possible (see, for instance, Refs. 54 and 58) that directly incorporate the effects of noise.

Provided that the population of neurons is large enough, Eq. (15) can be represented according to a probability density $ρ(θn,t)$ that evolves according to the Fokker–Planck equation using an Ito formulation,¹³ 

where

and $ρθ$ and $zθ$ denote the partial derivatives taken with respect to $θn$⁠. Because $θn∈S1$⁠, the boundary conditions of this partial differential equation (PDE) are periodic. Equation (17) is analogous to Eq. (3) and is valid when considering probability distributions instead of a finite population of neurons. For all simulations of (16), we will take $τ=1.68$⁠. This value is chosen so that the coupling yields phase cohesive solutions in the absence of input. Toward a formulation using the adaptive reduction strategy described in Sec. II C, we distinguish between the internal dynamics, and the external input, separating terms that depend on $Istim(t)$⁠. The PDE model (16) will be further analyzed using a finite difference approximation with 100 total grid points for the $θ$ axis. Letting $xθn∈R100$ represent the discretized solution of the phase distribution, after discretization the dynamics of Eq. (16) can be approximated according to

where $zd∈R100$ is a discretized version of $z(θn)$⁠, $Fθn$ represents the resulting dynamics that are solely a function of the state and coupling strength, $A$ is a matrix which implements the finite difference approximation of the derivative of $ρ(θn,t)$ with respect to $θn$⁠, and $diag(⋅)$ takes a vector as its argument and returns a diagonal matrix of compatible size.

Equation (18) is in the same general form as Eq. (5) [with $Istim(t)$ taking the place of $u(t)$], allowing for the implementation of the adaptive phase-amplitude reduction strategy described in Sec. II C. Here, we will take the synaptic coupling strength to be the adaptive parameter, allowing Eq. (18) to be analyzed in terms of the shadow system

where the coupling strength $g∈[0.225,0.455]$ is taken to be the adaptive parameter. For each value of $g$ considered, the model has a stable periodic orbit $xgγ$ with period ranging from 9.47 to 13.57 ms. Here, the variable $θ$ will be used to denote the population phase that describes the timing of the population-level oscillations. Note that $θ$⁠, the phase of the population oscillation, is distinct from $θn$⁠, the phase of an individual oscillator with respect to its unperturbed limit cycle. Each periodic orbit, $xgγ(θ)$⁠, will be aligned so that $θ=0$ corresponds to the moment that $s¯$ reaches is maximum value on the limit cycle. Implementing the adaptive reduction strategy discussed in Sec. II C and considering the equations governing the applied stimulus from (2), the dynamics of Eq. (19) can be represented according to

with $Istim=μfs(t+ζTs/2π)$ and $ν=0.1$⁠. Above, the structure of $fζ$ and $fμ$ has been rewritten to depend on the reduced order model parameters $θ$⁠, $g$⁠, and $μ$ rather than on time itself. Additionally, it is assumed that all but one of the isostable coordinates decay rapidly and can be ignored. Per Eq. (14), the associated probability distribution can be approximated in terms of the reduced order model parameters according to

Note that the parameter update rule for $g˙$ from Eq. (20) matches the general strategy from (13) that serves to drive $ψ1$ to smaller magnitude values. Figure 2 shows various characteristics of the reduced order model from (20). Panel (a) shows $ρ(θn,g)$ plotted at $θ=0$ for various values of $g$ on the appropriate limit cycles. As $g$ decreases, the coherence between the oscillators decreases, resulting in a less synchronous population oscillation. Panel (b) shows the magnitudes of the principal Floquet multiplier (blue line) as well as secondary and tertiary complex-conjugate Floquet multipliers. The relatively fast decay of the complex-conjugate eigenvalues allows model (19) to be represented using a single isostable coordinate. Taking $g0=0.35mS/cm2$⁠, in response to a the periodic input $Istim(t)=1.2sin⁡(2πt/5)$⁠, panel (c) shows the evolution of the adaptive parameter over time. Panel (d) compares two sample probability distributions associated with the reduced order Eq. (20) and full order Eq. (18) with solid and dashed lines, respectively, at two separate times. Times associated with larger (respectively, smaller) values of $g$ are associated with more (respectively less) phase cohesive oscillations. Panel (e) shows that the isostable coordinate $ψ1$ does indeed remain small over time as required by the adaptive reduction framework, and panel (f) shows the magnitude of the Kuramoto order parameter,

computed both from the output of the reduced order model (black line) and the full model simulations (dashed line). Compared to Eq. (4), Eq. (22) is valid in the limit that $N$ tends toward infinity. Overall, we observe good agreement between the full order and reduced order simulations.

In Sec. III, we illustrated that evolution of the probability distribution of phases within a large population of coupled neural oscillators can be represented in terms of a low order set of equations of form (20) by using the adaptive phase-amplitude reduction paradigm. By recognizing that $Ue$ from Eq. (20) is periodic in $t$⁠, we illustrate that further reduction in dimensionality can be realized using formal averaging techniques.^16,51

To begin, recalling from Eq. (2) that $Ts$ is the nominal period of the input $fs$⁠, we will work in a rotating reference frame defining

where $Ω≡2π/Ts$⁠. Note that we have incorporated $ζ$ in the definition of our rotating reference frame, which explicitly accounts for the slowly varying phase offset of the applied input.

Working in this rotating reference frame, Eq. (20) becomes

Recall from Eq. (19) that $Ue(xθn,g,g0,Istim(t))$ is periodic in time with period $Ts$⁠. We will assume that $Ue$⁠, $ψ1$⁠, and $ω(g)−Ω$ are $O(ϵ)$ terms so that Eq. (24) can be written in the general form $x˙=ϵFR(x,t)$ where $x=[ϕψ1gζμ]T$⁠, and $FR$ gives the dynamics and is periodic in $t$⁠. For such a system, the method of averaging^16,51 can be used to approximate the behavior of (24) according to

where

Using the averaging framework, Eq. (25) provides a close approximation to the unaveraged system of equations from (24). Note that the terms $A1$ through $A7$ do not explicitly depend on $ζ$⁠, which shifts the periodic terms of each integrand but does not influence the result when integrating over the entire period. Also, $A1$ through $A7$ are assumed to be $O(1)$ terms; even though $A1$ and $A4$ are multiplied by $1/ϵ$ , the term $Ue$ in the integrand is an $O(ϵ)$ term.

Further simplification of Eq. (25) is possible provided $ν$ is chosen large enough so that $ψ1$ rapidly converges to its steady state value. To this end, define

as the steady state value of $ψ1$ that results when $ϕ,g,ζ$⁠, and $μ$ are held fixed. Recognizing that the right hand sides of Eq. (25) are comprised of $O(ϵ)$ terms, $ϕ˙,g˙,μ˙$⁠, and $ζ˙$ are also $O(ϵ)$ terms. Taking the time derivative of $ψss$ one can show that

Next, defining $Δψ1≡ψ1−ψss$ and taking the time derivative, one finds

Because all periodic orbits used in the adaptive reduction are stable, $κ(g)<0$⁠. Likewise, $A5(ϕ,g)>0$ since the integrand from (26) is non-negative and $ν>0$ by definition. Considering that $κ(g)−νA5(ϕ,g)=O(1)$ provided that $Δψ1=O(ϵ2)$ at $t=0$⁠, $Δψ1=O(ϵ2)$ uniformly bounded in time and, consequently, $ψ1=ψss$ can be used as an approximation that is valid up to $O(ϵ2)$⁠. This allows for the elimination of the $ψ1$ dynamics from Eq. (25) to yield

Finally, provided $|νA5(ϕ,g)|≫|κ(g)|$⁠, one can write

Substituting (31) into Eq. (30) finally yields

where

Equation (32) will be used to investigate the response of the neural oscillator population to the near-periodic input described by Eq. (2). In the averaged setting, the phase offset and the magnitude of the input can be modulated slowly through the functions $A3$ and $A7$⁠, respectively. Strategies for designing these inputs will be considered momentarily. Note that in situations where $ζ˙=μ˙=0$⁠, i.e., when purely periodic inputs are applied, only two state variables are required to represent the perturbed dynamics. We emphasize that in the reduction from Eq. (25) to Eq. (32), the term $ν$ disappears provided that $|νA5(ϕ,g)|$ is sufficiently larger than $|κ(g)|$ for all choices of $ϕ$ and $g$⁠.

In Fig. 3, the temporal evolution of $ϕ$ is compared for simulations of the averaged model equations from (32) and the unaveraged model equations from (20) when using two different periodic sinusoidal inputs. In this example, we take $ζ=ζ˙=μ˙=0$ so that the input is fully periodic. Note here that for the input $Istim(t)=1.2sin⁡(2πt/5)$⁠, when computing the averaged equation $Ω$ is taken to be $2π/10$ in order to keep the magnitude of $ω(g)−Ω$ small. After computing the averaged functions $A1$–$A7$ using (26), $ϕ˙$ (respectively, $g˙$⁠) computed according to Eq. (32) is shown in panels (a) and (d) [resp., (b) and (e)] for the two different stimuli. Overall, the averaged equations from (32) accurately represent the unaveraged dynamics from Eq. (20) with results of individual simulations shown in panels (c) and (f). Nullclines are plotted with red and green dashed lines to highlight the salient dynamical behavior.

Results in Figs. 2 and 3 illustrate that both the reduced order equations obtained using the adaptive reduction strategy (20) and the averaged equations from (32) provide a good approximation to the dynamics of (16). In Sec. IV, we will use these reduced order equations to develop control strategies for manipulating the population-level oscillations with the goal of disrupting synchronization.

The results from Fig. 3 correspond to a situation where $fζ$ and $fμ$ from Eq. (2) are zero, that is, where the applied input is purely periodic. As we will illustrate here, by slowly changing the magnitude of the input and the relative offset of the rotating reference frame, $μ$ and $ζ$⁠, respectively, the cohesion between oscillators from (1) can be effectively manipulated.

Our goal is to identify functions $fζ$ and $fμ$ from Eq. (20) that can shift the state of the adaptive parameter in the reduced order equations from (20) from $g=g0$ to some $g=gtarg<g0$⁠. By shifting $g$ to smaller values, the resulting phase distribution will be less synchronized. We emphasize that in the implementation of the control algorithm, the actual strength of the synaptic coupling does not change. Rather, the goal is to apply input using $Istim(t)$ that drives the system to less synchronized states that are near the periodic orbit associated with the parameter $gtarg$⁠.

Intuitively, toward the identification of an appropriate control strategy, first consider the $g˙=0$ nullcline from panel (c) of Fig. 3 obtained for the averaged equations (32). Referring to this nullcline by $g=fn(ϕ)$⁠, for any $gtarg$ in the range of $fn(ϕ)$⁠, the control objective could be satisfied provided that $ϕ$ could be set to some value $ϕtarg$ for which $gtarg=fn(ϕtarg)$⁠. For instance, considering the $g˙=0$ nullcline from panel (c) of Fig. 3, if a control strategy can be found that sets $ϕ$ to a constant value of 5 rad, $g$ would tend toward 0.28. Likewise, considering the nullclines from panel (f) of Fig. 3, by actively controlling $ϕ$ to values between 4.7 and 5.6 rad, $g$ can be driven to any desired value less then $g0$⁠. In this case, however, the magnitude of the input, $μ$ would need to be decreased as $g$ is driven near its lower limit. With this in mind, we will consider a control strategy that seeks to drive $g$ and $ϕ$ to $gtarg$ and $ϕtarg$⁠, respectively. This control strategy will be implemented by mandating the following update rules for $ζ$ and $μ$ from Eq. (20):

where $K1$–$K4$ are appropriately chosen, positive constants. Here, in the averaged reference frame, i.e., using (32), this choice of $fζ$ and $fμ$ yields the relationships

where $ai$ is an extra variable used to represent the necessary terms from the integral control. Intuitively in the above equations, the term $K1(ϕ−ϕtarg)$ acts as a proportional controller that serves to update the offset of the rotating reference frame, $ζ$⁠, to manipulate $ϕ$⁠. The term $K2(g−gtarg)+aI$ serves as a proportional–integral controller that increases or decreases the magnitude of the applied input depending on whether $g$ is above or below its target value. The additional term, $−K4μ$⁠, is added to influence stability of the controller. Of course, the efficacy of the controller (34) depends on an adequate choice for $ϕtarg$ and $gtarg$⁠. These considerations will be discussed in the results to follow.

Implementation of the control strategy proposed in Eq. (35) requires knowledge of the functions $H1(ϕ,g,μ)$ and $H2(ϕ,g,μ)$⁠. If an underlying model of form (20) is known, these equations can be identified by first directly evaluating the integrals from Eq. (26) and then computing the necessary terms from Eq. (33). In the absence of an underlying dynamical model, the required functions $H1(ϕ,g,μ)$ and $H2(ϕ,g,μ)$ can still be estimated in real-time from observable data using a strategy described below.

Consider a general system of form (5) with a set of reduced order, averaged equations that take the form (32). As was the case for the reduction of the noisy oscillator population from (15), assume that $θ=0$ can be clearly identified from measured data [for instance, $θ=0$ for a given periodic orbit from Eq. (16) corresponds to the moment that $s¯$ reaches a local maximum]. Consider two successive traversals of the $θ=0$ level set, occurring at $t1$ and $t2$⁠. The associated coordinates in the rotating reference frame defined by Eq. (23) are

Toward identifying a strategy to estimate the unknown functions $H1$ and $H2$ from observable data, one can rewrite Eq. (32) as

Above, the arguments of $H1$ and $H2$ are evaluated at the fixed time $t2$⁠. Because $ϕ˙$⁠, $g˙$⁠, $μ˙$⁠, $H1$⁠, and $H2$ are $O(ϵ)$ terms, fixing time in the evaluation of $H1$ and $H2$ only contributes an additional $O(ϵ2)$ term when $t∈[t1,t2]$⁠. Neglecting the $O(ϵ2)$ and $O(ϵκ(g)/ν)$ terms, through direct integration of Eq. (37), one can write

Equation (38) can be used to obtain successive, real-time estimates of $H1$ and $H2$⁠, for instance, when implementing the control strategy suggested by Eq. (35), $ζ(t)$ is known exactly at all times because it is part of the controller, $ϕ(t1)$ and $ϕ(t2)$ can be computed according to (23) at successive crossings of the $θ=0$ level set, and $g(t2)$ and $g(t1)$ can be inferred from output data. Using this strategy, a reduced order model does not need to be computed before implementing the control strategy described by Eq. (35). Rather, the necessary terms of the reduced order equations can be computed in real-time as they are needed.

As a matter of practical implementation, $H1$ and $H2$ can be updated once per cycle and are assumed to be constant between updates. In order to mitigate small errors that would otherwise accumulate over time, the reduced order variables $ϕ$ and $g$ can also be updated once per cycle to reflect the measurements from data. The approximations from Eq. (38) require the right hand side of Eq. (32) to be small. As such, the accuracy of these approximations will degrade as larger magnitude inputs are used.

The model independent control strategy suggested in Sec. IV B relies on the estimation of the terms from Eq. (38) from observables and subsequent implementation of the controller (35). A step-by-step procedure for accomplishing this task is provided below.

• (Step 1)

  Identify an observable $y(t)∈R1$ to be used for the implementation of the control strategy. This observable should have clearly identifiable oscillations. Define $θ=0$ to correspond to a distinct feature of the oscillations that occurs once per cycle, for instance, the crossing of a Poincaré section or the attainment of a local maximum.

• (Step 2)

  Defining $t1$ and $t2$ to be the times of subsequent crossings of the $θ=0$ isochron, define a map from the observation of $y(t)$ on the interval $[t1,t2]$ to the variable $g(t2)$ that encodes for the oscillation amplitude. For instance, one could take $g(t2)=max(y(t))−min(y(t))$ on the interval $t=[t1,t2]$⁠.

• (Step 3)

  Choose a periodic stimulation $Istim(t)$ to be applied to the system. The period of $Istim(t)$ must be close to the period of oscillation of the observable chosen in step 1.

• (Step 4)

  Choose a set of constants $K1,K2,K3,K4$⁠, and $ϕtarg$ for the controller from Eq. (35). In a model independent setting, this will generally be a trial and error process. One can start with $ϕtarg$ and $K1$ to identify choices of $ϕtarg$ that effectively decreases the magnitude of oscillations. The terms $K2$⁠, $K3$⁠, and $K4$ can subsequently be tuned to appropriately modulate the stimulus magnitude.

• (Step 5)

  Choose a profile for $gtarg(t)$⁠, which sets the target for the oscillation amplitude.

• (Step 6)

  Implement the controller from Eq. (35). On a cycle-by-cycle basis, estimates of $H1$ and $H2$ can be updated according to Eq. (38), estimates of $ϕ(t1)$ and $ϕ(t2)$ can be computed according to (23) at successive crossings of the $θ=0$ level set, and $g(t2)$ and $g(t1)$ can be inferred from the mapping defined in Step 2.

• (Step 7)

  Adjust constants $K1,K2,K3,K4$⁠, and $ϕtarg$ as necessary. Performance is generally less sensitive to the choice of these constants when $gtarg$ is close to the nominal value of $g$ that results in the absence of $Istim$⁠. As $gtarg$ is decreased (i.e., as the oscillation amplitude decreases) parameter adjustments may be necessary to yield a viable controller.