Phase response curves and the role of coordinates

The “infinitesimal phase response curve” (PRC) is a common tool used to analyze phase resetting in the natural sciences in general and neuroscience in particular. We make the observation that the PRC with respect to a coordinate v actually depends on the choice of other coordinates. As a consequence, a complete delay embedding reconstruction of the dynamics using v which would allow phase to be computed still does not allow the v PRC to be computed. We give a coordinate-free definition of the PRC making this observation obvious. This leads to an experimental protocol: first collect an appropriate ensemble of measurements by intermittently controlling neuron voltage. Then, for any suitable current carrier dynamic postulated, we show how the ensemble can be used to compute the voltage PRC with that current carrier. The approach extends to many oscillators measured and controlled through a subset of their coordinates.

No datasets were generated or analysed during the current study.

1. In its general form, this is the same as “data-driven floquet analysis” representation of Revzen (2009) and Revzen and Kvalheim (2015), originally due to Floquet (1883).

2. Exponential convergence is sufficient for the existence of phase in oscillators, and necessary for the structural stability of the equations. Without structural stability the equations would be challenging to use in physical modeling.

3. Despite the fact that the closed 1-form defined by the right side of (9) is not exact on \({\mathbb {R}}^2{\setminus } \{{\textbf{0}}\}\), “\(d\theta \)” is still common notation for this 1-form. See for example Spivak (1971, p. 93).

4. Setting \(r = 1\) in \(\left. \frac{\partial {v}}{\partial {x}}\right| _{y}\) and viewing the resulting expression as a quadratic function of x with parameter y, the discriminant \(\Delta \) of this quadratic is \(\Delta = (e^\delta -1)^2C^2y^2-4(e^\delta -1)S^2\). Hence \(\Delta < 0\) with \(\delta > 0\) if and only if \((e^\delta -1)C^2y^2<4S^2\). This implies that, with \(\delta > 0\), \(\Delta < 0\) for all \((x,y)\in \{r=1\}\) if and only if \((e^\delta -1)C^2<4S^2.\) Since a quadratic has a real root if and only if its discriminant is nonnegative, the latter inequality is necessary and sufficient for the quantity \(\left. \frac{\partial {v}}{\partial {x}}\right| _{y}\) to be nowhere vanishing on \(\{r=1\}\). If it so happens that \(w\delta \not \in \frac{\pi }{2}+\pi {\mathbb {Z}}\), then this necessary and sufficient condition reads

   $$\begin{aligned} e^\delta -1< 4\tan ^2(w\delta ). \end{aligned}$$

   Since \(\tan ^2(w \delta ) = w^2 \delta ^2 + o(\delta ^5)\) as \(\delta \rightarrow 0\) while \(e^\delta - 1 = \delta + o(\delta )\), we see that, for this condition to hold, w must be very large if \(\delta \) is very small.

5. The notation \(\partial _z\), a vector, is not the partial derivative \(\frac{\partial }{\partial z}\); it is the derivation along z, i.e. given any path p(t), \(d/dt\,z(p(t)) = \langle \partial _z, \dot{p}(t) \rangle \). Thinking of z as scalar valued function, and ignoring issues of the presence or lack thereof of an underlying metric, \(\partial _z\) is the gradient of z, and the previous equation is an expression of the chain rule.

6. We note that many characterisations of neural oscillators in terms of phase do not apply a small voltage in order to estimate the PRC, but rather apply a spike (large voltage, short duration) to obtain a PRC which might be inconsistent with a first order approximation. Such an experiment is not considered here, although we note that what is required to calculate the induced phase change across the stability basin is some method of integrating the accumulated phase change resulting from travelling from the limit cycle to the perturbed state. Because the ensuing state change is typically “large”, i.e., of the order of the size of the limit cycle itself, infinitesimal approximations are unlikely to be accurate.

7. Even if we were assuming knowledge of the equations of motion, it seems likely to us that (17) is sufficiently complicated that constructing the delay coordinate transformation in closed form is intractable. We were able to construct the delay coordinate transformation in closed form for the toy example of Sect. 4.1 since its form was chosen for this very purpose.

Kvalheim was supported by ARO MURI W911NF-18-1-0327, the SLICE Multidisciplinary University Research Initiatives Program; and by ONR N00014-16-1-2817, a Vannevar Bush Faculty Fellowship sponsored by the Basic Research Office of the Assistant Secretary of Defense for Research and Engineering. Revzen was supported by ARO W911NF-14-1-0573, ARO MURI W911NF-17-1-0306, NSF CMMI 1825918, and the D. Dan and Betty Kahn Michigan-Israel Partnership for Research and Education. We thank John Guckenheimer for the discussion leading to the temporal 1-form approach to phase quantification.

Authors and Affiliations

1. Royal Veterinary College, London, UK

   Simon Wilshin

2. University of Maryland, Baltimore County, MD, USA

   Matthew D. Kvalheim

3. University of Michigan, Ann Arbor, MI, USA

   Shai Revzen

Contributions

All authors contributed to the writing of the manuscript. The formal proofs were undertaken by M.K. under supervision and with assistance of S.R. The numerical work was conducted by S.W under supervision and with the assistance of S.R.

Corresponding author

Correspondence to Shai Revzen.

Conflict of interest

The authors have no relevant financial or non-financial interests to disclose.

Communicated by Peter Thomas.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Appendix A. The Hopf oscillator satisfies Assumption 4

Consider the Hopf oscillator whose attractor is a trajectory around the unit circle with velocity 1.

Here \(\gamma _x(t) := \sin (t)\), leading to

as the non-autonomous y dynamics, which we expect to converge to \(y = \cos (t)\).

Changing to an error variable \(z := y-\cos (t)\), and shifting t by \(\pi /2\) for algebraic convenience gives

without loss of generality we can assume \(z>0\), since an initial negative z corresponds to a time shift of \(\pi \) in t, and we seek to prove a result for all t.

We need to show \(z\rightarrow 0\) for all initial z. Define \(L := \ln (z)\). We shall show that over a period of \(2\pi \), L decreases by at least some constant, thereby showing \(z\rightarrow 0\) exponentially. We do this by taking a separate bound for each half-cycle,

Appendix B. Estimating the period of the Fitzhugh-Nagumo system

Below is a listing of the python code used to calculate the period of the Fitzhugh-Nagumo oscillator. This oscillator is slow-fast and quite stiff, so some care is needed when integrating it; the default settings of the scipy odeint command do not produce reliable results.

It is useful to supply the Jacobian and to have a large number of intermediate steps. In addition, for the precision we required for the period in this work, a higher than default tolerance was needed.

Appendix C. Special cases estimating the infinitesimal PRC

There are multiple special cases which must be considered when estimating the partial derivatives used to calculate the infinitesimal phase response curve close to the limit cycle for this system.

We consider two models. First, we select a model where the dependent variable (the phase for the infinitesimal phase response curve, the un-transformed coordinates for the Jacobian) is a polynomial function of the independent variable (the transformed coordinate). We consider the linear and quadratic cases, and select from these models via an AIC. An example where the linear model was superior is included in Figure 7.

An example where a simple linear model was used to estimate the rate of change of the phase against change in dimensionless voltage at a fixed point on the limit cycle (black point in center of plot). We plotted the phases of the points at fixed delay voltages \(v_d\) across a range of voltages v (cyan crosses). We also plotted various model fits (dot-dash lines): linear (black), quadratic in v (red), and quadratic in phase \(\phi \) (blue). These imply potentially unequal derivatives which we indicate as lines with the associated slope (same color, solid thin lines). In this instance all fits give comparable estimates for the slope

Full size image

For some cases this model is deficient. There is a clear non-linear relationship between the independent variable and dependent variable, such that a quadratic in the dependent variable is a more appropriate model. Such cases can be readily identified by inspection, the non-linearity is obvious and there is frequently bi-modality in the distribution of the dependent variable. In such cases a quadratic in the dependent variable is instead used. This is illustrated in Figure 8.

An example where there is quadratic dependence of phase against delayed dimensionless voltage. The elements of this plot are as in Figure 7. The model with this quadratic dependence is clearly superior

Full size image

In some cases the transformation from the conventional coordinates of the Fitzhugh-Nagumo system to the delay coordinates contains a singularity which is very close to the limit cycle. This was handled in two ways.

In most cases the approach of this singularity was close, but not so close that a reasonable estimate of the gradient could not be obtained by simply excluding manually those points which came excessively close to the limit cycle. An example of such a case is illustrated in Figure 9.

An example where a singularity in the transformation to delay coordinates is present and close to the limit cycle. The elements of this plot are as in Figure 7. In this case, those points which are on the wrong side of the singularity can simply be excluded (right panel) and a better estimate of the partial derivative is obtained.

Full size image

For one test point the singularity came extremely close to the limit cycle and no reasonable estimate of the phase response was possible. Instead two additional points were selected near to this example and the phase response estimated at these locations instead. This is illustrated in Figure 10.

Here, the singularity in the transformation to delay coordinates is too close to the limit cycle, and no reasonable estimate of the partial derivative can be obtained, illustrated in the top panel. The elements of this plot are as in Figure 7. Instead, two points (bottom left and right panels) close by on the limit cycle are considered and estimates for the partial derivatives are obtained there

Full size image

Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.

Reprints and permissions