A Geometric Approach to Phase Response Curves and Its Numerical Computation Through the Parameterization Method

The phase response curve (PRC) is a tool used in neuroscience that measures the phase shift experienced by an oscillator due to a perturbation applied at different phases of the limit cycle. In this paper, we present a new approach to PRCs based on the parameterization method. The underlying idea relies on the construction of a periodic system whose corresponding stroboscopic map has an invariant curve. We demonstrate the relationship between the internal dynamics of this invariant curve and the PRC, which yields a method to numerically compute the PRCs. Moreover, we link the existence properties of this invariant curve as the amplitude of the perturbation is increased with changes in the PRC waveform and with the geometry of isochrons. The invariant curve and its dynamics will be computed by means of the parameterization method consisting of solving an invariance equation. We show that the method to compute the PRC can be extended beyond the breakdown of the curve by means of introducing a modified invariance equation. The method also computes the amplitude response functions (ARCs) which provide information on the displacement away from the oscillator due to the effects of the perturbation. Finally, we apply the method to several classical models in neuroscience to illustrate how the results herein extend the framework of computation and interpretation of the PRC and ARC for perturbations of large amplitude and not necessarily pulsatile.

This work has been partially funded by the Grants MINECO-FEDER MTM2015-65715-P, MDM-2014-0445, PGC2018-098676-B-100 AEI/FEDER/UE, the Catalan Grant 2017SGR1049, (GH, AP, TS), the MINECO-FEDER-UE MTM-2015-71509-C2-2-R (GH), and the Russian Scientific Foundation Grant 14-41-00044 (TS). GH acknowledges the RyC project RYC-2014-15866. TS is supported by the Catalan Institution for research and advanced studies via an ICREA academia price 2018. AP acknowledges the FPI Grant from project MINECO-FEDER-UE MTM2012-31714. We thank C. Bonet for providing us valuable references to prove Theorem 3.2 and A. Granados for numeric support. We also acknowledge the use of the UPC Dynamical Systems group’s cluster for research computing (https://dynamicalsystems.upc.edu/en/computing/).

Authors and Affiliations

1. Departament de Matemàtiques, Universitat Politècnica de Catalunya, Avda. Diagonal 647, 08028, Barcelona, Spain

   Alberto Pérez-Cervera, Tere M-Seara & Gemma Huguet

2. BGSMATH, Barcelona, Spain

   Alberto Pérez-Cervera, Tere M-Seara & Gemma Huguet

Corresponding author

Correspondence to Alberto Pérez-Cervera.

Communicated by Dr. Paul Newton.

Publisher's Note

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

Appendix: Numerical Algorithms

In this section, we review the numerical algorithms introduced in Haro et al. (2016) to compute the parameterization of an invariant curve \(\Gamma _A\) of a given map \(F:=F_A\) as well as the dynamics on the curve, i.e. \(f:=F|_{\Gamma _A}\). We present the algorithms in a format that is ready for numerical implementation, and we refer the reader to Haro et al. (2016), Pérez-Cervera et al. (2018) for more details.

The method to compute the invariant curve consists in looking for a map \(K:{\mathbb {T}} \rightarrow {\mathbb {R}}^2\) and a scalar function \(f:{\mathbb {T}} \rightarrow {\mathbb {T}}\) satisfying the invariance equation

In order to solve Eq. (63) by means of a Newton-like method, one needs to compute alongside the invariant normal bundle of K, denoted by N, and its linearized dynamics \(\Lambda _N\), which satisfy the invariance equation

Thus, the main algorithm provides a Newton method to solve Eqs. (63) and (64) altogether. More precisely, at step i of the Newton method one computes successive approximations \(K^i\),\(f^i\), \(N^i\) and \(\Lambda _N^i\) of K, f, N and \(\Lambda _N\), respectively. The algorithm is stated as follows:

Algorithm A.1

Main Algorithm to Solve Equations (63), (64). Given \(K(\theta )\), \(f(\theta )\), \(f^{-1}(\theta )\), \(N(\theta )\) and \(\Lambda _N(\theta )\), approximate solutions of Eqs. (63) and (64), perform the following operations:

1. 1.

   Compute the corrections \(\Delta K(\theta )\) and \(\Delta f(\theta )\) by using Algorithm A.2.

2. 2.

   Update \(K(\theta ) \leftarrow K(\theta ) + \Delta K(\theta ) \quad \quad f(\theta ) \leftarrow f(\theta ) + \Delta f(\theta )\).

3. 3.

   Compute the inverse function \(f^{-1}(\theta )\) using Algorithm A.3.

4. 4.

   Compute \(DK(\theta )\) and \(Df(\theta )\).

5. 5.

   Compute the corrections \(\Delta N(\theta )\) and \(\Delta _N (\theta )\) by using Algorithm A.5.

6. 6.

   Update \(N(\theta ) \leftarrow N(\theta ) + \Delta N(\theta ) \quad \quad \Lambda _N(\theta ) \leftarrow \Lambda _N(\theta ) + \Delta _N (\theta )\).

7. 7.

   Compute \(E = F \circ K - K \circ f\) and repeat steps 1–6 until E is smaller than the established tolerance.

Next, we provide the sub-algorithms for Algorithm A.1.

Algorithm A.2

Correction of the Approximate Invariant Curve. Given \(K(\theta )\), \(f(\theta )\), \(f^{-1}(\theta )\)\(N(\theta )\) and \(\Lambda _N(\theta )\), approximate solutions of Eqs. (63) and (64), perform the following operations:

1. 1.

   Compute \(E(\theta ) = F(K(\theta )) - K(f(\theta ))\).

2. 2.

   Compute \(P(f(\theta ))= \Big (DK(f(\theta )) | N (f(\theta ))\Big )\).

3. 3.

   Compute \(\eta (\theta ) = \begin{pmatrix} \eta _T(\theta ) \\ \eta _N(\theta ) \end{pmatrix} = -\,(P(f(\theta )))^{-1}E(\theta )\).

4. 4.

   Compute \(f^{-1}(\theta )\) using Algorithm A.3.

5. 5.

   Solve for \(\xi \) the equation \(\eta _N(f^{-1}(\theta ))= \Lambda _N(f^{-1}(\theta ))\xi (f^{-1}(\theta )) - \xi (\theta )\) by using Algorithm A.4.

6. 6.

   Set \(\Delta f(\theta ) \leftarrow - \eta _T(\theta )\).

7. 7.

   Set \(\Delta K(\theta ) \leftarrow N(\theta )\xi (\theta )\).

Algorithm A.3

Refine\({\mathbf{f}^{\mathbf{-1}}(\varvec{\theta })}\). Given a function \(f(\theta )\), its derivative \(Df(\theta )\) and an approximate inverse function \(f^{-1}(\theta )\), perform the following operations:

1. 1.

   Compute \(e(\theta ) = f(f^{-1}(\theta )) - \theta \).

2. 2.

   Compute \(\Delta f^{-1}(\theta ) = -\, \frac{e(\theta )}{Df(f^{-1}(\theta ))}\).

3. 3.

   Set \(f^{-1}(\theta ) \leftarrow f^{-1}(\theta ) + \Delta f^{-1}(\theta )\).

4. 4.

   Repeat steps 1–3 until \(e(\theta )\) is smaller than a fixed tolerance.

Algorithm A.4

Solution of a fixed point equation. Given an equation of the form \(B(\theta ) = A(\theta )\eta (g(\theta )) - \eta (\theta )\) with A, B, g known and \(\Vert A \Vert < 1\), perform the following operations:

1. 1.

   Set \(\eta (\theta ) \leftarrow B(\theta )\).

2. 2.

   Compute \(\eta (g(\theta ))\).

3. 3.

   Set \(\eta (\theta ) \leftarrow A(\theta )\eta (g(\theta )) + \eta (\theta )\).

4. 4.

   Repeat steps 2 and 3 until \(|A(\theta )\eta (g(\theta ))|\) is smaller than the established tolerance.

Algorithm A.5

Correction of the stable normal bundle. Given \(K(\theta )\), \(f(\theta )\), \(N(\theta )\) and \(\Lambda _N(\theta )\), approximate solutions of Eqs. (63) and (64), perform the following operations:

1. 1.

   Compute \(E_N(\theta ) = DF(K(\theta ))N(\theta ) - \Lambda _N(\theta )N(f(\theta ))\).

2. 2.

   Compute \(P(f(\theta )) = (DK(f(\theta )) \quad N(f(\theta )))\).

3. 3.

   Compute \(\zeta (\theta ) = \begin{pmatrix} \zeta _T(\theta ) \\ \zeta _N(\theta ) \end{pmatrix} = -(P(f(\theta )))^{-1}E_N(\theta )\).

4. 4.

   Solve for Q the equation \(Df^{-1}(\theta ) \zeta _T(\theta )=Df^{-1}(\theta ) \Lambda _N(\theta )Q(f(\theta )) - Q(\theta )\) by using Algorithm A.4.

5. 5.

   Set \(\Delta _N (\theta ) \leftarrow \zeta _N(\theta )\).

6. 6.

   Set \(\Delta N(\theta ) \leftarrow DK(\theta )Q(\theta )\).

Remark A.6

Since our functions are defined on \({\mathbb {T}}\), we will use Fourier series to compute the derivatives and composition of functions.

Of course, the main algorithm requires the knowledge of an approximate solution for Eqs. (63), (64). In our case, we can always use the limit cycle of the unperturbed system as an approximate solution for the invariant curve. However, for the normal bundle we cannot use the one obtained from the unperturbed limit cycle \(\Gamma _0\). The following algorithm provides an initial seed for Algorithm A.1.

Algorithm A.7

Computation of Initial Seeds. Given a planar vector field \({\dot{x}} = X(x)\), having an attracting limit cycle \(\gamma (t)\) of period T, perform the following operations:

1. 1.

   Compute the fundamental matrix \(\Phi (t)\) of the first variational equation along the periodic orbit.

2. 2.

   Obtain the characteristic multiplier \(\lambda \ne 1\) and its associated eigenvector \(v_{\lambda }\) from \(\Phi (T)\).

3. 3.

   Set \(N(\theta ) \leftarrow \Phi (T\theta )v_{\lambda }e^{\lambda \theta }\).

4. 4.

   Set \(\Lambda _N(\theta ) \leftarrow e^{\frac{\lambda T'}{T}}\).

5. 5.

   Set \(K(\theta ) \leftarrow \gamma (T\theta ), \quad \quad DK(\theta ) = TX(\gamma (T\theta ))\).

6. 6.

   Set \(f(\theta ) \leftarrow \theta + \frac{T'}{T}, \quad \quad f^{-1}(\theta ) \leftarrow \theta - \frac{T'}{T}, \quad \quad Df(\theta ) = 1\).

Given a family of maps \(F_A\) such that the solution for \(F_0\) is known (see Algorithm A.7), it is standard to set a continuation scheme to compute the solutions for the other values of A. Thus, assuming that the solution for \(A=A^{*}\) is known, one can take this solution as an initial seed to solve the equations for \(F_{A^*+h}\), with h small, using the Newton-like method described in Algorithm A.1. However, one can perform an extra step to refine the initial seed values \(K_{A^ {*}}\) and \(f_{A^*}\) for \(F_{A^{*}+h}\), described in Algorithm A.8.

Algorithm A.8

Refine an Initial Seed. Given \(K_A(\theta )\), \(f_A(\theta )\), \(N_A(\theta )\) and \(\Lambda _{N,A}(\theta )\), solutions of Eqs. (63), (64) for \(F=F_A\), perform the following operations:

1. 1.

   Compute \(E(\theta ) = \frac{\partial F_A}{\partial A}(K_A(\theta ))\).

2. 2.

   Compute \(\eta (\theta ) = \begin{pmatrix} \eta _T(\theta ) \\ \eta _N(\theta ) \end{pmatrix} = -\,(P(f_A(\theta )))^{-1}E(\theta )\).

3. 3.

   Solve for \(\xi \) the equation \(\xi (\theta ) = \Lambda _N(f^{-1}_A(\theta ))\xi (f^{-1}_A(\theta )) - \eta _N(f^{-1}_A(\theta ))\) by using Algorithm A.4.

4. 4.

   Set \(K_{A+h}(\theta ) \leftarrow K_A(\theta ) + N_A(\theta )\xi (\theta ) \cdot h, \quad \quad f_{A+h}(\theta ) \leftarrow f_A(\theta ) - \eta _T(\theta ) \cdot h\).

Remark A.9

The term \(\frac{\partial F}{\partial A}(K_A(\theta ))\) is computed using variational equations with respect to the amplitude.

The numerical continuation scheme for our problem is described in the following algorithm.

Algorithm A.10

Numerical Continuation. Consider a family of maps \(F_A\), such that \(F_0\) is the time-T map of a planar autonomous system having a hyperbolic attracting limit cycle of period T. Perform the following operations:

1. 1.

   Compute solutions \(K_0(\theta )\), \(f_0(\theta )\), \(N_0(\theta )\) and \(\Lambda _{N,0}(\theta )\) of Eqs. (63), (64) for \(F=F_0\) using Algorithm A.7.

2. 2.

   Set \(A=0\).

3. 3.

   Using \(K_A\), \(f_A\), \(N_A\), \(\Lambda _{N,A}\) compute an initial seed for \(F_{A+h}\) using Algorithm A.8

4. 4.

   Find solutions of Eqs. (63), (64) for \(F=F_{A+h}\) using Algorithm A.1.

5. 5.

   Set \(A \leftarrow A+h\).

6. 6.

   Repeat steps 3–5 until A reaches the desired value.

Reprints and permissions