Deep brain stimulation (DBS) is a commonly used treatment for medication resistant Parkinson’s disease and is an emerging treatment for other neurological disorders. More recently, phase-specific adaptive DBS (aDBS), whereby the application of stimulation is locked to a particular phase of tremor, has been proposed as a strategy to improve therapeutic efficacy and decrease side effects. In this work, in the context of these phase-specific aDBS strategies, we investigate the dynamical behavior of large populations of coupled neurons in response to near-periodic stimulation, namely, stimulation that is periodic except for a slowly changing amplitude and phase offset that can be used to coordinate the timing of applied input with a specified phase of model oscillations. Using an adaptive phase-amplitude reduction strategy, we illustrate that for a large population of oscillatory neurons, the temporal evolution of the associated phase distribution in response to near-periodic forcing can be captured using a reduced order model with four state variables. Subsequently, we devise and validate a closed-loop control strategy to disrupt synchronization caused by coupling. Additionally, we identify strategies for implementing the proposed control strategy in situations where underlying model equations are unavailable by estimating the necessary terms of the reduced order equations in real-time from observables.

High frequency deep brain stimulation (DBS) is a well-established treatment for Parkinson’s disease. In most clinical settings, DBS is implemented by injecting 130–180 Hz pulse trains of electrical current into an appropriate brain region. Recent evidence suggests that phase-specific adaptive deep brain stimulation (aDBS), whereby the application of stimulation is locked to a particular phase of tremor, may improve therapeutic efficacy and reduce side effects. 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. Here, we apply phase-amplitude-based coordinate reduction strategies to identify a low-dimensional model and subsequently devise and validate a closed-loop control strategy to manipulate the aggregate behavior of the neural population. This work can be viewed as potential theoretical underpinning for the phase-specific adaptive DBS strategies that have gained interest in recent years.

## I. INTRODUCTION

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,^{25} depression,^{31} tinnitus,^{56} and Tourette syndrome.^{53} 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,^{36} identifying energy optimal stimuli that achieve positive Lyapunov exponents and produce chaotic desynchronization,^{64,71} controlling neurons to a prespecified phase distribution,^{34} 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^{63} 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^{46} or the severity of tremor^{27} 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.

## II. COMPUTATIONAL MODEL, CONTROL OBJECTIVE, AND BACKGROUND ON MODEL ORDER REDUCTION TECHNIQUES

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\u2208RNa$ 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,\sigma T$, and $c2$ determine the dynamics of the synaptic variable, $Esyn=\u2212100$ mV is the reversal potential of the neurotransmitter so that the coupling is inhibitory, $\tau $ 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 $\eta i(t)$ is a zero-mean white noise process associated with neuron $i$ for which $D=1\mu A2/cm4$ sets the intensity. Here, the collection $\eta 1(t),\u2026,\eta 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,^{17} for concreteness, we will use a model of thalamic neurons from^{50} 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=\u221220mV,\sigma 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.

### A. Description of the primary control objective

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 $\mu $ and phase offset $\zeta $. Both $\zeta $ and $\mu $ are allowed to change slowly according to the relations $f\zeta $ and $f\mu $, 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\zeta $ and $f\mu $ 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 $\tau =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 $\theta 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\u22481$ 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)\u22480.9$ after the coupling is turned on indicating a high degree of phase cohesiveness.

### B. Phase-amplitude reduction methods

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\u2208RN$ is the system’s state, $F$ sets the dynamics, $p0\u2208RM$ is a collection of nominal parameters, and $u(t)\u2208R$ is an external input. Suppose that Eq. (5) admits a stable, $T$-periodic limit cycle $x\gamma $ 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 $\theta \u2208[0,2\pi )$ is defined on the limit cycle and scaled so that $d\theta /dt=2\pi /T=\omega $ 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)\u2208x\gamma $, 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 $||\u22c5||$ is some vector norm. Intuitively, initial conditions on the same isochron have the same asymptotic behavior. Defining the phase $\theta (x)$ using isochrons, one can arrive at a reduced order model by restricting attention to a close neighborhood of the periodic orbit,

where $z(\theta )\u2261\u27e8\u2202\theta \u2202x,\u2202F\u2202u\u27e9$ is known as the phase response curve, where each partial derivative is evaluated at $x\gamma (\theta )$ and $u=0$. Note that in the transformation to (7), Taylor expansion is used to obtain $x\u02d9=F(x,p0,0)+\u2202F\u2202uu(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(\u03f5)$ where $0<\u03f5\u226a1$.

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.^{30} 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 $\Delta x=x\u2212x\gamma (\theta )$, one can write $\Delta x\u02d9=J(\theta )\Delta x+O(||\Delta x||2)$, where $J(\theta )$ denotes the Jacobian of the vector field evaluated at $x\gamma (\theta )$. Letting $\Upsilon $ be the fundamental matrix associated with this linear equation [i.e., that solves $\Delta x(T)=\Upsilon \Delta x(0)$], define the left eigenvectors, right eigenvectors, and eigenvalues of $\Upsilon $ as $wj$, $vj$, and $\lambda j$, respectively. Also, define $\kappa j=log\u2061(\lambda j)/T$ as the Floquet exponent associated with $\lambda j$. These Floquet exponents will be sorted so that $\kappa N=0$ (corresponding to the $\lambda N=1$ eigenvalue of the fundamental matrix) and that $|Re(\kappa j)|\u2264|Re(\kappa j+1)|$ for all others. Provided $\lambda 1$ is not a defective eigenvalue,^{66} the $\psi 1$ isostable coordinate can be defined in the basin of attraction of the periodic orbit as follows:

where $\chi (t,x)$ represents the flow of (5) when $u(t)=0$, $x0\Gamma $ is the intersection of the periodic orbit and the $\theta =0$ isochron, and $t\Gamma $ gives the time of the $kth$ transversal of the $\theta =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.

### C. Adaptive phase-amplitude reduction

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\gamma $ exists and are continuously differentiable with respect to both $\theta $ and $p$. In reference to each periodic orbit, a generalized phase, $\theta (x,p)$ and a set of slowly decaying generalized isostable coordinates $\psi 1(x,p),\u2026,\psi \beta (x,p)$ can be considered. Assuming that $\theta (x,p)$ and each $\psi 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,\u2026,\beta $. 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, $\omega (p)$ is the natural frequency of the unperturbed periodic orbit $xp\gamma $, $Z(\theta ,p)$ is the gradient of the phase coordinate evaluated on $xp\gamma $, $D(\theta ,p)=\u2202\theta \u2202p$ and captures how changes to the adaptive parameter set influence the generalized phase, $\kappa j(p)$ is the Floquet exponent associated with the $jth$ isostable coordinate, $Ij(\theta ,p)$ is the gradient of the isostable coordinate $\psi j$ with respect to the state evaluated on $xp\gamma $, $Qj(\theta ,p)=\u2202\psi j\u2202p$ 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\gamma $. 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\u2208R1$, it often works well in practice to choose

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

## III. ANALYSIS OF NEAR-PERIODIC FORCING IN AN OSCILLATORY POPULATION OF SYNAPTICALLY COUPLED NEURONS

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\zeta $ and $f\mu $ from (2) on the efficacy of the control law.

### A. Representation of the noisy population of oscillators using the Fokker–Planck equation

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,\u2026,N$. Here, $\theta jn$ is the phase of the $jth$ neuron. Because the periodic orbit is strongly attracting, functions $V(\theta in)$ and $s(\theta in)$ are used to represent the transmembrane voltage and synaptic variable associated with the $ith$ neuron as a function of the phase.^{22} Additionally, $\omega n=2\pi /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^{26} (cf., Ref. 2). In the formulation above, $g0$ and $Istim$ are assumed to be small enough so that $si(t\u2212\tau )$ is well approximated by $s(\theta i\u22122\pi \tau /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 $\rho (\theta n,t)$ that evolves according to the Fokker–Planck equation using an Ito formulation,^{13}

where

and $\rho \theta $ and $z\theta $ denote the partial derivatives taken with respect to $\theta n$. Because $\theta n\u2208S1$, 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 $\tau =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 $\theta $ axis. Letting $x\theta n\u2208R100$ represent the discretized solution of the phase distribution, after discretization the dynamics of Eq. (16) can be approximated according to

where $zd\u2208R100$ is a discretized version of $z(\theta n)$, $F\theta 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 $\rho (\theta n,t)$ with respect to $\theta n$, and $diag(\u22c5)$ takes a vector as its argument and returns a diagonal matrix of compatible size.

### B. Representation of the temporal evolution of the probability distribution using the adaptive reduction paradigm

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\u2208[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\gamma $ with period ranging from 9.47 to 13.57 ms. Here, the variable $\theta $ will be used to denote the population phase that describes the timing of the population-level oscillations. Note that $\theta $, the phase of the population oscillation, is distinct from $\theta n$, the phase of an individual oscillator with respect to its unperturbed limit cycle. Each periodic orbit, $xg\gamma (\theta )$, will be aligned so that $\theta =0$ corresponds to the moment that $s\xaf$ 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=\mu fs(t+\zeta Ts/2\pi )$ and $\nu =0.1$. Above, the structure of $f\zeta $ and $f\mu $ has been rewritten to depend on the reduced order model parameters $\theta $, $g$, and $\mu $ 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\u02d9$ from Eq. (20) matches the general strategy from (13) that serves to drive $\psi 1$ to smaller magnitude values. Figure 2 shows various characteristics of the reduced order model from (20). Panel (a) shows $\rho (\theta n,g)$ plotted at $\theta =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\u2061(2\pi 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 $\psi 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,

### C. Analysis of reduced order model dynamics using formal averaging techniques

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 $\Omega \u22612\pi /Ts$. Note that we have incorporated $\zeta $ 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\theta n,g,g0,Istim(t))$ is periodic in time with period $Ts$. We will assume that $Ue$, $\psi 1$, and $\omega (g)\u2212\Omega $ are $O(\u03f5)$ terms so that Eq. (24) can be written in the general form $x\u02d9=\u03f5FR(x,t)$ where $x=[\varphi \psi 1g\zeta \mu ]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 $\zeta $, 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/\u03f5$ , the term $Ue$ in the integrand is an $O(\u03f5)$ term.

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

as the steady state value of $\psi 1$ that results when $\varphi ,g,\zeta $, and $\mu $ are held fixed. Recognizing that the right hand sides of Eq. (25) are comprised of $O(\u03f5)$ terms, $\varphi \u02d9,g\u02d9,\mu \u02d9$, and $\zeta \u02d9$ are also $O(\u03f5)$ terms. Taking the time derivative of $\psi ss$ one can show that

Next, defining $\Delta \psi 1\u2261\psi 1\u2212\psi ss$ and taking the time derivative, one finds

Because all periodic orbits used in the adaptive reduction are stable, $\kappa (g)<0$. Likewise, $A5(\varphi ,g)>0$ since the integrand from (26) is non-negative and $\nu >0$ by definition. Considering that $\kappa (g)\u2212\nu A5(\varphi ,g)=O(1)$ provided that $\Delta \psi 1=O(\u03f52)$ at $t=0$, $\Delta \psi 1=O(\u03f52)$ uniformly bounded in time and, consequently, $\psi 1=\psi ss$ can be used as an approximation that is valid up to $O(\u03f52)$. This allows for the elimination of the $\psi 1$ dynamics from Eq. (25) to yield

Finally, provided $|\nu A5(\varphi ,g)|\u226b|\kappa (g)|$, one can write

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 $\zeta \u02d9=\mu \u02d9=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 $\nu $ disappears provided that $|\nu A5(\varphi ,g)|$ is sufficiently larger than $|\kappa (g)|$ for all choices of $\varphi $ and $g$.

In Fig. 3, the temporal evolution of $\varphi $ 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 $\zeta =\zeta \u02d9=\mu \u02d9=0$ so that the input is fully periodic. Note here that for the input $Istim(t)=1.2sin\u2061(2\pi t/5)$, when computing the averaged equation $\Omega $ is taken to be $2\pi /10$ in order to keep the magnitude of $\omega (g)\u2212\Omega $ small. After computing the averaged functions $A1$–$A7$ using (26), $\varphi \u02d9$ (respectively, $g\u02d9$) 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.

## IV. A CONTROL STRATEGY TO MANIPULATE THE PHASE COHESION OF THE OSCILLATOR POPULATION

The results from Fig. 3 correspond to a situation where $f\zeta $ and $f\mu $ 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, $\mu $ and $\zeta $, respectively, the cohesion between oscillators from (1) can be effectively manipulated.

### A. A controller for the input magnitude and phase offset

Our goal is to identify functions $f\zeta $ and $f\mu $ 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\u02d9=0$ nullcline from panel (c) of Fig. 3 obtained for the averaged equations (32). Referring to this nullcline by $g=fn(\varphi )$, for any $gtarg$ in the range of $fn(\varphi )$, the control objective could be satisfied provided that $\varphi $ could be set to some value $\varphi targ$ for which $gtarg=fn(\varphi targ)$. For instance, considering the $g\u02d9=0$ nullcline from panel (c) of Fig. 3, if a control strategy can be found that sets $\varphi $ 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 $\varphi $ 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, $\mu $ 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 $\varphi $ to $gtarg$ and $\varphi targ$, respectively. This control strategy will be implemented by mandating the following update rules for $\zeta $ and $\mu $ 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\zeta $ and $f\mu $ 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(\varphi \u2212\varphi targ)$ acts as a proportional controller that serves to update the offset of the rotating reference frame, $\zeta $, to manipulate $\varphi $. The term $K2(g\u2212gtarg)+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, $\u2212K4\mu $, is added to influence stability of the controller. Of course, the efficacy of the controller (34) depends on an adequate choice for $\varphi targ$ and $gtarg$. These considerations will be discussed in the results to follow.

### B. Real-time estimation of reduced order terms for model independent control

Implementation of the control strategy proposed in Eq. (35) requires knowledge of the functions $H1(\varphi ,g,\mu )$ and $H2(\varphi ,g,\mu )$. 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(\varphi ,g,\mu )$ and $H2(\varphi ,g,\mu )$ 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 $\theta =0$ can be clearly identified from measured data [for instance, $\theta =0$ for a given periodic orbit from Eq. (16) corresponds to the moment that $s\xaf$ reaches a local maximum]. Consider two successive traversals of the $\theta =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 $\varphi \u02d9$, $g\u02d9$, $\mu \u02d9$, $H1$, and $H2$ are $O(\u03f5)$ terms, fixing time in the evaluation of $H1$ and $H2$ only contributes an additional $O(\u03f52)$ term when $t\u2208[t1,t2]$. Neglecting the $O(\u03f52)$ and $O(\u03f5\kappa (g)/\nu )$ 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), $\zeta (t)$ is known exactly at all times because it is part of the controller, $\varphi (t1)$ and $\varphi (t2)$ can be computed according to (23) at successive crossings of the $\theta =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 $\varphi $ 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.

### C. List of steps to implement the model independent control strategy

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)\u2208R1$ to be used for the implementation of the control strategy. This observable should have clearly identifiable oscillations. Define $\theta =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 $\theta =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))\u2212min(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 $\varphi 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 $\varphi targ$ and $K1$ to identify choices of $\varphi 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 $\varphi (t1)$ and $\varphi (t2)$ can be computed according to (23) at successive crossings of the $\theta =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 $\varphi 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.

## V. RESULTS

Here, we present results for the implementation of the control strategy from Sec. IV applied to both the Fokker–Planck model for the temporal evolution of the probability distribution of a large population of neurons from Eq. (16) and for a full model of $N=1000$ coupled neurons from Eq. (1). As described in Sec. III B, the Fokker–Planck model can be represented in a reduced order form using the adaptive reduction from Eq. (20). Consequently, the terms $H1(\varphi ,g,\mu )$ and $H2(\varphi ,g,\mu )$ can be computed directly from the terms of the reduced order equations using the averaging framework. For the full model of $N=1000$ coupled neurons, direct computation of $H1$ and $H2$ is not possible and the real-time estimation strategy from Sec. IV B will be used.

### A. Control of probability distributions of the Fokker–Planck equation

Here, we consider results of the control strategy from Eq. (35) applied to modify the phase cohesion of the Fokker–Planck model from Eq. (16). Starting with Eq. (16), the terms of the adaptive reduction from (20) are first computed taking $g0=0.35mS/cm2$. Recall that for the allowable range of coupling strengths $g\u2208[0.225,0.455]$, the associated periodic orbits have periods $T\u2208[9.47,13.57]$ ms. For this application, we consider the nominal input to be a periodic series of pulses,

Here, $Ts=10$ ms and the input is a series of six pulses shown taking $\mu =1$ in panel (b) of Fig. 4. Using this input, the terms of the averaged, reduced order equation of the form (32) are determined by directly evaluating the integrals from Eq. (26) and using them to compute the relationships from Eq. (33). Panel (a) of Fig. 4 shows level sets of $H2(\varphi ,g,\mu )=0$ shown at different values of $\mu $. Relative to the baseline at $g0=0.35mS/cm2$, $g\u02d9$ tends to decrease when $\varphi \u2208[4.2,6.2]$ and tends to increase for $\varphi \u2208[0.4,4]$. With this in mind, recalling that our goal is to decrease the phase cohesion by driving the adaptive parameter, $g$, to lower values, to implement the control strategy from Eq. (35), we choose $\varphi targ=5$ rad. Three different parameter sets $C1,C2,$ and $C3$ are chosen to implement the control strategy taking ${K1,K2,K3,K4}={4,35,0,0}$, ${4,35,1,0}$, and ${4,35,1,0.2}$, respectively. Parameter set $C1$ uses only two proportional controllers to determine the update rules for $\zeta $ and $\mu $. Parameter set $C2$ adds an integral controller to the update rule for $\mu $, and parameter set $C3$ adds an additional term that is proportional to $\mu $. The results to follow are qualitatively similar when modifying the non-zero values in these parameter sets.

We investigate the efficacy of these three parameter sets by considering fixed points of

i.e., with the same dynamics as Eq. (35) but neglecting the $\zeta $ variable. Here, $\zeta $ is ignored to account for the fact that $\varphi $ might not fully reach its target value. Fixed points of Eq. (40) are obtained with a Newton iteration for various choices of $gtarg$ and plotted against eigenvalues with the largest real component, $\lambda c$, of the linearization about the fixed point. For parameter sets $C1$ and $C2$, stability is lost for lower values of $gtarg$ as $Real(\lambda c)$ crosses from negative to positive values. Panel (d) illustrates the control strategy using parameter sets $C1,C2,$ and $C3$ for simulations of the reduced order equations from (35) taking $gtarg(t)=0.24+0.11exp\u2061(\u22120.01t)$. The parameter set $C3$ is able to follow this target accurately while $C1$ and $C2$ display unstable behavior for small values of $gtarg$. The emergence of this instability is predicted by the eigenvalue plots from panel (c).

Figure 4 considers results for the proposed control strategy applied directly to the reduced order equations from (35). Next, we consider the behavior of the proposed control strategy applied to the Fokker–Planck equation from (16). Once again, we use the pulsatile input specified by Eq. (39) with $Ts=10$ ms and take $\varphi targ=5$ rad. Like in the previous example, $H1$ and $H2$ from Eq. (35) are obtained by directly evaluating the terms of (33) associated with the averaged, reduced order model dynamics from Eq. (25). Recall that these reduced order model dynamics are approximations to the full order behavior of the Fokker–Planck equation. In order to mitigate the accumulation of small errors over time, at every crossing of $\theta =0$ [corresponding to the moment that $s\xaf$ from Eq. (17) achieves a local maximum] both $\varphi $ and $\mu $ are inferred and updated based on the measurements of $s\xaf(t)$. For each of these updates, $\varphi $ is inferred from Eq. (23) and $\mu $ is inferred by comparing the value of $s\xaf$ to the maximum values of $s\xaf$ attained for each periodic orbit $xg\gamma $. Results are shown in Fig. 5. Parameter sets $C1,C2$, and $C3$ are identical to those from Fig. 4 with the exception of taking $K2=20$ instead of 35. Parameters $\varphi $, $\zeta $, and $\mu $ are shown over the course of each simulation in panels (a)–(c) and the adaptive parameter $g$ is shown in panel (d). In this example, $gtarg(t)=0.35$ for $t<100$ and $gtarg(t)=0.232+0.118exp\u2061(\u22120.01(t\u2212100))$ for $t\u2265100$. Results are qualitatively similar to those from Fig. 4 with parameter sets $C1$ and $C2$ yielding controllers with unstable behavior.

The implementation of the proposed control strategy from the examples in Figs. 4 and 5 is predicated on the existence of a collection of periodic orbits that change continuously with respect to the adaptive parameter. In these examples, when the coupling strength $g$ is taken to be the adaptive parameter, periodic orbits do not exist when taking $g<0.225$. This represents a challenge when the overall goal is desynchronization, since the Kuramoto order parameter associated with the periodic orbit that emerges for $g\u22480.225$ is still relatively large. In order to address this limitation, in the example to follow, we will assume that a less phase cohesive set of periodic orbits, $xp\gamma $, exists for some adaptive parameter, $p\u2208R$. Once again, these orbits will be aligned so that $\theta =0$ corresponds to the moment that $s\xaf$ from Eq. (17) reaches a local maximum. We will assume the existence of some model observable $y(t)$ such that $y(t)=p(t)$. Ultimately, we will employ the real-time estimation strategy described in Sec. IV B to identify the reduced order terms used in the control strategy from Eq. (35). Note that in this example, the adaptive parameter set will be defined implicitly. Because we will be using the model observable to infer the value of the adaptive parameter, we do not need an explicit definition for the adaptive parameter.

To proceed, we will consider simulations of the Fokker–Planck equation with the observable,

Here, this observable is defined to be the difference between the maximum and minimum value of $s\xaf$ on the periodic orbit $xp\gamma (\theta )$. This definition allows for the inference from measurements of $sD(t)$ once per cycle. Note that $s\xaf$ still needs to be measured continuously in order to evaluate Eq. (41). Using this strategy, the Fokker–Planck equation from (16) is simulated and the controller from Eq. (35) is implemented by estimating $H1$ and $H2$ on a cycle-by-cycle basis. In this example, $sD$ takes the place of the adaptive parameter $g$ from (35). For the control parameters, we use $\varphi targ=5$ and take $sD,targ(t)=0.31$ for $t<100$ and $sD,targ(t)=0.01+0.30exp\u2061(\u22120.005(t\u2212100))$ for $t\u2265100$. These values of $sD,targ$ are chosen to be close to the nominal value that results in the absence of control for $t<100$ and to gradually decrease toward lower values (which correspond to less synchronous solutions) over time. Controller parameters $K1$, $K2$, $K3$, and $K4$ are taken to be 4, 30, 1, and 1, respectively. Qualitatively similar results can be obtained using other values for $K1$ through $K4$. Results are shown in Fig. 6. Panels (a)–(c) show $\varphi ,\zeta $, and $\mu $ over the course of the simulation. These traces are qualitatively similar to those obtained using the $C3$ parameter set from Fig. 5. The average value of the synaptic variable and the order parameter from Eqs. (17) and (22), respectively, are shown in panels (d) and (e). Both $s\xaf$ and $r$ oscillate rapidly relative to the timescale shown giving the appearance of a thick line. Over the course of the simulation, the order parameter is driven below the minimum values observed in the results from Fig. 5. Panel (f) shows the inferred value of $sD(t)$ during the simulation, which gradually decays toward smaller values, as mandated by the chosen value of $sD,targ(t)$. We emphasize that we are able to drive the probability distribution to less synchronous values without performing any *a priori* computations to obtain a reduced order model.

### B. Control of a large population of conductance-based neurons

Next, we consider the results of the proposed control strategy applied directly to the population of $N=1000$ noisy, conductance-based neurons characterized by Eq. (1) using the same values of coupling and synaptic time delay used in the results from Fig. 1. In this example, the value of the input, $Istim$ is chosen to be identical to Eq. (39) with the exception of taking $Ts=10.6$ ms, which is closer to the nominal period of 10.9 ms for the coupled oscillations.

Like in the example presented in Fig. 6, we assume that a set of periodic orbits $xp\gamma $ exists for some adaptive parameter $p\u2208R$. We will use the observable $sD$ as defined in Eq. (41) as a proxy for the adaptive parameter and employ the real-time estimation strategy described in Sec. IV B to identify the reduced order terms $H1$ and $H2$ used in the control strategy from Eq. (35). We define $\theta =0$ for each of these orbits to correspond to the moment that $sD$ crosses 0.10 with a positive slope. Once per cycle, when this $\theta =0$, crossing is detected $\varphi $ can be updated using to Eq. (23). Likewise, $sD$ can be updated by computing the difference between the maximum and minimum value of $s\xaf(t)$ on the previous cycle. In-between cycles, the reduced order state variables can be updated according to the reduced order equations from (35) (with $sD$ taking the place of $g$).

As a first step in the implementation, the controller from Eq. (35) is applied taking ${K1,K2,K3,K4}={4,0,0,0}$. This choice of parameters allows for direct control of $\varphi $ (through manipulation of $\zeta $) but does not allow for the control input $\mu $ to change. As such, by setting $\varphi targ$ to be a particular value, $sD$ will approach an associated region in phase space for which $s\u02d9D=0$. Consequently, an approximation for the nullclines of $s\u02d9D=0$ can be obtained by fixing a particular value of $\mu $ and slowly sweeping through different values of $\varphi \u2208[0,2\pi )$. These results are shown in Fig. 7. Taking $\mu =2$, panels (a) and (b) show the value of the order parameter computed according to Eq. (4) and the inferred value of $sD$, respectively; Qualitatively, similar shapes are seen in both traces. Panel (c) shows the nullclines that emerge for various choices of $\mu $. For $\mu $ large enough, setting $\varphi \u22485$ completely disrupts the synchronous oscillations and no nullclines can be identified. Note that simply setting $\varphi =5$ and taking $\mu $ large enough is not a viable long-term desynchronization strategy—when oscillations are lost no information can be obtained for the reduced order model allowing oscillations to re-emerge in the system due to the nominal synaptic coupling between neurons.

Notice from panel (a) of Fig. 7 that the most prominent desynchronization (judging by the order parameter) occurs near values of $\varphi =5$. As such, in a full implementation of the desynchronizing controller from Eq. (35) we choose $\varphi targ=5$ rad with constants ${K1,K2,K3,K4}={4,6,0.18,0.30}$. Results are shown in Fig. 8 for a simulation taking $sD,targ(t)=0.34$ for $t<100$ and $sD,targ(t)=0.08+0.26exp\u2061(\u22120.01(t\u2212100))$ for $t\u2265100$. These values of $sD,targ$ are chosen to be close to the nominal value that results in the absence of control when $t<100$ and to gradually decrease toward lower values (which correspond to less synchronous solutions) over time. For $t<100$, $\mu $ is set to zero and the controller is turned off. The controller is turned on for $t\u2265100$. The black line in panel (a) of Fig. 8 shows the value of $\varphi (t)$ inferred from simulations with the red line corresponding to the target value. Instantaneous updates to $\varphi $ are made once per cycle when the threshold corresponding to $\theta =0$ is crossed, which explains the appearance of short vertical lines in the black trace. Panel (b) shows the actual value of $sD(t)$ superimposed over $sD,targ(t)$. Panels (c) and (d) show the value of the phase offset and the control amplitude, respectively. When the controller is initially turned on, a relatively large value of the input magnitude, $\mu $, is required to initiate the desynchronization; as $sD$ is driven to smaller values, comparatively less control input is required. Panel (e) shows a representative sample of the of voltage traces from the neurons in the oscillator population (colored lines) with the average value of transmembrane voltage shown in black. Disruption of the nominal synchronization is clearly observed. The average value of the synaptic current and the order parameter are shown in panels (f) and (g), respectively. In the absence of input, $r=0.9$ indicating a highly synchronized population. Soon after the input is turned on and transient behavior dies out, the order parameter drops to near 0.2 where it remains indefinitely. Qualitatively, similar results to those shown in Fig. 8 can be obtained by choosing different values of $K1$, $K2$, $K3$, and $K4$.

Robustness of the control algorithm is also considered by varying model parameters. For the same model used to obtain the results from Fig. 8, representative simulations for different choices of the parameter $K1$ are shown in panels (a)–(f) of Fig. 9. When changing one parameter at a time and keeping the rest identical to the parameter set used in Fig. 8, the controller is able to track the specified value of $sD,targ$ for $K1\u2208[3.3,9.0]$, $K2\u2208[3.7,19.5]$, $K3\u2208[0.01,0.29]$, $K4\u2208[0.19,0.39]$, and $\varphi targ\u2208[4.6,5.4]$. These results provide preliminary evidence that the proposed controller is robust to changes in the parameters; however, a more rigorous study of controller robustness would be warranted in future studies. Additionally, alternative implementations of the proposed control framework are also possible. Consider, for instance, taking ${K1,K2,K3,K4}={4,6,0,0}$ and letting the input magnitude to saturate with limits $\mu \u2208[0.1,5]$. Implementing the controller with these parameters does drive the order parameter to lower values, although occasional short-lasting resynchronization sometimes occurs (results are not shown for this alternative parameter set).

Finally, this model independent control strategy is implemented using a modified stimulus,

Instead of a series of only positive pulses, Eq. (42) mandates an alternating set of positive and negative pulses. For a constant value of $\mu $, this stimulation protocol would yield a charge-balanced stimulus, which is a practical requirement in clinical settings to mitigate Faradaic reactions that damage both the surrounding brain tissue and the DBS probe.^{21,33} For a nonstatic choice of $\mu $, the input will be nearly charge balanced. Note that Eq. (42) does not describe a standard biphasic pulse,^{6} but nonetheless will have a substantially smaller accumulation of charge than the monophasic pulses from Eq. (39). Results are shown in Fig. 10 using the pulses mandated by Eq. (42) and taking $\varphi targ=0.5$ rad with all other parameters taken identical to those from Fig. 8. Once again, the controller is turned on at $t=100$ ms. Panels (a)–(c) show $\varphi $, $sD$, and $Istim(t)$ over the course of a representative simulation. Panels (d)–(f) show voltage traces, $s\xaf$, and the order parameter, respectively, over the first 400 ms of simulation. Despite the use of a different input, the results from Fig. 10 are qualitatively similar to those from Fig. 8.

## VI. DISCUSSION AND CONCLUSION

In this work, we develop a closed-loop control framework that can be used to manipulate the behavior of large populations of coupled neurons using near-periodic inputs, i.e., inputs with the general form given in Eq. (2) that are nominally periodic but with the incorporation of a slowly varying phase offset and a slowly changing amplitude. This work can be viewed as potential theoretical underpinning for the phase-specific aDBS application strategies that have gained interest in recent years;^{7,18,49} here, the slowly varying phase offset can be used to coordinate the timing of the DBS stimulus relative to some measured observable. In this work, our primary control objective is to engender desynchronization; however, other control objectives could conceivably be considered using this strategy. Additionally, this general framework could potentially be used in other applications, for instance, in conjunction with non-invasive stimulation therapies^{45} by using accelerometer data from tremor patients as an observable used to define the phase and amplitude of the oscillation.

We employ a series of model transformations and reductions in the analysis performed in this work. Starting with a large population of conductance-based neurons from Eq. (1), we represent the behavior in the large $N$ limit using the Fokker–Planck Eq. (16). The adaptive phase-amplitude reduction paradigm^{67,68} is used to analyze the temporal evolution of the Fokker–Planck equation in terms of a low-order reduction and formal averaging techniques^{16,51} are used to ultimately arrive at an analytically tractable reduced order model of the form (32). A relatively simple controller of form (35) is subsequently proposed to manipulate the state variables of the reduced order model. As illustrated in the results, by simply manipulating the phase offset and the magnitude of the periodic input, the variables of the reduced order model can be controlled successfully.

In comparison to other closed-loop DBS strategies, the proposed method is most similar to phase-triggered methods.^{3,9,19,62,63} Such methods seek to apply a control input at an appropriate phase of oscillation in order to reduce the degree of synchronization in the overall population. In contrast to these other phase-triggered methods, here we use a rotating reference frame [i.e., defined in Eq. (23)] that considers the offset, $\varphi $, between the population phase $\theta $ and the phase of the periodic input. By appropriately manipulating $\varphi $, the population oscillation can be driven to less synchronous states. Also, we note the differences between the proposed method and delayed-feedback approaches detailed in Refs. 40, 47, and 52, which can generally be implemented by applying a feedback signal $Istim(t)=K(y(t\u2212\tau )\u2212y(t))$, where $\tau $ and $K$ are a time delay and a proportional gain, respectively, and $y(t)$ is some observable. With an appropriate observable that provides a good representation of the population level behavior, it is often possible to find some combination of $\tau $ and $K$ that can yield desynchronization. Extensions that use nonlinear delayed-feedback^{39} or charge-balanced pulses^{41} are also possible. In contrast to the method presented in this work, delayed-feedback approaches do not require an *a priori* specification of the applied DBS waveform.

While the terms of the averaged, reduced order model derived in this work can be computed numerically according to the integrals from Eq. (26), they do result after a long series of transformations and reductions. Errors associated with these transformations will compound at each step. Additionally, in general, the underlying model equations can either be unavailable or can be too complicated to numerically compute the terms of the reduction. For this reason, the real-time strategy for estimation of the terms of the reduced order model as described in Sec. IV B would be essential for practical implementation of the proposed control strategy. Indeed, using this strategy, coupled oscillations in the full order model from Eq. (1) can be disrupted without knowledge of the underlying model equations. In the examples considered in this work, we use the average value of the time-delayed synaptic variable as defined in Eq. (3) as a system observable. In experimental applications, the proposed control strategy could potentially be implemented using accelerometer measurements from hand tremors or LFP measurements.

There are many limitations of the proposed control strategy that are not addressed here. For the purposes of this work, parameter sets for the proposed controller are identified by first determining an adequate choice of $\varphi targ$ for which the application of the near-periodic input has a desynchronizing influence. The remaining parameters were obtained using a trial-and-error approach. In future work, it would be of interest to provide a more careful analysis of this control strategy to investigate questions related to controllability and optimality. An additional limitation of the proposed control strategy is that it requires measurable oscillations to be present at all times. As a consequence of this assumption, synchronization cannot be completely abolished. For instance, when $sD,targ$ is set to values below 0.08 in the results shown in Fig. 8, oscillations can occasionally become difficult to identify subsequently degrading the real-time estimates of $\varphi $, $sD$, $H1$, and $H2$. When this happens, the neural population can transiently resynchronize until accurate estimates for these terms can be reobtained and the controller can regain its efficacy. Future work will investigate strategies to mitigate this issue, for instance, by investigating strategies to stabilize the fully desynchronized state. In the examples presented in this work, it is assumed that all but one of the isostable (amplitude) coordinates in Eq. (20) decay rapidly so that they can be ignored. As illustrated in panel (b) of Fig. 2, this assumption is explicitly validated for the discretized version of the Fokker–Planck equation from Eq. (19) for the model used in this study. Numerical illustrations of the proposed control strategy are performed using a single population of identical neurons (1) that do not contain a realistic description of the overall brain circuit that is relevant to clinical DBS. It is likely that models with more complicated coupling structures and models with multiple interacting subpopulations may require the inclusion of additional amplitude coordinates, thereby complicating the analysis. Investigation of the proposed control strategy in a more realistic model would be warranted. Additionally, numerical illustrations in this work use the average value of the synaptic variable as an observable used to define the phase and amplitude of oscillations. In practical applications, such information is not readily available. This would necessitate the use of other biomarkers, such as beta band LFP power. Additional issues created by both measurement noise and delays between the application of stimulation and the effect on system observables would likely need to be considered in practical applications.

## ACKNOWLEDGMENTS

This material is based upon work supported by the National Science Foundation (NSF) (Grant No. CMMI-2140527).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflict of interest to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX: NEURAL MODEL IONIC CURRENTS

In order to model the ionic currents $Iion$ for the population of neurons described by Eq. (1), a three-dimensional model that replicates the behavior of tonically firing thalamic neurons from Ref. 50 is used. The model equations in the absence of coupling, noise, and external stimulus are

where $Iion(V,r,h)$ is comprised of five ionic currents,

each given by

Auxiliary functions are given below:

Additional constants are $C=1\mu F/cm2,gL=0.05mS/cm2,EL=\u221270mV,gNa=3mS/cm2,ENa=50mV,gK=5mS/cm2,EK=\u221290mV,gT=5mS/cm2,andET=0mV$. With these specific parameters, the neuron settles to a tonically firing limit cycle in steady state with a period of 8.40 ms.

## REFERENCES

*2019 IEEE 58th Conference on Decision and Control*(IEEE, 2019), pp. 4717–4722.