Computing isochrons
Ioschron geometry can be very complicated, even for planar models. Biological models that motivate the use of isochrons, often fall into the class of slow-fast systems. Slow-fast systems are characterized as having a variables with at least two different timescales, making some variables fast than others. In the case of biological models, changes in the cellular currents are typically much slower than of the membrane voltage. However, large timescale separations can lead to sensitivity in the models, making the isochrons difficult to resolve with conventional integration-based methods. Notably, Arthur Winfree, who was a theoretical biologist, computed the isochrons of the FitzHugh-Nagumo model, a model that describes an ionic model. However, he encountered the phase sensitivity, which lead to breaks in his isochrons, which meant he could not compute the isochrons as single connected curves. We set out to complete these isochrons.
Our approach to computing isochrons uses numerical continuation.
We first compute the linear approximation of the isochron.
The linear approximation E(gamma_0) is then extended as the set of initial points of a family of orbit segments that reach the linear approximation after time equal to one period of the periodic orbit. Each orbit segment in this family reaches the linear approximation at a different distance eta from gamma_0, hence, this distance parameterizes the family. We construct a two-point boundary-value problem where this parametrised family is a solution. The family can then be computed with numerical continuation. The limit cycle Gamma provides a first solution with zero distance from gama_0. This distance parameter is then increased, computing the corresponding orbit segments, where the initial points of the orbit segments trace out a one-dimensional extension of the isochron.
This has the advantage that each isochron was computed directly as a single parameterized curve, where the step size could be adapted automatically when encountering areas of sensitivity.
However, the method required knowledge of the phase point gamma_theta at to compute isochrons at phases other than theta=0. Then the isochron could be computed as described above.
I adapted this method such that all the isochrons could be computed using a single linear approximation at gamma_0. This still required that the end point of the family of orbit segments reached the linear approximation, but by increasing the time of the orbit segments, the initial points of the family, i.e. the isochron extension, were computed at a different phase. Specifically, the isochron at phase theta, required a time of T=(2-theta)T_Gamma. This meant that the initial solution, the periodic orbit, needed extending backwards in time. Since integration in backwards time is numerically delicate, our approach was to first shift the periodic orbit solution such that the base point was at gamma_theta, and extend the periodic orbit in forward time.
In this way the expensive linear approximation step of the method was performed only once, from which all isochrons could be computed.
With this adaptation, the isochron geometry could be efficiently computed and we could complete Winfree's puzzle. We were able to show that the isochron geometry exhibited hooked turns, where portions of isochron accumulated in small regions of each other, leading to areas of high phase sensitivity. We showed how these turns were directly related to the slow manifold in the area - a feature typical of slow-fast systems. This connection and the adaptation to the numerical method formed the basis of the article Solving Winfree's Puzzle: