Inferring the phase response curve from observation of a continuously perturbed oscillator

Phase response curves are important for analysis and modeling of oscillatory dynamics in various applications, particularly in neuroscience. Standard experimental technique for determining them requires isolation of the system and application of a specifically designed input. However, isolation is not always feasible and we are compelled to observe the system in its natural environment under free-running conditions. To that end we propose an approach relying only on passive observations of the system and its input. We illustrate it with simulation results of an oscillator driven by a stochastic force.

(a)  Figure S1. A reconstruction for a phase oscillator driven by white noise (analogous to Fig. 1 in the main text, (a, c) in blue differences of the Wiener process ∆W instead of signal p). The driving strength is ε Z = 1. (a) Figure S2. An example of a reconstruction with strong driving, Morris-Lecar neuron (similar to Fig. 1 in the main text).

2/7
The Hindmarsh-Rose neuron: example on bursting To further support the argument that the data our method requires is general we perform a reconstruction on a bursting neuron. We use the Hindmarsh-Rose neuronal model 2 : the parameters are: I = 1.28, a = 1, b = 3, c = 1, d = 5, r = 0.0006, s = 4 and x R = −1.6. In this regime the neuron is bursting periodically. For the reconstruction we use parameters ε Z = 0.01, τ = 0.01, t sim = 1000, dt = 0.0001 and N = 15. The beginning of the period is determined by the time of the first spike in a burst (this is arbitrary, any of the spikes could be considered for this role). See Fig. S4 for a reconstruction depiction. The PRC of this bursting cycle has a very sharp feature around the phase ϕ ≈ 1 where the bursting stops. This is hard to capture with our method since we take the Fourier representation of the PRC and steep features require several harmonics to be expressed. Taking many harmonics can be impractical since each harmonic yields an unknown constant needed to be optimized, which can make the data requirements considerably big. With N = 15 Fourier harmonics that feature is to some degree smeared in phase, but still apparent.
(a) Figure S4. An example of a reconstruction using bursts as the events with equal phase, Hindmarsh-Rose neuron (similar to Fig. 1 in the main text). Parameters are ε Z = 0.01 and τ = 0.01.

3/7
The van der Pol oscillator: varying the number of Fourier harmonics N We explore the effect of using different numbers of Fourier harmonics N. The PRC of the oscillator for the chosen parameters can not be captured well by a single harmonic term but we also know that using more terms means more unknowns to fit and therefore a necessity for more data to constrain them. By increasing the number of Fourier harmonics N while keeping the length of the time series t sim constant we therefore expect the error to have a minimum, see

Choice of a proper Poncaré section
In the case of a forced oscillator, the correct Poincaré section for determining instants of same phase corresponds to an isochrones surface. Here we show the isochrone structure of the two oscillators used in the main text, van der Pol and Morris-Lecar, see Fig. S6. The isochrones were computed using the backward integration method 3 . (a)

4/7
Now suppose we only measure one variable x(t) but want to embed our signal in 2 dimensions and then determine the periods using a Poincaré section at an arbitrary angle α. This is done by rotating the embedded limit cycle and calculating the threshold value accordingly, a depiction can be seen in Fig. S7. Figure S7. A schematic depiction of period determination with the use of signal embedding and inclined surfaces of section, as used for Fig. 6 in the main text. In (a) the measured signal x(t) and in (b) the proxi variablex(t) (in our case we used the first derivative). In (c) the embedded signal (with green) as well as the chosen surface of section depicted with a straight red line. In (d) the embedded signal rotated by α so that the chosen surface of section becomes paralel to the horizontal axis. In (e) the vertical projection of the rotated embedded signal which is to be thresholded.

5/7
Using a similar search to that performed in the main text in Fig. 6, we can estimate the isochrone structure in the vicinity of the limit cycle, see Fig. S8. First, we embed our signl x(t) in two dimensions (see Fig. S7 above). Then, considering the embedded signal in polar cooridnates, we average the radius variable over the angle variable with finite binning to obtain an approximation of the unpertubed limit cycle. On the obtained limit cycle approximation we choose points uniformly distributed in angle and for each point try several Poincaré sections at different angles. The angle corresponding to the lowest error should closely match the local isochrone.