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.
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 reset28,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 experimental7,18,49 and computational studies3,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 oscillations63 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 power46 or the severity of tremor27 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 framework67,68 in conjunction with formal averaging techniques16,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, sets the number of neurons in the population, is a constant membrane capacitance, is the transmembrane voltage of neuron (in mV), is a collection of auxiliary variables (gating variables, ionic concentrations, etc.) with dynamics governed by , is a collection of ionic currents, is a synaptic variable where , and determine the dynamics of the synaptic variable, mV is the reversal potential of the neurotransmitter so that the coupling is inhibitory, is a constant time delay, is a constant conductance that sets the synaptic coupling strength, is an externally applied current that is applied identically to each neuron, and is a zero-mean white noise process associated with neuron for which sets the intensity. Here, the collection is independent and identically distributed. The specific forms of and 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 from50 to specify the ionic dynamics. These model equations are given in the Appendix. Terms governing the synaptic variables are taken to be . 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 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 ; the equations governing will take the general form
where is a -periodic input with amplitude and phase offset . Both and are allowed to change slowly according to the relations and , respectively. The input is nearly -periodic, where the magnitude and phase shift are allowed to slowly vary in time. Further structure for the terms and will be considered in the analysis to follow.
Figure 1 shows the behavior of a population of synaptically coupled neurons from Eq. (1) taking the synaptic coupling time delay to be ms. Initially, the coupling strength is set to zero allowing noise to fully desynchronize the population. At ms, the coupling strength is set to 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 is the asymptotic phase of neuron [as defined by the isochrons according to Eq. (8)]. Intuitively, takes values between 0 and 1 and gives a sense of the overall level of phase coherence. 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 . 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 starts near zero and grows to 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 is the system’s state, sets the dynamics, is a collection of nominal parameters, and is an external input. Suppose that Eq. (5) admits a stable, -periodic limit cycle when .
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 is defined on the limit cycle and scaled so that when . 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 and taking , for any initial condition , the isochron corresponding to is defined to be the set of all such that
where gives the solution of (5) with initial condition and is some vector norm. Intuitively, initial conditions on the same isochron have the same asymptotic behavior. Defining the phase using isochrons, one can arrive at a reduced order model by restricting attention to a close neighborhood of the periodic orbit,
where is known as the phase response curve, where each partial derivative is evaluated at and . Note that in the transformation to (7), Taylor expansion is used to obtain as an intermediate step. Phase reduction is a tremendously powerful tool that can be used to analyze oscillations of the -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 where .
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 , one can write , where denotes the Jacobian of the vector field evaluated at . Letting be the fundamental matrix associated with this linear equation [i.e., that solves ], define the left eigenvectors, right eigenvectors, and eigenvalues of as , , and , respectively. Also, define as the Floquet exponent associated with . These Floquet exponents will be sorted so that (corresponding to the eigenvalue of the fundamental matrix) and that for all others. Provided is not a defective eigenvalue,66 the isostable coordinate can be defined in the basin of attraction of the periodic orbit as follows:
where represents the flow of (5) when , is the intersection of the periodic orbit and the isochron, and gives the time of the transversal of the 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
One can verify that Eq. (9) is mathematically equivalent to Eq. (5). Here, is taken to be the nominal dynamical behavior and is an effective additive input. It will be assumed that for some nominal range of system parameters , periodic solutions exists and are continuously differentiable with respect to both and . In reference to each periodic orbit, a generalized phase, and a set of slowly decaying generalized isostable coordinates can be considered. Assuming that and each are continuously differentiable in a neighborhood of the limit cycle for all and , one can transform to phase and isostable coordinates via the chain rule
for . 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, is the natural frequency of the unperturbed periodic orbit , is the gradient of the phase coordinate evaluated on , and captures how changes to the adaptive parameter set influence the generalized phase, is the Floquet exponent associated with the isostable coordinate, is the gradient of the isostable coordinate with respect to the state evaluated on , and captures how changes to the adaptive parameter set influence the isostable coordinates, and 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 in order keep the isostable (amplitude) coordinates small, thereby keeping the state close to the periodic orbit . To accomplish this, must be chosen appropriately. Some general design heuristics for are proposed in Ref. 68. For instance, when only a single non-truncated isostable coordinate is considered and , it often works well in practice to choose
taking . For higher dimensional reductions, must be designed on an ad hoc basis. Provided an appropriate can be found so that each remains an 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 and 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 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 . Here, is the phase of the neuron. Because the periodic orbit is strongly attracting, functions and are used to represent the transmembrane voltage and synaptic variable associated with the neuron as a function of the phase.22 Additionally, . As part of the transformation to phase coordinates, it is assumed that is small enough so that terms proportional to the square of the noise intensity can be truncated26 (cf., Ref. 2). In the formulation above, and are assumed to be small enough so that is well approximated by for all . 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 that evolves according to the Fokker–Planck equation using an Ito formulation,13
and and denote the partial derivatives taken with respect to . Because , 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 . 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 . The PDE model (16) will be further analyzed using a finite difference approximation with 100 total grid points for the axis. Letting represent the discretized solution of the phase distribution, after discretization the dynamics of Eq. (16) can be approximated according to
where is a discretized version of , represents the resulting dynamics that are solely a function of the state and coupling strength, is a matrix which implements the finite difference approximation of the derivative of with respect to , and 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 taking the place of ], 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 is taken to be the adaptive parameter. For each value of considered, the model has a stable periodic orbit 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 , the phase of an individual oscillator with respect to its unperturbed limit cycle. Each periodic orbit, , will be aligned so that corresponds to the moment that 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 and . Above, the structure of and has been rewritten to depend on the reduced order model parameters , , 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 from Eq. (20) matches the general strategy from (13) that serves to drive to smaller magnitude values. Figure 2 shows various characteristics of the reduced order model from (20). Panel (a) shows plotted at for various values of on the appropriate limit cycles. As 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 , in response to a the periodic input , 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 are associated with more (respectively less) phase cohesive oscillations. Panel (e) shows that the isostable coordinate 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 tends toward infinity. Overall, we observe good agreement between the full order and reduced order simulations.
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 from Eq. (20) is periodic in , we illustrate that further reduction in dimensionality can be realized using formal averaging techniques.16,51
To begin, recalling from Eq. (2) that is the nominal period of the input , we will work in a rotating reference frame defining
where . 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 is periodic in time with period . We will assume that , , and are terms so that Eq. (24) can be written in the general form where , and gives the dynamics and is periodic in . For such a system, the method of averaging16,51 can be used to approximate the behavior of (24) according to
Using the averaging framework, Eq. (25) provides a close approximation to the unaveraged system of equations from (24). Note that the terms through 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, through are assumed to be terms; even though and are multiplied by , the term in the integrand is an term.
Further simplification of Eq. (25) is possible provided is chosen large enough so that rapidly converges to its steady state value. To this end, define
as the steady state value of that results when , and are held fixed. Recognizing that the right hand sides of Eq. (25) are comprised of terms, , and are also terms. Taking the time derivative of one can show that
Next, defining and taking the time derivative, one finds
Because all periodic orbits used in the adaptive reduction are stable, . Likewise, since the integrand from (26) is non-negative and by definition. Considering that provided that at , uniformly bounded in time and, consequently, can be used as an approximation that is valid up to . This allows for the elimination of the dynamics from Eq. (25) to yield
Finally, provided , one can write
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 and , respectively. Strategies for designing these inputs will be considered momentarily. Note that in situations where , 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 is sufficiently larger than for all choices of and .
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 so that the input is fully periodic. Note here that for the input , when computing the averaged equation is taken to be in order to keep the magnitude of small. After computing the averaged functions – using (26), (respectively, ) 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 and 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.
A. A controller for the input magnitude and phase offset
Our goal is to identify functions and from Eq. (20) that can shift the state of the adaptive parameter in the reduced order equations from (20) from to some . By shifting 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 that drives the system to less synchronized states that are near the periodic orbit associated with the parameter .
Intuitively, toward the identification of an appropriate control strategy, first consider the nullcline from panel (c) of Fig. 3 obtained for the averaged equations (32). Referring to this nullcline by , for any in the range of , the control objective could be satisfied provided that could be set to some value for which . For instance, considering the nullcline from panel (c) of Fig. 3, if a control strategy can be found that sets to a constant value of 5 rad, 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, can be driven to any desired value less then . In this case, however, the magnitude of the input, would need to be decreased as is driven near its lower limit. With this in mind, we will consider a control strategy that seeks to drive and to and , respectively. This control strategy will be implemented by mandating the following update rules for and from Eq. (20):
where – are appropriately chosen, positive constants. Here, in the averaged reference frame, i.e., using (32), this choice of and yields the relationships
where is an extra variable used to represent the necessary terms from the integral control. Intuitively in the above equations, the term acts as a proportional controller that serves to update the offset of the rotating reference frame, , to manipulate . The term serves as a proportional–integral controller that increases or decreases the magnitude of the applied input depending on whether is above or below its target value. The additional term, , is added to influence stability of the controller. Of course, the efficacy of the controller (34) depends on an adequate choice for and . 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 and . 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 and 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 can be clearly identified from measured data [for instance, for a given periodic orbit from Eq. (16) corresponds to the moment that reaches a local maximum]. Consider two successive traversals of the level set, occurring at and . The associated coordinates in the rotating reference frame defined by Eq. (23) are
Toward identifying a strategy to estimate the unknown functions and from observable data, one can rewrite Eq. (32) as
Above, the arguments of and are evaluated at the fixed time . Because , , , , and are terms, fixing time in the evaluation of and only contributes an additional term when . Neglecting the and terms, through direct integration of Eq. (37), one can write
Equation (38) can be used to obtain successive, real-time estimates of and , for instance, when implementing the control strategy suggested by Eq. (35), is known exactly at all times because it is part of the controller, and can be computed according to (23) at successive crossings of the level set, and and 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, and 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 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 to be used for the implementation of the control strategy. This observable should have clearly identifiable oscillations. Define 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 and to be the times of subsequent crossings of the isochron, define a map from the observation of on the interval to the variable that encodes for the oscillation amplitude. For instance, one could take on the interval .
- (Step 3)
Choose a periodic stimulation to be applied to the system. The period of must be close to the period of oscillation of the observable chosen in step 1.
- (Step 4)
Choose a set of constants , and for the controller from Eq. (35). In a model independent setting, this will generally be a trial and error process. One can start with and to identify choices of that effectively decreases the magnitude of oscillations. The terms , , and can subsequently be tuned to appropriately modulate the stimulus magnitude.
- (Step 5)
Choose a profile for , which sets the target for the oscillation amplitude.
- (Step 6)
Implement the controller from Eq. (35). On a cycle-by-cycle basis, estimates of and can be updated according to Eq. (38), estimates of and can be computed according to (23) at successive crossings of the level set, and and can be inferred from the mapping defined in Step 2.
- (Step 7)
Adjust constants , and as necessary. Performance is generally less sensitive to the choice of these constants when is close to the nominal value of that results in the absence of . As is decreased (i.e., as the oscillation amplitude decreases) parameter adjustments may be necessary to yield a viable controller.
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 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 and can be computed directly from the terms of the reduced order equations using the averaging framework. For the full model of coupled neurons, direct computation of and 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 . Recall that for the allowable range of coupling strengths , the associated periodic orbits have periods ms. For this application, we consider the nominal input to be a periodic series of pulses,
Here, ms and the input is a series of six pulses shown taking 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 shown at different values of . Relative to the baseline at , tends to decrease when and tends to increase for . With this in mind, recalling that our goal is to decrease the phase cohesion by driving the adaptive parameter, , to lower values, to implement the control strategy from Eq. (35), we choose rad. Three different parameter sets and are chosen to implement the control strategy taking , , and , respectively. Parameter set uses only two proportional controllers to determine the update rules for and . Parameter set adds an integral controller to the update rule for , and parameter set adds an additional term that is proportional to . 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 variable. Here, is ignored to account for the fact that might not fully reach its target value. Fixed points of Eq. (40) are obtained with a Newton iteration for various choices of and plotted against eigenvalues with the largest real component, , of the linearization about the fixed point. For parameter sets and , stability is lost for lower values of as crosses from negative to positive values. Panel (d) illustrates the control strategy using parameter sets and for simulations of the reduced order equations from (35) taking . The parameter set is able to follow this target accurately while and display unstable behavior for small values of . 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 ms and take rad. Like in the previous example, and 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 [corresponding to the moment that from Eq. (17) achieves a local maximum] both and are inferred and updated based on the measurements of . For each of these updates, is inferred from Eq. (23) and is inferred by comparing the value of to the maximum values of attained for each periodic orbit . Results are shown in Fig. 5. Parameter sets , and are identical to those from Fig. 4 with the exception of taking instead of 35. Parameters , , and are shown over the course of each simulation in panels (a)–(c) and the adaptive parameter is shown in panel (d). In this example, for and for . Results are qualitatively similar to those from Fig. 4 with parameter sets and 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 is taken to be the adaptive parameter, periodic orbits do not exist when taking . This represents a challenge when the overall goal is desynchronization, since the Kuramoto order parameter associated with the periodic orbit that emerges for 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, , exists for some adaptive parameter, . Once again, these orbits will be aligned so that corresponds to the moment that from Eq. (17) reaches a local maximum. We will assume the existence of some model observable such that . 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 on the periodic orbit . This definition allows for the inference from measurements of once per cycle. Note that 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 and on a cycle-by-cycle basis. In this example, takes the place of the adaptive parameter from (35). For the control parameters, we use and take for and for . These values of are chosen to be close to the nominal value that results in the absence of control for and to gradually decrease toward lower values (which correspond to less synchronous solutions) over time. Controller parameters , , , and are taken to be 4, 30, 1, and 1, respectively. Qualitatively similar results can be obtained using other values for through . Results are shown in Fig. 6. Panels (a)–(c) show , and over the course of the simulation. These traces are qualitatively similar to those obtained using the 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 and 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 during the simulation, which gradually decays toward smaller values, as mandated by the chosen value of . 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 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, is chosen to be identical to Eq. (39) with the exception of taking 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 exists for some adaptive parameter . We will use the observable 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 and used in the control strategy from Eq. (35). We define for each of these orbits to correspond to the moment that crosses 0.10 with a positive slope. Once per cycle, when this , crossing is detected can be updated using to Eq. (23). Likewise, can be updated by computing the difference between the maximum and minimum value of on the previous cycle. In-between cycles, the reduced order state variables can be updated according to the reduced order equations from (35) (with taking the place of ).
As a first step in the implementation, the controller from Eq. (35) is applied taking . This choice of parameters allows for direct control of (through manipulation of ) but does not allow for the control input to change. As such, by setting to be a particular value, will approach an associated region in phase space for which . Consequently, an approximation for the nullclines of can be obtained by fixing a particular value of and slowly sweeping through different values of . These results are shown in Fig. 7. Taking , panels (a) and (b) show the value of the order parameter computed according to Eq. (4) and the inferred value of , respectively; Qualitatively, similar shapes are seen in both traces. Panel (c) shows the nullclines that emerge for various choices of . For large enough, setting completely disrupts the synchronous oscillations and no nullclines can be identified. Note that simply setting and taking 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 . As such, in a full implementation of the desynchronizing controller from Eq. (35) we choose rad with constants . Results are shown in Fig. 8 for a simulation taking for and for . These values of are chosen to be close to the nominal value that results in the absence of control when and to gradually decrease toward lower values (which correspond to less synchronous solutions) over time. For , is set to zero and the controller is turned off. The controller is turned on for . The black line in panel (a) of Fig. 8 shows the value of inferred from simulations with the red line corresponding to the target value. Instantaneous updates to are made once per cycle when the threshold corresponding to is crossed, which explains the appearance of short vertical lines in the black trace. Panel (b) shows the actual value of superimposed over . 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, , is required to initiate the desynchronization; as 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, 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 , , , and .
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 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 for , , , , and . 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 and letting the input magnitude to saturate with limits . 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 , 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 , 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 rad with all other parameters taken identical to those from Fig. 8. Once again, the controller is turned on at ms. Panels (a)–(c) show , , and over the course of a representative simulation. Panels (d)–(f) show voltage traces, , 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 therapies45 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 limit using the Fokker–Planck Eq. (16). The adaptive phase-amplitude reduction paradigm67,68 is used to analyze the temporal evolution of the Fokker–Planck equation in terms of a low-order reduction and formal averaging techniques16,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, , between the population phase and the phase of the periodic input. By appropriately manipulating , 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 , where and are a time delay and a proportional gain, respectively, and 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 and that can yield desynchronization. Extensions that use nonlinear delayed-feedback39 or charge-balanced pulses41 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 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 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 , , , and . 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.
This material is based upon work supported by the National Science Foundation (NSF) (Grant No. CMMI-2140527).
Conflict of Interest
The authors have no conflict of interest to disclose.
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 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 is comprised of five ionic currents,
each given by
Auxiliary functions are given below:
Additional constants are . With these specific parameters, the neuron settles to a tonically firing limit cycle in steady state with a period of 8.40 ms.