By Jeff Moehlis and Dan Wilson
The traditional boundaries between engineering and the life sciences are rapidly blurring as interdisciplinary researchers explore fundamental and practical questions in biology and medicine. A particularly fruitful avenue for such explorations uses classical phase reduction techniques and a related technique we recently proposed, called isostable reduction, to rigorously derive simplified equations. Optimal control theory can then be applied to these simplified equations to find novel potential treatments for medical ailments such as Parkinson’s disease and cardiac arrhythmias.
Many physical, biological, and technological systems produce rhythmic oscillations. A powerful classical technique for the analysis of such oscillators is the rigorous reduction to phase models, with a single variable describing the phase of the oscillation with respect to some reference state. Through reduction to phase models, one can understand the dynamics of high dimensional and analytically intractible models in a more convenient form (see, e.g., [9, 2, 4, 1, 3]),
\[\frac{d \theta}{dt} = \omega + Z(\theta) u(t), \:\:\:\:\: (1)\]
where \(\theta\) is the phase variable for the oscillator, \(\omega\) is the natural frequency in the absence of external forcing, and \(Z (\theta)\) is the phase response curve (PRC), which can in principle and—often in practice—be measured experimentally. Such a reduction is based on isochrons for the system, which are foliations of phase space that extend the notion of the phase of a stable periodic orbit to the basin of attraction of the periodic orbit; see the top panel of Figure 1, and a presentation with videos illustrating isochrons, isostables, and other aspects of this article. Each point in the basin of attraction lies on only one isochron, and two points on the same isochron converge to the periodic orbit with the same phase.
Figure 1. Left: Black lines show isochrons for a dynamical system with a stable periodic orbit (shown as a gray closed curve). Right: Isostables for a system with a stable fixed point (shown as a gray dot). Trajectories with different initial conditions on the same isochron/isostable are shown as dashed red lines, with white dots showing the solutions at particular times. Adapted from [7].
While an equation of the form \((1)\) is relatively simple to work with, it retains the essence of the system under study, allowing for progress to be made on difficult problems. For example, consider the hypothesis that pathological synchronization of spiking neurons in the basal ganglia-cortical loop within the brain is a factor contributing to tremors exhibited by patients with Parkinson’s disease, along with the established treatment option called deep brain stimulation in which a neurosurgeon implants an electrode that can inject a current into the brain tissue. These suggest that it might be useful to design a single electrical stimulus, which desynchronizes a population of neural oscillators. Using only the dynamic equations describing each neuron in terms of its full nonlinear dynamics, it would be difficult to make any analytical progress, but with the phase-reduced models for the neurons, one can show that exponential separation will occur with a stimulus \(u(t)\) which yields a positive Lyapunov exponent [7]
\[\Lambda(\tau) = \frac{1}{\tau} \int_0^\tau Z'(\theta(s)) u(s) ds. \:\:\:\:\: (2)\]
The Lyapunov exponent provides a control objective for which optimal control theory can be used: find the stimulus which maximizes the Lyapunov exponent while simultaneously minimizing the power used by the stimulus [7]. The latter goal is desirable because it allows the battery which generates the electrical stimulus to have a longer life. Numerical results from this approach are shown in Figure 2 for a population of thalamic neurons. In [7], we also showed that this methodology is robust for weak neuron coupling and heterogeneities.
Figure 2. Desynchronization of a neural population. (Top) Black lines show voltage traces for 100 simulated coupled neurons, and red line shows average voltage. When the average voltage exceeds the threshold defined by the horizontal purple line, the optimal control stimulus is applied (bottom) until peaks in voltage no longer exceed the blue line. Adapted from [7].
For systems which have a stable fixed point rather than a stable periodic orbit, one can define isostables [5], which are sets of points in phase space that approach the fixed point together and are analogous to isochrons for asymptotically periodic systems; see the bottom panel of Figure 1. In [8] we show that one can perform an isostable reduction for such systems, which leads to an equation of the form
\[\frac{d \psi}{dt} = \kappa + I(\psi) u(t),\:\:\:\:\: (3)\]
where \(\psi\) is the scalar isostable coordinate similar in nature to the phase in a periodic system, \(\kappa\) is a constant which describes the rate of change of the isostable coordinate in the absence of external forcing, \(I(\psi)\) is the isostable response curve, which is analogous to the PRC, and \(u(t)\) is an external stimulus. We demonstrate the utility of isostable reduction by considering the problem of finding an electrical stimulus, \(u(t)\), which can eliminate the cardiac arrhythmia known as alternans, the beat-to-beat alternation of electrochemical cardiac dynamics at a constant rate of pacing, which has been implicated as a possible precursor to more serious cardiac arrhythmias. In [8], we applied optimal control theory to \((3)\) to find a stimulus which eliminates alternans by stabilizing a periodic heartbeat while simultaneously minimizing the power associated with the stimulus. Isostable reduction allows such an approach to be used even for high-dimensional models of cardiac activity, as shown in Figure 3.
Figure 3. Control of alternans. (Top) Voltage trace showing action potentials for a 13-dimensional model for cardiac cell dynamics subjected to an applied control stimulus (middle) found by applying optimal control theory to the equation obtained through isostable reduction. (Bottom) Successive values of the action potential duration (APD) when the control is off (not shaded) and on (shaded). Adapted from [8].
Perhaps the greatest strength of both reduction strategies is their experimental applicability to living systems. In an experimental setting, the full right hand side of the dynamical equations are usually unknown, but phase response curves can still be measured using direct methods [3, 6]. We propose that with a similar experimental protocol, isostable response curves could be measured in living tissue as well. As technology continues to evolve at a rapid pace, scientists are able to observe and record more and more of the dynamical behavior associated with debilitating diseases. Both phase reduction for systems with a stable periodic orbit and isostable reduction for systems with a stable fixed point offer tremendous promise for making sense of this dynamical behavior and for pointing to new treatment options.
This article is based on a talk given by Professor Moehlis at the 2015 SIAM Conference on Applications of Dynamical Systems.
Acknowledgments: This research was supported by National Science Foundation grants NSF-1264535 and NSF-1363243.
References
[1] Brown, E., Moehlis, J., & Holmes, P. (2004). On the phase reduction and response dynamics of neural oscillator populations. Neural Comp., 16, 673–715.
[2] Guckenheimer, J. (1975). Isochrons and phaseless sets. J. Math. Biol., 1, 259–273.
[3] Izhikevich, E. M. (2007). Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. London: MIT Press.
[4] Kuramoto, Y. (1984). Chemical Oscillations, Waves, and Turbulence. Berlin: Springer.
[5] Mauroy, A., Mezic, I., & Moehlis, J. (2013). Isostables, isochrons, and Koopman spectrum for the action-angle representation of stable fixed point dynamics. Physica D, 261, 19-30.
[6] Netoff, T., Schwemmer, M. A., & Lewis T. J. (2012). Experimentally estimating phase response curves of neurons: theoretical and practical issues in Phase Response Curves in Neuroscience, (pp 95–129). New York, NY: Springer.
[7] Wilson, D. & Moehlis, J. (2014).Optimal chaotic desynchronization for neural populations. SIAM J. Appl. Dyn. Syst., 13, 276-305.
Jeff Moehlis is a professor and Dan Wilson is a graduate student in the Department of Mechanical Engineering at the University of California, Santa Barbara.This article is based on a talk given by Professor Moehlis at the 2015 SIAM Conference on Applications of Dynamical Systems.