Abstract
In this paper, we calculate magnitude-constrained optimal stimuli for desynchronizing a population of neurons by maximizing the Lyapunov exponent for the phase difference between pairs of neurons while simultaneously minimizing the energy which is used. This theoretical result informs the way optimal inputs can be designed for deep brain stimulation in cases where there is a biological or electronic constraint on the amount of current that can be applied. By exploring a range of parameter values, we characterize how the constraint magnitude affects the Lyapunov exponent and energy usage. Finally, we demonstrate the efficacy of this approach by considering a computational model for a population of neurons with repeated event-triggered optimal inputs.
1 Introduction
As deep brain stimulation (DBS) emerges as an effective therapy for a wide range of neurological disorders, theoretical perspectives are growing to help inform effective stimulation protocols. DBS technology involves surgically implanting electrodes to deliver electrical stimuli over time to specific brain regions with wires connecting the implanted electrodes to an implantable pulse generator, which can be programmed to set DBS parameters (; ; ). Novel technologies enable the optimization of DBS stimulation parameters including the input pulse’s shape, amplitude, frequency, and interstimulus interval (; ). Computational analysis of the interaction between an applied electrical stimulus and the spiking dynamics of neural populations can inform the effective design of DBS parameters to enhance clinical outcomes while respecting device engineering constraints.
Among other conditions, deep brain stimulation has become an effective treatment for Parkinson’s disease, where pathological synchronization of the basal ganglia-cortical loop is associated with dopaminergic denervation of the striatum, which over time leads to motor impairment including tremors, bradykinesia, and akinesia (; ). In particular, the parkinsonian low-dopamine state is observed to be related to excessive synchronization in the beta frequency band (15–30 Hz) in the subthalamic nucleus (STN) of the basal ganglia (; ). It has been proposed that the symptoms of parkinsonian resting tremors are caused by excessive synchronization in populations of neurons firing at a similar frequency to that of the tremor (). Deep brain stimulation works as a therapeutic intervention by modulating such synchronization patterns, often ameliorating motor impairment for patients living with Parkinson’s disease. DBS is also used in the treatment of essential tremor (ET), epilepsy, Tourette’s syndrome, obsessive compulsive disorder, and treatment-resistant depression.
The most commonly offered protocol for DBS therapy is continuous high-frequency stimulation (; ). Despite its clinical success, conventional high-frequency DBS faces several limitations including diminishing efficacy over time, stimulation-induced side effects, and high energy consumption necessitating battery replacements. Furthermore, as open-loop stimulation, where the device continuously applies electrical inputs as long as it is on, its static parameters do not adapt to the dynamic nature of disease symptoms, limiting its long-term effectiveness. In recent years, there has been growing interest in developing alternative stimulation paradigms that can effectively disrupt pathological synchrony while minimizing energy consumption and side effects. Several desynchronization methods have been proposed and tested on patients, including coordinated reset stimulation (; ), adaptive deep brain stimulation (aDBS) (), and phase-specific stimulation approaches (). Such techniques often leverage theoretical frameworks from nonlinear dynamics and control theory to design stimuli that can efficiently desynchronize neural populations.
In particular, coordinated reset uses multiple electrode implants which deliver identical impulses separated by a time delay between implants (; ; ; ; ; ). This leads to clustering behavior for the neural populations, in which each cluster fires at different times, giving (partial) desynchronization of the dynamics. This approach has achieved preliminary clinical success (; ). In adaptive deep brain stimulation (aDBS), closed-loop systems monitor biomarkers in real time and can initiate changes in stimulation parameters (; ). The goal of aDBS is to improve clinical outcomes by designing control signals based on neural data to deliver stimulation only when needed. For instance, in Parkinson’s disease, a potential biomarker is the amplitude of the pathological beta rhythms, with stimulation becoming active if this exceeds some prescribed threshold. On demand stimulation through aDBS can prolong device battery life, thereby extending device lifetimes and prolonging the interval between surgeries in clinical treatment protocols. Finally, for phase-specific stimulation the inputs occur at a particular dynamical phase in order to disrupt synchrony (; ; cf. ).
In parallel, there have been a number of computational studies exploring the mechanisms by which DBS might be working (e.g., ; ). There have also been theoretical and computational studies exploring different strategies for desynchronizing neural populations (), with approaches including delayed feedback control (; ; ; ), phase randomization through optimal phase resetting (; ; ), phase distribution control (; ), cluster control (; ; ; ), machine learning and data-driven approaches (; ), and chaotic desynchronization (; ; ). Each of these approaches has advantages and disadvantages based on the control objective and what is known about and what can be measured for the neural dynamics ().
In this paper, we focus on chaotic desynchronization, for which an energy-optimal stimulus exponentially desynchronizes a population of neurons. This approach relies on phase reduction methods, which have proven particularly valuable in analyzing and controlling neural oscillators (; ). These methods allow for the simplification of complex neuronal dynamics into phase models, where the behavior of an oscillating neuron can be characterized by its phase and response to perturbations, captured by the phase response curve (PRC). Unlike previous studies of chaotic desynchronization, here we include a constraint on stimulus magnitude. Such constraints are important engineering considerations for the practical applications of DBS, as there can exist both biological limitations on the maximum electrical stimulation that can be safely applied to brain tissue as well as electronic limitations on the current that stimulation devices can reliably store and deliver over time.
Specifically, in this paper we calculate magnitude-constrained optimal stimuli that maximize the Lyapunov exponent for the phase difference between pairs of neurons while simultaneously minimizing energy consumption. The Lyapunov exponent quantifies the exponential divergence rate of initially close trajectories, making it an appropriate measure of desynchronization efficiency. By systematically exploring different constraint magnitudes, we characterize the tradeoff between maximum allowable stimulus amplitude, desynchronization efficacy, and energy usage. We set up the optimal control problem in Section 2.1. In Section 2.2 we describe several canonical phase response curves representing different types of neuronal dynamics: Sinusoidal, SNIPER, Hodgkin-Huxley, and Reduced Hodgkin-Huxley models. In Section 3.1, we investigate our approach for each PRC, computing the optimal stimulus under various magnitude constraints and evaluating its performance in desynchronizing initially synchronized neurons. In Section 3.2, we validate our approach using computational simulations of neural populations with coupling and noise, demonstrating that our magnitude-constrained optimal stimuli can effectively desynchronize neural populations. Finally, a discussion of our results is given in Section 4. Overall, this work provides a theoretical foundation for designing energy-efficient DBS protocols that respect hardware and biological constraints while effectively disrupting pathological neural synchronization. We respectfully present this study as a tribute to the pioneering work of Hermann Haken on the control of complex systems.
2 Methods
2.1 Optimal control problem
We present a procedure for finding an energy-optimal stimulus which maximizes the Lyapunov exponent associated with the phase difference between a pair of neurons, while accounting for a constraint on the stimulus magnitude. This approach is based on the phase reduction of neural oscillators in the presence of an input (see, for example, ), and only requires knowledge of a neuron’s phase response curve (PRC). We note that the PRC can in principle be measured experimentally (), or can be calculated numerically if the model is known (; ). In particular, we consider the following set of equations:where is the phase of the neuron; , where is the period of the neuron in the absence of stimulus; and is the control stimulus. Note that here we are assuming that the neurons are identical (having the same and ), and for simplicity we assume that the neurons are the same distance from the electrode so they receive the same stimulus . Neuron fires an action potential when crosses through 0.
Following , we suppose that the neurons are nearly synchronized . Defining , we obtain
Linearizing about , the solution to Equation 2 is , where
Here can be viewed as the finite time Lyapunov exponent (), which characterizes the exponential growth or decay of the phase difference . A positive Lyapunov exponent will correspond to the divergence of nearby trajectories, and hence desynchronization. We note that this Lyapunov exponent corresponds to phase difference direction, so it is not directly related to the non-trivial Floquet multipliers which describe transverse stability of the periodic orbits for the neurons (). We formulate the control problem in terms of the cost functionwhere the goal is to maximize the Lyapunov exponent while minimizing the energy used, where the energy is the integral of the square of the control stimulus . Here is the time that we choose for the duration of the stimulus, is a parameter that scales the importance of the Lyapunov exponent term relative to the energy term. Generalizing the formulation in , here we consider a magnitude constraint on the control stimulus given by Equation 5:
In order to account for this constraint, we use a Hamiltonian formulation for the optimal control problem (), with the Hamiltonian given in Equation 6:where is the Lagrange multiplier or co-state for the system. From Hamilton’s equations,
This defines a two-point boundary value problem which must be solved subject to the boundary conditions and . The latter boundary condition ensures that the phase at time is the same as what it would’ve been in the absence of stimulus. The function in these equations will be found using Pontryagin’s minimum principle (), which states that should be chosen as the extremum of the Hamiltonian, subject to the constraints. If there is no constraint on the magnitude of , the optimal control stimulus is the solution to , giving Equation 9:where is the optimal unconstrained input. With constraints, Pontryagin’s minimum principle gives the following expression for the optimal magnitude-constrained input :
In particular, the optimal might or might not be the solution to , because the extremum may be reached at a constraint boundary. To summarize, we solve the two-point boundary value problem Equations 7, 8 using Equation 11. This is done numerically using a shooting method, and the optimal stimulus is given by Equation 11.
2.2 Example phase response curves
The PRC quantifies the effect of an external stimulus on the phase of a periodic orbit. In this paper, we consider four example PRCs: Sinusoidal, SNIPER, Hodgkin-Huxley, and Reduced Hodgkin-Huxley, as shown in Figure 1.
FIGURE 1
The Sinusoidal PRCshown in Figure 1A with , is a special case of the PRC which is found for periodic orbits close to a supercritical Hopf bifurcation. (Recall that when a parameter is on one side of a supercritical Hopf bifurcation there is a stable fixed point and no periodic orbit, and when the parameter is on the other side of the supercritical Hopf bifurcation there is an unstable fixed point and a stable periodic orbit). More generally, the PRC for a periodic orbit close to a supercritical Hopf bifurcation is a phase-shifted form of the PRC (Equation 12); see . This is an example of a Type II PRC (), in which the PRC takes both positive and negative values.
Periodic orbits can also arise from a SNIPER bifurcation, which stands for Saddle-Node Infinite Period bifurcation; this is also often called a SNIC bifurcation, which stands for Saddle-Node Invariant Circle bifurcation. Here, for a parameter on one side of the bifurcation there is a stable fixed point and a saddle fixed point that lie on an invariant circle. As the parameter is varied, these fixed points annihilate in a saddle-node bifurcation, and when the parameter is on the other side of the bifurcation there is a stable periodic orbit whose period approaches infinity as the bifurcation is approached. For a periodic orbit near a SNIPER bifurcation, the PRC is approximately given by Equation 13 (; ):shown in Figure 1B with . This is an example of a Type I PRC (; ), in which the PRC takes only non-negative values.
The Hodgkin-Huxley equations are a well-studied conductance-based model for neural activity, and were developed to describe the dynamics for a squid giant axon (). Mathematically, they are a four-dimensional set of coupled ordinary differential equations for the voltage across the neural membrane and three gating variables associated with the flow of ions across the membrane. The full equations are given in the Supplementary Appendix. We chose a baseline current value so that the Hodgkin-Huxley equations have a stable periodic orbit, and then found the PRC for this periodic orbit using XPP (). For computational convenience, we approximate this as a Fourier series with the first ten and terms to give the PRC shown in Figure 1C.
Finally, the Reduced Hodgkin-Huxley equations are an approximation to the full Hodgkin-Huxley equations (; ). Mathematically, they are two-dimensional set of coupled ordinary differential equations for the voltage across the neural membrane and one gating variable. The equations are given in the Supplementary Appendix. We chose a baseline current value so that the Reduced Hodgkin-Huxley equations have a stable periodic orbit, and then found the PRC for this periodic orbit using XPP (). For computational convenience, we approximate this as a Fourier series with the first two hundred and terms to give the PRC shown in Figure 1D. The PRCs for the Hodgkin-Huxley and Reduced Hodgkin-Huxley models are examples of Type II PRCs.
3 Results
3.1 Results for pairs of neurons
In this section, we consider the dynamics of a pair of neurons satisfying Equation 1, where is chosen to be the optimal control stimulus for different values of the constraint . For simplicity, we will take , so that the duration of the control stimulus is equal to the period of the neuron in the absence of stimulus. By design, we expect that application of one cycle of the optimal stimulus will cause the phase difference to increase, at least when the initial value for is small. We will consider an event-based approach for which multiple cycles of the optimal control stimulus are applied, where a new cycle of the control stimulus is triggered when , that is, when Neuron 1 fires an action potential. We will see that this leads to growing phase difference .
Figure 2 shows results for the Sinusoidal PRC with , and , corresponding to . In particular, Figure 2A shows the calculated optimal control stimuli for the unconstrained case, , and . In this case, adding a magnitude constraint gives an optimal control stimulus which appears to be very similar to the unconstrained stimulus “chopped off” at the constraint; we will discuss this further below. We observe that the unconstrained input has the highest efficacy of desynchronization as measured by the Lyapunov exponent, and also the highest energy utilization. The input with a constraint gives desynchronization results that are similar while achieving a significant reduction in energy usage.
FIGURE 2
To investigate the efficacy of desynchronization between Neurons 1 and 2, we ran simulations of the phases of Neurons 1 and 2 over a full cycle of the control stimulus, with initial conditions and . Figure 2B shows the time series traces of these two neurons for the Sinusoidal PRC. We see growing desynchronization over one cycle of the control stimulus, as measured by the phase difference between the two traces.
Next, we apply successive control stimuli to pairs of neurons with the event-based approach described above. We compute the phase difference between the pair of neurons, which is observed to grow exponentially over multiple cycles of the control stimulus; see Figure 2C. This is also evident in Figure 2D, where we observe that approximately grows linearly with . To quantify this, we estimate the Lyapunov exponent based on the line of best fit to versus over multiple event-triggered control stimuli. We used the first half of the time interval for these fits, since there is saturation in these traces in the later part of time interval as the small approximation used to obtain Equation 3 no longer holds. These estimated Lyapunov exponents are plotted at varying values of the constraint in Figure 2E, along with estimates obtained by plugging the computed stimulus into Equation 3 and numerically evaluating the integral over one cycle of the stimulus according to Equation 14:
Table 1 compares and for several values of ; good agreement is found between these approaches. Finally, Figure 2F shows the energyused over one cycle of the control stimulus for a range of values of .
TABLE 1
| 0.5 | 0.114 | 0.126 |
| 0.35 | 0.0937 | 0.102 |
| 0.2 | 0.0583 | 0.0621 |
For the Sinusoidal PRC, comparison of the Lyapunov exponent estimated from slope of the line of best fit for versus with the Lyapunov exponent calculated from the integral formulation.
Similarly, Figure 3 shows results for the SNIPER PRC with , , and , corresponding to , and Table 2 compares and for several values of . Moreover, Figure 4 shows results for the Hodgkin-Huxley PRC with and period , which was obtained numerically; Table 3 compares the Lyapunov exponent estimates. Finally, Figure 5 shows results for the Reduced Hodgkin-Huxley PRC with and period , which was obtained numerically; Table 4 compares the Lyapunov exponent estimates. In all of these examples, the optimal input gives exponential divergence of the phases of the neurons. As becomes smaller, the Lyapunov exponent becomes smaller while staying positive, and the energy associated with the input stimulus is reduced. We observe that the numerically calculated optimal inputs without the magnitude constraint resemble for each PRC. This was first noticed in , where it was attributed to the numerical observations that the optimal input is weak enough that , and dominates in Equation 10, so . This approximation is explored in more detail in .
FIGURE 3
TABLE 2
| 0.5 | 0.122 | 0.123 |
| 0.35 | 0.103 | 0.101 |
| 0.2 | 0.0626 | 0.0618 |
For the SNIPER PRC, comparison of the Lyapunov exponent estimated from slope of the line of best fit for versus with the Lyapunov exponent calculated from the integral formulation.
FIGURE 4
TABLE 3
| 0.4 | 0.0239 | 0.0243 |
| 0.3 | 0.0221 | 0.0222 |
| 0.2 | 0.0163 | 0.0172 |
For the Hodgkin-Huxley PRC, comparison of the Lyapunov exponent estimated from slope of the line of best fit for versus with the Lyapunov exponent calculated from the integral formulation.
FIGURE 5
TABLE 4
| 2.5 | 0.219 | 0.227 |
| 1.5 | 0.160 | 0.172 |
| 0.5 | 0.0655 | 0.0625 |
For the Reduced Hodgkin-Huxley PRC, comparison of the Lyapunov exponent estimated from slope of the line of best fit for versus with the Lyapunov exponent calculated from the integral formulation.
Here we make the new observation that in some cases the magnitude-constrained optimal input resembles the unconstrained input simply “chopped off” at the constraint, i.e., , where here would be the optimal input found by solving Equations 7, 8 using Equation 10, i.e., without any magnitude constraint. That said, we note that is not necessarily a good approximation to the optimal input. Solutions to the two-point boundary value problem have the property of modifying a neuron’s phase to go from to in the time . This must be the case for the unconstrained optimal input and for the constrained optimal input . Because is different from during the time intervals for which the constraint is applied, but otherwise the same, we do not expect it to exactly take from 0 to in the time . However, numerically we find for some examples that looks very similar to . This appears to be because the product is very small in these examples. But, it is clear from Figure 5 that for the constraint the optimal can be significantly different from .
When is large, the Lyapunov exponent term in Equation 4 dominates the energy usage term. Therefore, in the limit of large the optimal inputs approach what is known as Bang-bang control, as seen in Figure 6. For Bang-bang control, the controller alternates between maximum and minimum inputs which switch at optimal times (). In particular, the control is driven by the sign of : if we take , and if we take . Thus, at every time the integrand is maximized subject to the constraint, so the value of the Lyapunov exponent is maximized. The unconstrained stimulus magnitudes rise with increasing , so in the large limit the magnitude constraint becomes more important. Figure 6 shows the approach to Bang-bang control for increasing (from left to right panel for each row) on the optimal input for each PRC.
FIGURE 6
3.2 Results for population-level simulations of neurons
While the results in the previous section illustrate that the stimuli with and without magnitude constraints can give positive Lyapunov exponents for the phase difference betweeen pairs of neural oscillators, we are also interested in how such inputs perform for a larger population of coupled neural oscillators. In this section, we consider a population of Reduced Hodgkin Huxley neurons with all-to-all electrotonic coupling, and independent additive noise for each neuron. The governing equations are:where is the number of neurons, is the common input for all neurons, is the coupling strength, and is intrinsic noise modeled as zero-mean Gaussian white noise with variance , where . These values are chosen so that there is a balance between the synchronizing influence of the coupling and the desynchronizing influence of the noise. Expressions for and are given in the Supplementary Appendix. We suppose that at all neurons have the same and corresponding to the neuron at the peak of its action potential. Each neuron receives the same input , but a different realization of noise . It is useful to think of the noise as having a desynchronizing effect, and the coupling as having a synchronizing effect.
To test our magnitude-constrained optimal inputs , we use an event-based control scheme similar to , , , . In particular, when the average voltage for the neurons crosses a threshold, we input one cycle of the pre-computed optimal stimulus, which will have a desynchronizing influence. Another control input occurs if the previous input has finished and the average voltage again crosses threshold, for example, due to the synchronizing influence of the coupling. According to this control logic, each simulation generates a control input for the full 350 time window. Figures 7–9 respectively show results for population-level simulations for the optimal unconstrained input, optimal constrained input with , and optimal constrained input with . Comparing the results using event-based control with the network’s behavior without control, it is apparent that all of these stimuli are able to keep the neural population desynchronized. We can interpret these results as follows: if the population is too synchronized (its average voltage goes above the control activation threshold), a cycle of input is applied. If this does not sufficiently desynchronize the population, another cycle of input is applied. If the population is sufficiently desynchronized, no input is needed. However, the coupling eventually leads to a level of synchronization which is above the control activation threshold, triggering another cycle of input. For it is apparent that more cycles of the stimulus are needed to achieve and maintain desynchronized dynamics; see Figure 9.
FIGURE 7
FIGURE 8
FIGURE 9
This motivates us to better understand how the average amount of energy required to keep the neural population desynchronized with this event-based scheme depends on the magnitude constraint. To investigate this, we consider 100 different population-level simulations for different values of . For each simulation, the neurons start out completely synchronized. Results are shown in Figure 10 for the average energywhere denotes the average over the 100 simulations, and the integration is over the first 350 . The error bars represent the standard deviation over the 100 iterations of the population-level simulations. The variability is due to the different realizations of noise for each simulation. We see that the average energy needed to desynchronize the population increases monotonically with , until the constraint no longer has any effect on the computed input.
FIGURE 10
We used the average voltage as a measure of synchronization because this is easily defined, and one expects this to be related to the local field potential, which can be measured experimentally. However, a more common measure of synchronization is the Kuramoto order parameter (), whose amplitude iswhere is the number of neurons, and is the phase of the neuron in the sense of isochrons (). In particular, the phase is defined at all points in the basin of attraction of the periodic orbit; this is important because to calculate the Kuromoto order parameter we need to know the phase even if noise, coupling, and/or a control input causes the state of the neuron to be off the periodic orbit. We define (which is equivalent to ) to be the phase at which the neuron spikes, and parametrize the isochrons so that in the absence of noise, coupling, and control input, where is the period of the neuron. We can estimate the phase of a neuron by setting its current state as the initial condition at for the Reduced Hodgkin-Huxley equations for a single neuron in the absence of noise, coupling and control input. We integrate these equations forward in time until a voltage spike occurs at time . Because the spike corresponds to and the phase advances according to , the phase corresponding to the state is given by Equation 16:
To obtain the order parameter at a given time, we estimate the phases of all neurons at that time, and use these in Equation 15. Higher values of correspond to greater synchronization.
Figure 11 shows results from this order parameter calculation for the plots shown in Figure 7 (no constraint on the magnitude of ), Figure 8 (for ), and Figure 9 (for ), at intervals. We see that, broadly speaking, the control inputs tend to reduce , and that increases when there are no control inputs, due to the synchronizing effect of the coupling. The results for the no constraint case and for are quite similar except for later times. The fact that increases more rapidly in the no constraint case than for is due to the variability in results due to noise. More interestingly, we see that the order parameter tends to have higher values for ; because the control input is more constrained, a single cycle of input has a smaller effect on the order parameter.
FIGURE 11
4 Discussion
Motivated by deep brain stimulation treatment of neurological disorders including Parkinson’s disease, there has been much recent interest in using control theory to design optimal input stimuli which desynchronize neural populations. One such approach used chaotic desynchronization to achieve this goal in an energy-optimal fashion (). In this paper, we generalized this optimal chaotic desynchronization methodology by including a magnitude constraint on the input. In particular, we showed how to calculate the optimal stimuli which satisfy such a contraint, and demonstrated that such inputs lead to exponential desynchronization of pairs of neurons and effective dysynchronization of populations of coupled, noisy neurons. This approach is based on a phase-reduction of the dynamics of a neuron, with the phase response curve quantifying the effect of an external input on a neuron’s phase. Based on the knowledge of a neuron’s phase response curve, we calculated the inputs that are optimally desynchronizing while minimizing energy utilization, using methods from optimal control. The addition of a magnitude constraint allows for the design of optimal stimulation inputs with a maximum amplitude that is customizable based on biological and electronic considerations. Interestingly, while these constrained inputs use less energy, they still achieve population-level desynchronization.
We note an extension of the current paper which could be of interest in future work: incorporation of both magnitude and charge-balance constraints on the control stimulus. This follows from the observation that non-charge-balanced stimuli, such as those considered in this paper, can cause harmful Faradaic reactions that may damage the DBS electrode or neural tissue (). This has motivated the use of a charge-balance constraint for optimal control design (; ). However, this presents additional challenges because it increases the dimension of the two-point boundary value problem which must be solved numerically. Our approach could be applied to other types of neurons, including those for common neurostimulation targets such as the subthalamic nucleus or the thalamus, even for periodically bursting neurons, as long as the phase response curve can be determined. If it is not possible to obtain this from electrophysiological measurements (), or if there is not an accurate mathematical model which would allow numerical techniques to be used , , an approach such as that described in , which can estimate the phase response curve based on aggregate measurements such as the local field potential, could be used. Moreover, we expect that similar population-level control results would be found for other types of coupling, such as synaptic coupling and/or heterogeneous coupling, provided that the coupling strength is not too strong.
We imagine that the results from this paper can be useful to the neuroscience community in cases where there are biological or electronic hardware considerations which limit the allowed input magnitude for a stimulus. With deep brain stimulation becoming an increasingly adopted therapeutic technique for treatment of neurological disorders, this research extends ongoing research efforts to theoretically inform the optimal design of deep brain stimulation protocols.
Statements
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: (https://github.com/michaelzimet/MgOptChaoticDesync) (DOI: 10.5281/zenodo.15595745).
Author contributions
MZ: Software, Writing – review and editing, Methodology, Investigation, Writing – original draft, Conceptualization. FR: Conceptualization, Writing – review and editing, Software. JM: Writing – review and editing, Conceptualization, Methodology, Supervision, Software.
Funding
The author(s) declare that no financial support was received for the research and/or publication of this article.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declare that no Generative AI was used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnetp.2025.1646391/full#supplementary-material
References
1
AbouzeidA.ErmentroutB. (2009). Type-II phase resetting curve is optimal for stochastic synchrony. Phys. Rev. E80, 011911. 10.1103/PhysRevE.80.011911
2
AdamchicI.HauptmannC.BarnikolU. B.PawelczykN.PopovychO.BarnikolT. T.et al (2014). Coordinated reset neuromodulation for Parkinson’s disease: proof-of-concept study. Mov. Disord.29, 1679–1684. 10.1002/mds.25923
3
AslM. M.VahabieA.-H.ValizadehA.TassP. A. (2022). Spike-timing-dependent plasticity mediated by dopamine and its role in Parkinson’s disease pathophysiology. Front. Netw. Physiology2, 817524. 10.3389/fnetp.2022.817524
4
BrownE.MoehlisJ.HolmesP. (2004). On the phase reduction and response dynamics of neural oscillator populations. Neural Comput.16, 673–715. 10.1162/089976604322860668
5
CagnanH.PedrosaD.LittleS.PogosyanA.CheeranB.AzizT.et al (2017). Stimulating at the right time: phase-specific deep brain stimulation. Brain140, 132–145. 10.1093/brain/aww286
6
CagnanH.DenisonT.McIntyreC.BrownP. (2019). Emerging technologies for improved deep brain stimulation. Nat. Biotechnol.37, 1024–1033. 10.1038/s41587-019-0244-6
7
DanzlP.HespanhaJ.MoehlisJ. (2009). Event-based minimum-time control of oscillatory neuron models: phase randomization, maximal spike rate increase, and desynchronization. Biol. Cybern.101, 387–399. 10.1007/s00422-009-0344-3
8
ErmentroutB. (1996). Type I membranes, phase resetting curves, and synchrony. Neural Comput.8, 979–1001. 10.1162/neco.1996.8.5.979
9
ErmentroutB. (2002). Simulating, analyzing, and animating dynamical systems: a guide to XPPAUT for researchers and students. Philadelphia: Society for Industrial and Applied Mathematics.
10
FreyJ.CagleJ.JohnsonK. A.WongJ. K.HilliardJ. D.ButsonC. R.et al (2022). Past, present, and future of deep brain stimulation: hardware, software, imaging, physiology and novel approaches. Front. Neurol.13, 825178. 10.3389/fneur.2022.825178
11
GuckenheimerJ.HolmesP. J. (1983). Nonlinear oscillations, dynamical systems and bifurcations of vector fields. New York: Springer-Verlag.
12
HammondC.BergmanH.BrownP. (2007). Pathological synchronization in Parkinson’s disease: networks, models and treatments. Trends Neurosci.30, 357–364. 10.1016/j.tins.2007.05.004
13
HanselD.MatoG.MeunierC. (1995). Synchrony in excitatory neural networks. Neural Comp.7, 307–337. 10.1162/neco.1995.7.2.307
14
HodgkinA. L.HuxleyA. F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol.117, 500–544. 10.1113/jphysiol.1952.sp004764
15
HoltA. B.WilsonD.ShinnM.MoehlisJ.NetoffT. I. (2016). Phasic burst stimulation: a closed-loop approach to tuning deep brain stimulation parameters for Parkinson’s disease. PLoS Comput. Biol.12, e1005011. 10.1371/journal.pcbi.1005011
16
KeenerJ.SneydJ. (1998). Mathematical physiology. New York: Springer.
17
Khaledi-NasabA.KromerJ. A.TassP. A. (2022). Long-lasting desynchronization of plastic neuronal networks by double-random coordinated reset stimulation. Front. Netw. Physiology2, 864859. 10.3389/fnetp.2022.864859
18
KirkD. (1998). Optimal control theory. New York: Dover Publications.
19
KraussJ. K.LipsmanN.AzizT.BoutetA.BrownP.ChangJ. W.et al (2021). Technology of deep brain stimulation: current status and future directions. Nat. Rev. Neurol.17, 75–87. 10.1038/s41582-020-00426-z
20
KubotaS.RubinJ. E. (2018). Numerical optimization of coordinated reset stimulation for desynchronizing neuronal network dynamics. J. Comput. Neurosci.45, 45–58. 10.1007/s10827-018-0690-z
21
KuramotoY. (1984). Chemical oscillations, waves, and turbulence. Berlin: Springer.
22
LozanoA. M.LipsmanN. (2013). Probing and regulating dysfunctional circuits using deep brain stimulation. Neuron77, 406–424. 10.1016/j.neuron.2013.01.020
23
LückenL.YanchukS.PopovychO. V.TassP. A. (2013). Desynchronization boost by non-uniform coordinated reset stimulation in ensembles of pulse-coupled neurons. Front. Comput. Neurosci.7, 63. 10.3389/fncom.2013.00063
24
LysyanskyB.PopovychO. V.TassP. A. (2011). Desynchronizing anti-resonance effect of m: n on-off coordinated reset stimulation. J. Neural Eng.8, 036019. 10.1088/1741-2560/8/3/036019
25
ManosT.Diaz-PierS.TassP. A. (2021). Long-term desynchronization by coordinated reset stimulation in a neural network model with synaptic and structural plasticity. Front. Physiology12, 716556. 10.3389/fphys.2021.716556
26
MatchenT.MoehlisJ. (2018). Phase model-based neuron stabilization into arbitrary clusters. J. Comput. Neurosci.44, 363–378. 10.1007/s10827-018-0683-y
27
MatchenT.MoehlisJ. (2021). Leveraging deep learning to control neural oscillators. Biol. Cybern.115, 219–235. 10.1007/s00422-021-00874-w
28
MerrillD.BiksonM.JefferysJ. (2005). Electrical stimulation of excitable tissue: design of efficacious and safe protocols. J. Neurosci. Methods141, 171–198. 10.1016/j.jneumeth.2004.10.020
29
MoehlisJ. (2006). Canards for a reduction of the Hodgkin-Huxley equations. J. Math. Biol.52, 141–153. 10.1007/s00285-005-0347-1
30
MoehlisJ.ZimetM.RajabiF. (2025). Nearly optimal chaotic desynchronization of neural oscillators.
31
MongaB.MoehlisJ. (2019). Phase distribution control of a population of oscillators. Phys. D.398, 115–129. 10.1016/j.physd.2019.06.001
32
MongaB.FroylandG.MoehlisJ. (2018). “Synchronizing and desynchronizing neural populations through phase distribution control,” in 2018 Annual American Control Conference (ACC), 2808–2813. 10.23919/ACC.2018.8431114
33
MongaB.WilsonD.MatchenT.MoehlisJ. (2019). Phase reduction and phase-based optimal control for biological systems: a tutorial. Biol. Cybern.113, 11–46. 10.1007/s00422-018-0780-z
34
NabiA.MoehlisJ. (2009). “Charge-balanced optimal inputs for phase models of spiking neurons,” in Proceedings of the ASME 2009 Dynamic Systems and Control Conference, 685–687. 10.1115/DSCC2009-2541
35
NabiA.MirzadehM.GibouF.MoehlisJ. (2013). Minimum energy desynchronizing control for coupled neurons. J. Comput. Neurosci.34, 259–271. 10.1007/s10827-012-0419-3
36
NajeraR. A.MahavadiA. K.KhanA. U.BoddetiU.BeneV. A. D.WalkerH. C.et al (2023). Alternative patterns of deep brain stimulation in neurologic and neuropsychiatric disorders. Front. Neuroinform.17, 1156818. 10.3389/fninf.2023.1156818
37
NetoffT.SchwemmerM. A.LewisT. J. (2012). “Experimentally estimating phase response curves of neurons: theoretical and practical issues,” in Phase response curves in neuroscience. Editors SchultheissN. W.PrinzA. A.ButeraR. J. (New York, NY: Springer), 95–129.
38
NeumannW.-J.HornA.KühnA. A. (2007). Insights and opportunities for deep brain stimulation as a brain circuit intervention. Trends Neurosci.46, 472–487. 10.1016/j.tins.2023.03.009
39
OehrnC. R.CerneraS.HammerL. H.ShcherbakovaM.YaoJ.HahnA.et al (2024). Chronic adaptive deep brain stimulation versus conventional stimulation in Parkinson’s disease: a blinded randomized feasibility trial. Nat. Med.30, 3345–3356. 10.1038/s41591-024-03196-z
40
PopovychO. V.TassP. A. (2014). Control of abnormal synchronization in neurological disorders. Front. Neurology5, 268. 10.3389/fneur.2014.00268
41
PopovychO. V.HauptmannC.TassP. A. (2006). Control of neuronal synchrony by nonlinear delayed feedback. Biol. Cybern.95, 69–85. 10.1007/s00422-006-0066-8
42
PopovychO. V.LysyanskyB.RosenblumM.PikovskyA.TassP. A. (2017). Pulsatile desynchronizing delayed feedback for closed-loop deep brain stimulation. PLoS One12, e0173363. 10.1371/journal.pone.0173363
43
QinY.NobiliA. M.BassettD. S.PasqualettiF. (2023). Vibrational stabilization of cluster synchronization in oscillator networks. IEEE Open J. Control Syst.2, 439–453. 10.1109/OJCSYS.2023.3331195
44
RajabiF.GibouF.MoehlisJ. (2025). Optimal control for stochastic neural oscillators. Biol. Cybern.119, 9. 10.1007/s00422-025-01007-3
45
RosenblumM. G.PikovskyA. S. (2004a). Controlling synchronization in an ensemble of globally coupled oscillators. Phys. Rev. Lett.92, 114102. 10.1103/PhysRevLett.92.114102
46
RosenblumM. G.PikovskyA. S. (2004b). Delayed feedback control of collective synchrony: an approach to suppression of pathological brain rhythms. Phys. Rev. E70, 041904. 10.1103/PhysRevE.70.041904
47
RubchinskyL. L.ParkC.WorthR. M. (2012). Intermittent neural synchronization in Parkinson’s disease. Nonlinear Dyn.68, 329–346. 10.1007/s11071-011-0223-z
48
Sandoval-PistoriusS. S.HackerM. L.WatersA. C.WangJ.ProvenzaN. R.de HemptinneC.et al (2023). Advances in deep brain stimulation: from mechanisms to applications. J. Neurosci.43, 7575–7586. 10.1523/JNEUROSCI.1427-23.2023
49
SantanielloS.McCarthyM. M.MontgomeryJr., E. B.GaleJ. T.KopellN.SarmaS. V. (2015). Therapeutic mechanisms of high-frequency stimulation in Parkinson’s disease and neural restoration via loop-based reinforcement. Proc. Natl. Acad. Sci.112, E586–E595. 10.1073/pnas.1406549111
50
SpiliotisK.StarkeJ.FranzD.RichterA.KöhlingR. (2022). Deep brain stimulation for movement disorder treatment: exploring frequency-dependent efficacy in a computational network model. Biol. Cybern.116, 93–116. 10.1007/s00422-021-00909-2
51
TassP. A. (2003). A model of desynchronizing deep brain stimulation with a demand-controlled coordinated reset of neural subpopulations. Biol. Cybern.89, 81–88. 10.1007/s00422-003-0425-7
52
TassP. A. (2006). Phase resetting in medicine and biology: stochastic modelling and data analysis. Berlin: Springer.
53
VuM.SinghalB.ZengS.LiJ.-S. (2024). Data-driven control of oscillator networks with population-level measurement. Chaos34, 033138. 10.1063/5.0191851
54
WilsonD. (2020). Optimal open-loop desynchronization of neural oscillator populations. J. Math. Biol.81, 25–64. 10.1007/s00285-020-01501-1
55
WilsonD.MoehlisJ. (2014a). Locally optimal extracellular stimulation for chaotic desynchronization of neural populations. J. Comput. Neurosci.37, 243–257. 10.1007/s10827-014-0499-3
56
WilsonD.MoehlisJ. (2014b). Optimal chaotic desynchronization for neural populations. SIAM J. Appl. Dyn. Syst.13, 276–305. 10.1137/120901702
57
WilsonD.MoehlisJ. (2015a). Clustered desynchronization from high-frequency deep brain stimulation. PLoS Comput. Biol.11, e1004673. 10.1371/journal.pcbi.1004673
58
WilsonD.MoehlisJ. (2015b). Determining individual phase response curves from aggregate population data. Phys. Rev. E92, 022902. 10.1103/PhysRevE.92.022902
59
WilsonD.MoehlisJ. (2022). Recent advances in the analysis and control of large populations of neural oscillators. Annu. Rev. Control54, 327–351. 10.1016/j.arcontrol.2022.05.002
60
WilsonC. J.Beverlin IIB.NetoffT. (2011). Chaotic desynchronization as the therapeutic mechanism of deep brain stimulation. Front. Syst. Neurosci.5, 50. 10.3389/fnsys.2011.00050
61
WinfreeA. (1967). Biological rhythms and the behavior of populations of coupled oscillators. J. Theor. Biol.16, 15–42. 10.1016/0022-5193(67)90051-3
Summary
Keywords
optimal control, desynchronization, deep brain stimulation, Lyapunov exponent, network physiology
Citation
Zimet M, Rajabi F and Moehlis J (2025) Magnitude-constrained optimal chaotic desynchronization of neural populations. Front. Netw. Physiol. 5:1646391. doi: 10.3389/fnetp.2025.1646391
Received
13 June 2025
Accepted
29 September 2025
Published
21 October 2025
Volume
5 - 2025
Edited by
Eckehard Schöll, Technical University of Berlin, Germany
Reviewed by
Alexander Neiman, Ohio University, United States
Louis M Pecora, Naval Research Laboratory, United States
Updates
Copyright
© 2025 Zimet, Rajabi and Moehlis.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jeff Moehlis, moehlis@ucsb.edu
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.