Isochron

Post-publication activity

Curator: Kresimir Josic Figure 1: This animation shows the evolution of two points lying on the same isochron. The black point lies on the periodic orbit, and the red point approaches it asymptotically, resulting in oscillations having the same phase.

An isochron (Greek for equal time) of a dynamical system is a set of initial conditions resulting in oscillations having the same phase.

Definition

Consider an autonomous system of ordinary differential equations

$\frac{d x}{dt} = {F}(x), \qquad x \in R^n, \qquad (n \ge 2)$

having a hyperbolic periodic orbit $$\gamma(t)$$ of period $$T\ .$$ For each point $$x$$ in the basin of attraction of the periodic orbit there exists a unique $$\theta(x)$$ such that

$\lim_{t \rightarrow \infty} |x(t) - \gamma(t+\theta(x))| = 0,$

where $$x(t)$$ is a trajectory (solution) of the dynamical system above starting with the initial point $$x\ .$$ The value $$\theta(x)$$ is called the asymptotic (or latent) phase of $$x\ .$$ This notion allows us to assign a phase to each point in the basin of attraction of a periodic orbit, so that $$\theta(x)$$ is a function on this domain, taking values in $$[0, T)\ .$$ (The phase is often normalized by $$T/(2\pi)$$ so that the function takes values in $$[0,2\pi)$$).

The collection of all points in the basin of attraction of $$\gamma$$ with the same asymptotic phase is called an isochron; in other words, an isochron is a level set of the function $$\theta(x)\ .$$ For an exponentially stable periodic orbit of an autonomous system of $$n$$ ordinary differential equations, each isochron is an $$(n-1)$$-dimensional hypersurface. The flow of two initial conditions on the same isochron is shown in the figure above.

More intuitively, consider the forward orbit, $$x(t)\ ,$$ of an initial condition $$x$$ and the collection of points $$\{x(kT)\}_{k=0}^{\infty}\ .$$ These points are obtained by observing the orbit only at times that are an integer multiple of the period of $$\gamma\ ,$$ so that they are defined by a Poincaré map. All points in this sequence lie on an isochron (Winfree 2001). Therefore, the isochron is a particular Poincaré section for which the time of the first return for any point equals $$T\ .$$ Points at which isochrons of a periodic orbit cannot be defined form the phaseless set.

An excellent discussion of isochrons with historical references can be found in Winfree (2001).

Existence and properties

Isochrons can be shown to exist for any stable hyperbolic periodic orbit. In this case the isochrons are codimension one manifolds as smooth as the vector field, and transversal to the periodic orbit $$\gamma\ .$$ Their union covers an open neighborhood of $$\gamma\ .$$ This can be proved directly by using the Implicit Function Theorem (Guckenheimer 1975, Coddington and Levinson 1955), and is also a consequence of results on normally hyperbolic invariant manifolds (Wiggins 1994).

An example

This example taken from Winfree (2001) illustrates some properties of isochrons. Consider the planar differential equation in polar coordinates

$\frac{dR}{dt} = (1-R)R^2$

$\frac{d\phi}{dt} = R.$

This system has an attracting periodic orbit $$\gamma = \{R = 1\}\ .$$ As above, the function $$\theta(\phi,R)$$ defines the asymptotic phase, and isochrons are the level sets of this function -- that is, 1-dimensional sets of points with the same asymptotic phase. Since the vector field is symmetric under rotations around the origin, the family of isochrons must be invariant under that symmetry. Consistent with this symmetry, now assume that the function $$\theta(\phi,R)$$ has the form $$\theta(\phi,R)=\phi - f(R)\ ,$$ for some function $$f$$ to be determined. Then each isochron is a graph in polar coordinates. In particular, the isochron corresponding to asymptotic phase $$c$$ is the set of points $$\{\phi,R\}$$ satisfying $$c=\phi - f(R)\ .$$

On the limit cycle $$\frac{d\phi}{dt} = 1\ .$$ Moreover, all points on an isochron share the same phase value, and therefore

$\frac{d\theta(x)}{dt} = 1\ .$

It follows that

$\frac{d\theta}{dt} = 1 = \frac{d\phi}{dt} - \frac {df}{dR} \frac{dR}{dt} \qquad \Longrightarrow \qquad \frac {df}{dR}= \frac{\frac{d\phi}{dt} - 1}{\frac{dR}{dt} }\qquad \Longrightarrow \qquad f(R) = \frac{1}{R} + const.$

Setting the constant to -1, we find that each isochron is therefore defined by $$\phi = c + \frac{1}{R} - 1\ .$$ The constant is set to -1 so that the value of the asymptotic phase agrees with the angular coordinate on the limit cycle. Figure 2: Four representative isochrons of the periodic orbit (red) considered in the example.

In this example isochrons are defined everywhere in the plane, except at the origin, which is the only element of the phaseless set. Note that any neighborhood of the phaseless set in this example intersects all isochrons, a situation which can be expected generally (Guckenheimer 1975).

If we define a new phase coordinate $$\psi = \phi - f(R)$$ away from the origin, which coincides precisely with the asymptotic phase, then the original differential equation takes the simpler form

$\frac{dR}{dt} = (1-R)R^2$ $\frac{d\psi}{dt} = 1.$

This illustrates how the asymptotic phase can be used to find a change of coordinates in which the dynamics of the phase coordinate is decoupled from the remaining coordinates.

Applications

Isochrons have found a number of applications. As noted in the previous example, they can be used to effectively reduce the dimension of the equation in the neighborhood of a periodic orbit. They are also useful in:

• Extending the notion of phase of a periodic orbit to a neighborhood of the periodic orbit. It is conventional to take the phase as advancing linearly in time, with values in the interval $$[0, T)\ ,$$ $$[0,2 \pi)\ ,$$ as in the discussion above, or in $$[0,1)\ .$$
• Obtaining an expression governing the behavior of the phase difference between weakly coupled oscillators (in particular, see the Kuramoto approach in that entry).
• Giving a sense of how long a trajectory spends in different regions of phase space: for example, a trajectory moves slowly through regions of phase space where isochrons spaced equally in phase lie close to each other.

Computing Isochrons

Except for simple examples such as that given above, isochrons cannot be calculated analytically. There are two main numerical methods for finding isochrons.

Numerical methods

1. Using forward integration. One chooses an initial condition $$x_0$$ on the periodic orbit, and integrates forward for some long time $$\tau\ .$$ Then one takes an initial condition $${ x_{test}}$$ in the basin of attraction of the periodic orbit, and integrates forward for the same time $$\tau\ .$$ Let $$\Phi_t(x)$$ denote the flow of the dynamical system, i.e., the trajectory starting with $$x\ .$$ If $$|\Phi_\tau(x_0) - \Phi_\tau(x_{test})| < \epsilon$$ for some small tolerance $$\epsilon\ ,$$ then $$x_{test}$$ approximately lies on the isochron defined by $$x_0\ .$$ If not, one considers a different $$x_{test}\ ,$$ perhaps from an algorithm which predicts a better choice.
2. Using backward integration. One chooses an initial condition $$x_0$$ on the periodic orbit, and integrates backward for some time $$\tau\ .$$ Then one takes an initial condition $$x_{test}$$ which is close to $$x_0\ ,$$ that is with $$|x_0 - x_{test}| < \epsilon$$ for some small tolerance $$\epsilon\ ,$$ and integrates backward for the same time $$\tau\ .$$ The points $$\Phi_{-\tau}( x_0)$$ and $$\Phi_{-\tau}( x_{test})$$ approximately lie on the same isochron.

The second method is typically much more computationally efficient, although one must be careful about the instability associated with integrating backward in time. Ten isochrons equally spaced in phase calculated using this method for a realistic planar Type I neuron model (a Hindmarsh-Rose neuron with equations as given in Appendix C of Brown et al. (2004) and a baseline current of $$I^b = 7.6 \mu A/cm^2$$) are shown as dashed lines in the figure on the right. Figure 3: Example of isochrons equally spaced in phase for a planar Type I neuron model.

Twenty isochrons equally spaced in phase calculated using this method for a realistic planar Type II neuron model (Hodgkin-Huxley Model with standard parameters and baseline current $$I^b = 10 \mu A/cm^2\ ,$$ with the reduction $$m \rightarrow m_\infty$$ and $$h = 0.8 - n\ ;$$ see Keener and Sneyd (1998)) are shown as dashed lines in the second figure on the right. Figure 4: Example of isochrons equally spaced in phase for a planar Type II neuron model.

Local methods

Local approximations to isochrons in the neighborhood of a base point $$x^\gamma$$ on the limit cycle $$\gamma$$ can be developed using an alternative approach; while this approach typically also requires numerical solution, in some cases the local approximations (which depend only on the vector field linearized around $$\gamma$$) can be obtained analytically even though the global" isochrons cannot. The idea is to calculate the gradient $$\nabla \theta \left( x^\gamma \right)$$ of the asymptotic phase at the base point; the isochron, being a set of constant asymptotic phase, is then tangent to the $$n-1$$ dimensional plane normal to this gradient. We now describe two methods for computing the gradient $$\nabla \theta\left(x^\gamma \right)\ .$$

Direct method

By definition

$\frac{\partial \theta}{\partial x_i} \left( x^\gamma \right) = \lim_{\Delta x_i \rightarrow 0} \frac{ \Delta \theta}{\Delta x_i}$

where $$\Delta \theta = \left[ \theta(x^\gamma + \Delta x_i \hat{i}) - \theta(x^\gamma) \right]$$ is the change in $$\theta(x)$$ resulting from a perturbation $$x_i \rightarrow x_i + \Delta x_i$$ from the base point $$x^\gamma$$ on $$\gamma$$ in the direction of the $$i^{th}$$ coordinate. Since $$\dot{\theta}=2 \pi /T$$ everywhere in the neighborhood of $$\gamma\ ,$$ the difference $$\Delta \theta$$ is preserved under the flow; thus, it may be measured in the limit as $$t \rightarrow \infty\ ,$$ when the perturbed trajectory has collapsed back to the limit cycle $$\gamma\ .$$ That is, $$\frac{\partial \theta}{\partial x_i}\left( x^\gamma \right)$$ can be found by comparing the phases of solutions in the infinite-time limit starting on and infinitesimally shifted from base points on $$\gamma$$ (Winfree 2001, Glass and Mackey 1998, Brown et al. 2004).

The figure on the right illustrates the direct method of computing $$\frac{\partial \theta}{\partial x_i}\left( x^\gamma \right)$$ at the base point $$x^\gamma$$ by taking the limit of $$\Delta \theta / \Delta x_i$$ for vanishingly small perturbations $$\Delta x_i\ .$$ One can calculate $$\Delta \theta$$ in the limit $$t \rightarrow \infty$$ numerically, or, in some cases, analytically (Brown et al. 2004). Figure 5: The direct method for computing $$\frac{\partial \theta}{\partial x_i}\left( x^\gamma \right)$$ at the base point $$x^\gamma$$ is to take the limit of $$\Delta \theta / \Delta x_i$$ for vanishingly small perturbations $$\Delta x_i\ .$$ One can calculate $$\Delta \theta$$ in the limit $$t \rightarrow \infty$$ numerically, or, in some cases, analytically (Brown et al. 2004).

Another technique for finding $$\nabla \theta \left( x^\gamma \right)$$ involves solving the associated adjoint problem (Ermentrout and Kopell 1991, Hoppensteadt and Izhikevich 1997, Ermentrout 2002); this procedure is automated in the program XPP (Ermentrout 2002) and is equivalent to the direct method discussed above. We note that this equivalence is implicit in the calculation of coupling functions presented in Hoppensteadt and Izhikevich (1997) and Ermentrout (2002); see also the appendix in Brown et al. (2004).

In this approach, one solves for $$\nabla \theta \left( x^\gamma (t) \right)\ ,$$ where the $$t$$ dependence indicates parametrization of the limit cycle by time. The underlying equation is (Ermentrout and Kopell 1991, Hoppensteadt and Izhikevich 1997, Ermentrout 2002):

$\frac{d \nabla \theta \left( x^\gamma (t) \right) }{d t} = - DF^T( x^\gamma (t)) \, \nabla \theta \left( x^\gamma (t) \right)$

Here the matrix $$D F^T( x^\gamma (t))$$ is the transpose (i.e., adjoint) of the (real) matrix $$D F( x^\gamma (t))\ .$$ This equation is subject to the condition

$\nabla \theta \left( x^\gamma (0) \right) \cdot F(x^\gamma(0)) = 1.$

Since $$\nabla \theta \left( x^\gamma (t) \right)$$ evolves in $$R^n\ ,$$ this condition supplies only one of $$n$$ required initial conditions. The rest arise from requiring that the solution be $$T$$-periodic (Hoppensteadt and Izhikevich 1997, Ermentrout 2002, Ermentrout and Kopell 1991).

Note that the adjoint equation and the following condition correspond to equations (9.16) and (9.17) of Hoppensteadt and Izhikevich (1997), with the identification of $$\nabla \theta = Q\ .$$ This is the adjoint problem that XPP solves to numerically find the PRC $$Q_{\rm XPP}\ .$$

Generalizations

Isochrons are the leaves of the invariant foliation of the stable manifold of a periodic orbit $$\gamma\ .$$ In this discussion it was assumed that $$\gamma$$ is stable so that the stable manifold contains some open neighborhood of $$\gamma\ .$$ As long as $$\gamma$$ is hyperbolic the stable manifold of $$\gamma$$ can be foliated in a similar way, so that isochrons can be defined even when $$\gamma$$ is unstable. If $$\gamma$$ is not stable, its stable manifold will not contain a neighborhood of $$\gamma\ ,$$ so that isochrons in this case fail to foliate an open neighborhood of $$\gamma\ .$$

Isochrons can also be defined in the case of excitable systems (Rabinovitch et al. 1994, Ichinose et al. 1998, Coombes and Osbaldestin 2000).

Deterministic chaos and stochastic dynamics can show a rhythmic component, and such types of irregular oscillations have been described by phase models. Without a periodic orbit, the standard definition of phase via isochrones is however not applicable. To address this problem, generalizations of isochrones to the domain of irregular oscillations have been proposed for chaotic oscillations (Beck and Josić 2003, Schwabedal et al. 2012), and stochastic oscillations (Schwabedal and Pikovsky 2013). Here, isochrones are constructed as sections in state space, showing constant first return times: for chaotic oscillations the condition of return times is fulfilled approximately, for stochastic oscillations in an average sense. The approach does not make reference to an unstable periodic orbit and thus allows for applications of isochrones to a wide range of oscillatory systems.