Phase reduction and phase-based optimal control for biological systems: a tutorial

Abstract

A powerful technique for the analysis of nonlinear oscillators is the rigorous reduction to phase models, with a single variable describing the phase of the oscillation with respect to some reference state. An analog to phase reduction has recently been proposed for systems with a stable fixed point, and phase reduction for periodic orbits has recently been extended to take into account transverse directions and higher-order terms. This tutorial gives a unified treatment of such phase reduction techniques and illustrates their use through mathematical and biological examples. It also covers the use of phase reduction for designing control algorithms which optimally change properties of the system, such as the phase of the oscillation. The control techniques are illustrated for example neural and cardiac systems.

Appendices

Appendix A: Models

In this appendix, we give details of the mathematical models used in the main text.

Thalamic neuron model

The thalamic neuron model is given as

$$\begin{aligned} \dot{v}= & {} \frac{-I_L-I_\mathrm{Na}-I_K-I_T+I_\mathrm{b}}{C_m}+u(t),\\ \dot{h}= & {} \frac{h_{\infty }-h}{\tau _h},\\ \dot{r}= & {} \frac{r_{\infty }-r}{\tau _r}, \end{aligned}$$

where

$$\begin{aligned} h_\infty= & {} 1/(1+\exp ((v+41)/4)),\\ r_\infty= & {} 1/(1+\exp ((v+84)/4)),\\ \alpha _h= & {} 0.128\exp (-(v+46)/18),\\ \beta _h= & {} 4/(1+\exp (-(v+23)/5)),\\ \tau _h= & {} 1/(\alpha _h+\beta _h),\\ \tau _r= & {} (28+\exp (-(v+25)/10.5)),\\ m_\infty= & {} 1/(1+\exp (-(v+37)/7)),\\ p_\infty= & {} 1/(1+\exp (-(v+60)/6.2)),\\ I_L= & {} g_L(v-e_L),\\ I_\mathrm{Na}= & {} g_\mathrm{Na}({m_\infty }^3)h(v-e_\mathrm{Na}),\\ I_K= & {} g_K((0.75(1-h))^4)(v-e_K),\\ I_T= & {} g_T(p_\infty ^2)r(v-e_T),\\ C_m= & {} 1, \qquad g_L = 0.05, \qquad e_L = -70, \\ g_\mathrm{Na}= & {} 3, \qquad e_\mathrm{Na} = 50,\\ g_K= & {} 5,\qquad e_K = -90,\qquad g_T = 5,\\ e_T= & {} 0, \qquad I_\mathrm{b} = 5. \end{aligned}$$

Morris–Lecar model

The Morris–Lecar model is given as

$$\begin{aligned} C_M\dot{v}= & {} I_\mathrm{b}-g_L(v-E_L)-g_Kn(v-E_K)\\&-g_\mathrm{Ca}m_\infty (v)(v-E_\mathrm{Ca}),\\ \dot{n}= & {} \phi (n_\infty (v)-n)/\tau _n(v),\\ m_\infty (v)= & {} 0.5\left( 1+\tanh \left( \frac{v-v_1}{v_2}\right) \right) ,\\ \tau _n(v)= & {} \frac{1}{\left( \cosh \left( \frac{v-v_3}{2v_4}\right) \right) },\\ n_\infty (v)= & {} 0.5\left( 1+\tanh \left( \frac{v-v_3}{v_4}\right) \right) ,\\ \phi= & {} 0.067 , \quad g_\mathrm{Ca}=4,\quad g_K=8, \quad g_L=2, \\ E_\mathrm{Ca}= & {} 120, \quad E_k=-84, \\ E_l= & {} -60, \quad v_1=-1.2, \quad v_2=18, \quad v_3=12, \\ v_4= & {} 17.4, \quad C_M=20. \end{aligned}$$

Fox–McHarg–Gilmour (FMG) model

The FMG model [26] describes the electrophysiological behavior of a canine ventricular myocytes with behavior governed by various potassium, sodium, and calcium currents. The transmembrane voltage dynamics are governed by the flow of ionic current across the cell membrane

$$\begin{aligned} \dot{V} = -I_{\mathrm{ion}} - I_{\mathrm{stim}}. \end{aligned}$$

(122)

Here $I_{\mathrm{stim}}$ represents a stimulus current used to elicit action potentials, and $I_{\mathrm{ion}}$ is the total membrane current density, both of which have units of $\upmu \hbox {A}/\upmu \hbox {F}$. $I_{\mathrm{ion}}$ is comprised of 13 individual currents, which are determined by 13 total state variables. A full description of the equations and nominal parameters is given in [26].

Appendix B: Derivation of the Euler–Lagrange equations

The Euler–Lagrange equations can be derived using the methods of analytical mechanics [30] or optimal control theory [47]. Suppose we want to find the function q(t) which extremizes the functional

$$\begin{aligned} C[q(t)] = \int _0^T {{\mathcal {L}}}(q,\dot{q}) \hbox {d}t. \end{aligned}$$

The functional derivative has the property that for any function v(t),

$$\begin{aligned} \frac{\delta C}{\delta q} \cdot v= & {} \lim _{\epsilon \rightarrow 0} \frac{C[q+\epsilon v] - C[q]}{\epsilon } \\= & {} \lim _{\epsilon \rightarrow 0} \frac{1}{\epsilon } \int _0^T \left\{ {{\mathcal {L}}}(q+\epsilon v, \dot{q} + \epsilon \dot{v}) - {{\mathcal {L}}}(q,\dot{q}) \right\} \hbox {d}t \\= & {} \lim _{\epsilon \rightarrow 0} \frac{1}{\epsilon } \int _0^T \left\{ {{\mathcal {L}}}(q,\dot{q}) + \frac{\partial {{\mathcal {L}}}}{\partial q} \epsilon v + \frac{\partial {{\mathcal {L}}}}{\partial \dot{q}} \epsilon \dot{v} \right. \\&\left. \quad + \cdots - {{\mathcal {L}}}(q,\dot{q}) \right\} \hbox {d}t \\= & {} \int _0^T \left( \frac{\partial {{\mathcal {L}}}}{\partial q} v + \frac{\partial L}{\partial \dot{q}} \dot{v} \right) \hbox {d}t \\= & {} \int _0^T \left[ \frac{\partial {{\mathcal {L}}}}{\partial q} - \frac{d}{dt} \left( \frac{\partial {{\mathcal {L}}}}{\partial \dot{q}} \right) \right] v dt + \left[ \frac{\partial {{\mathcal {L}}}}{\partial \dot{q}} v \right] _0^\mathrm{T}, \end{aligned}$$

where the last equation follows from integration by parts applied to the second term of the previous line. If we suppose that $v(0) = v(T) = 0$, then

$$\begin{aligned} \frac{\delta C}{\delta q} \cdot v = \int _0^T \left[ \frac{\partial {{\mathcal {L}}}}{\partial q} - \frac{d}{dt} \left( \frac{\partial {{\mathcal {L}}}}{\partial \dot{q}} \right) \right] v \hbox {d}t. \end{aligned}$$

For this to be an extremum, it must hold for any v. Therefore, we obtain the Euler–Lagrange equation

$$\begin{aligned} \frac{\hbox {d}}{\hbox {d}t} \left( \frac{\partial {{\mathcal {L}}}}{\partial \dot{q}} \right) = \frac{\partial {{\mathcal {L}}}}{\partial q}. \end{aligned}$$

(123)

To incorporate constraints into this approach, it is instructive to first consider the simpler optimization problem where one wants to find an extremum of the function $\mathbf{f}(\mathbf{x})$ subject to the constraint that $\mathbf{g}(\mathbf{x}) = \mathbf{c}$. Recognizing that at an extremum the level surface of $\mathbf{f}(\mathbf{x})$ must be tangent to the surface defined by $\mathbf{g}(\mathbf{x}) = \mathbf{c}$, we see that

$$\begin{aligned} \nabla \mathbf{f} = \lambda \nabla \mathbf{g} \end{aligned}$$

(124)

for some scalar $\lambda$, which is called the Lagrange multiplier. Finding the extremum of $\mathbf{f}$ while simultaneously satisfying (124) is equivalent to finding the extremum of

$$\begin{aligned} \mathbf{F}(\mathbf{x},\lambda ) \equiv \mathbf{f}(\mathbf{x}) - \lambda [\mathbf{g}(\mathbf{x}) - \mathbf{c}]. \end{aligned}$$

Indeed,

$$\begin{aligned} \nabla \mathbf{F} = 0 \quad\Rightarrow & {} \quad \nabla \mathbf{f} = \lambda \nabla \mathbf{g}, \\ \frac{\partial \mathbf{F}}{\partial \lambda } = 0 \quad\Rightarrow & {} \quad \mathbf{g}(\mathbf{x}) = \mathbf{c}. \end{aligned}$$

In optimal control problems, it is necessary to minimize or maximize the cost function subject to the constraint that the dynamics must satisfy the appropriate evolution equation. For example, consider the energy-optimal phase control problem of designing the input u(t) such that the cost function ${{\mathcal {G}}}[u(t)]$ given by (28) is minimized, subject to the constraint that the solution must satisfy (27). We can rewrite (27) as

$$\begin{aligned} \frac{\hbox {d} \theta }{\hbox {d}t} - \omega - Z(\theta ) u(t) = 0. \end{aligned}$$

Integrating this, we obtain

$$\begin{aligned} {{\mathcal {K}}}[u(t)] \equiv \int _0^{T_1} \left[ \frac{\hbox {d} \theta }{\hbox {d}t} - \omega - Z(\theta ) u(t) \right] \hbox {d}t = 0. \end{aligned}$$

We therefore want to find u(t) such that the following hold:

$$\begin{aligned} \frac{\delta {{\mathcal {G}}}[u(t)]}{\delta u} = 0, \quad {{\mathcal {K}}}[u(t)] = 0. \end{aligned}$$

By analogy with the above constrained optimization example, the level surfaces of ${{\mathcal {G}}}[u(t)]$ will be tangent to surfaces for which ${{\mathcal {K}}}[u(t)] = 0$. Thus,

$$\begin{aligned} \frac{\delta {{\mathcal {G}}}[u(t)]}{\delta u} = -\lambda (t) \frac{\delta {{\mathcal {K}}}[u(t)]}{\delta u} \end{aligned}$$

(125)

for some scalar function $\lambda (t)$. [Without loss of generality, we have inserted a minus sign on the right-hand side of (125).]. Rearranging (125) gives

$$\begin{aligned} \frac{\delta }{\delta u} \left\{ {{\mathcal {G}}}[u(t)] + \lambda {{\mathcal {K}}}[u(t)] \right\} = 0, \end{aligned}$$

which is identical to the cost function (29). We now use (123) with q(t) taken in turn to be u(t), $\lambda (t)$, and $\theta (t)$, which gives Eqs. (30–32) in the main text.

Appendix C: Solving a two-point boundary value problem

Consider a general two-point boundary value problem

$$\begin{aligned} \dot{y} = f(t,y), \quad y \in {\mathbb {R}}^n, \quad 0\le t \le b, \end{aligned}$$

(126)

with the linear boundary condition

$$\begin{aligned} B_0y(0)+B_by(b)=a, \quad B_0,\ B_b \in {\mathbb {R}}^{n\times n}. \end{aligned}$$

To solve such a boundary value problem, we integrate Eq. (126) with the initial guess $c=y(0)$ and calculate the function g(c):

$$\begin{aligned} g(c)=B_0c+B_by(b)-a, \end{aligned}$$

where y(b) is the solution at time b with the initial condition c. If we had chosen the correct initial condition c, g(c) would be 0. Based on the current guess $c^\nu$, and the $g(c^\nu )$ value, we choose the next initial condition by the Newton Iteration as

$$\begin{aligned} c^{\nu +1}=c^\nu -\left( \left. \frac{\partial g}{\partial c}\right| _{c^\nu }\right) ^{-1}g(c^\nu ). \end{aligned}$$

(127)

We compute the Jacobian $J=\left. \frac{\partial g}{\partial c}\right| _{c^\nu }$ numerically as

$$\begin{aligned} J_i=\frac{g^+-g^-}{2\epsilon }, \end{aligned}$$

where

$$\begin{aligned} g^+= & {} g\left( c^\nu +e_i \epsilon \right) ,\\ g^-= & {} g\left( c^\nu -e_i \epsilon \right) , \end{aligned}$$

$J_i$ is the $i\mathrm{th}$ column of J, $\epsilon$ is a small number, and $e_i$ is a column vector with 1 in the $i\mathrm{th}$ position and 0 elsewhere.