Inverse Approach to Chronotaxic Systems for Single-variable Time Series

Following the development of a new class of self-sustained oscillators with a time-varying but stable frequency, the inverse approach to these systems is now formulated. We show how observed data arranged in a single-variable time series can be used to recognize such systems. This approach makes use of time-frequency domain information using the wavelet transform as well as the recently developed method of Bayesian-based inference. In addition, a set of methods, named phase fluctuation analysis, is introduced to detect the defining properties of the new class of systems by directly analyzing the statistics of the observed perturbations.We apply these methods to numerical examples but also elaborate further on the cardiac system.


I. INTRODUCTION
A new class of oscillatory systems which are characterized by a time-varying but stable frequency has been recently developed [1].These chronotaxic systems provide a way of modeling stable dynamics that had not previously been considered.In particular, they offer a novel way of generating oscillations with time-dependent frequencies that are resistant to external perturbations.As shown in the previous work [1], chronotaxic systems contain a deterministic element that is shielded from noise and other external effects but they can still appear stochastic to other analytical frameworks.The theory therefore greatly simplifies the description of complex dynamics and also reveals stability where there may have previously been thought to be none.
As we begin to apply the theory of chronotaxic systems to the real world, in this paper we now ask the important question: How can we detect a chronotaxic system from observations?By answering this question we intend to formulate the analysis of chronotaxic systems in the inverse approach.
The inverse approach to dynamical systems has been an active field since the introduction of time and frequency domain techniques [2] as well as phase space methods [3][4][5].The stochastic approach to dynamical systems has also seen much activity [6].Such methods reveal properties of an underlying system from the information contained in its observable variables, allowing classification of the dynamics.
However, one of the greatest challenges that remains is the case of inverse problems of nonautonomous systems, i.e., those that are inherently time dependent, which includes chronotaxic systems.Such systems are widely prevalent in nature, from life [7] to climate dynamics [8], where the situation is inevitably one that is thermodynamically open and subject to time-dependent effects.They also commonly appear as a consequence of observing only a single variable in coupled systems, which effectively turns the other variables into a timedependent influence.Being able to analyze single-variable data of nonautonomous systems therefore has wide-reaching implications.
The development of methods that can be applied to timevarying dynamics has been an extremely active field in recent years.Of particular note are wavelet-based methods such as wavelet phase coherence [9], bispectrum [10], and also the synchrosqueezed transform [11].Methods to fit nonautonomous models to data have also been investigated [12], including a novel technique based on Bayesian inference [13,14] which is able to track time-dependent system parameters [15,16].Moreover, these methods have been successfully applied to many areas including blood flow dynamics [17,18], aging [19,20], neuroscience [21], and climate science [8].
The following will detail the procedure for detecting chronotaxicity for the nontrivial case where only a singlevariable time series of a system is observed.In Secs.II and III the background to this work is presented, including an overview of how chronotaxic systems are defined and their relation to established techniques used in the inverse approach to dynamical systems.Sections IV A-IV D detail the framework of analysis for chronotaxic systems where a singlevariable time series is observed, as well as the restrictions for when different parts of this framework can be applied.In Sec.V A the application of the methods is demonstrated for numerically generated phase oscillators.Lastly, the framework is applied to a physiological recording of the electrical activity of the heart in Sec.V B and the summary, conclusions, and current outlook are presented in Sec.VI.

II. MODEL OF CHRONOTAXIC SYSTEMS
Following Refs.[1,22], chronotaxic systems are defined as a subclass of nonautonomous dynamical systems, which means their dynamics is governed by an independent variable or time-dependent component.This is achieved by the use of unidirectionally coupled differential equations, known as a skew-product flow [23][24][25], master-slave configuration [26], or as drive and response systems [27], ẋ = g(x,p).
Here x(t) contains the dynamics of the observed nonautonomous system.One of the most important properties of such a system is that its trajectory depends not only on the current time t and initial conditions x 0 as in autonomous systems, but also depends on the initial time t 0 .
The structure defined in (1) is required in order for an oscillatory system x to have stable dynamics, which are typically characterized by an amplitude and phase.In terms of amplitude, stable dynamics can be obtained through an autonomous limit cycle, which corresponds to a stationary stable fixed point in the amplitude.This can be expressed by the action of a negative Lyapunov exponent, where the initial displacement from the limit cycle and subsequent amplitude perturbations exponentially decay.As a result, unperturbed trajectories converge to the limit cycle over time.Similarly, higher-dimensional systems which are nonchaotic can also have stable amplitudes as they are decomposable into individual limit cycles [28].
Besides the amplitude, the oscillatory dynamics discussed above is described by the coordinate along the limit cycle which is known as the phase.In autonomous systems the phase of a system on a limit cycle is not stable, i.e., its Lyapunov exponent is zero.The natural frequency of the system describes the rate at which the unperturbed phase changes, but the phase can be pushed in either direction by an external force and the underlying rate of change will not be affected.
This means that while the natural frequency itself does not change, the actual frequency of the oscillations can be slowed if perturbations push the system in the opposite direction to the change in phase, or increased if perturbations push the system in the same direction.
In order to make the phase dynamics stable and to consequently have a stable frequency, there has to be a steady state x A on the limit cycle which trajectories converge to.In other words, the Lyapunov exponent of the phase dynamics has to be negative rather than zero.However, once the system has reached x A it would remain stationary, meaning that the oscillations would cease.Therefore, to maintain oscillations the system must also be time dependent (nonautonomous) so that the steady state moves around the limit cycle.In this case the frequency of the oscillations is fixed to the rate at which x A (t) moves.Furthermore, the frequency can change over time but still remains stable.
In general, x A (t) is a time-dependent point attractor or steady state for all trajectories, satisfying the condition of invariance the condition that all trajectories converge to it in the pullback limit [23,24], and in the forward limit, The motion of the point attractor x A (t) along the limit cycle results in a special object-a chronotaxic limit cycle.The above overview of the theory of chronotaxic systems sufficiently describes the aspects related to the inverse approach.However, extensions to this theory are covered in Ref. [22], including a discussion of its applications.

III. PREVIOUS INVERSE APPROACHES AND NONAUTONOMOUS SYSTEMS
Following the brief introduction of chronotaxic systems we now focus on tackling such systems in the inverse approach.Since chronotaxic systems are dominated by deterministic dynamics, the conventional approach to their analysis would involve the reconstruction of the attractor in phase space and the properties of the system would then be determined using the ergodic theory [29].However, these methods fall down when applied to chronotaxic systems as they are nonautonomous.
One of the most important issues in analyzing any nonautonomous system in this way is that many of the common timeindependent analytical methods do not give useful information.In most cases the frequency distribution of the system is nonstationary, which means the Fourier transform contains a mixture of the power spectra at different times that becomes difficult to interpret.In addition, other representations such as probability distributions and phase space embeddings give a time-averaged view of the system that makes it appear high dimensional and stochastic [30].This also excludes the possibility of using probabilistic approaches based on Granger causality which are otherwise very useful when applied to stationary data [31].
The usual way of dealing with nonstationary data is to perform the same analyses but only on part of the time series by using a moving time window.However, this creates a number of issues.In the case of the Fourier transform the window size limits the frequency resolution of the Fourier spectrum, which has a huge impact on the detection of low-frequency oscillations (an effect known as the Gabor limit [2]).In phase space, the system may be unable to explore all parts of the attractor within a given window, resulting in an incomplete reconstruction.
For the problem of analyzing nonstationary data in the frequency domain the solution is to use the wavelet transform, which is defined in the following section.However, in the case of phase space analysis, the methodology used to reconstruct a time-dependent attractor can result in "blurring," where trajectories either overlap or do not match up correctly.One way in which this happens is from the fact that the data used in the embedding is taken from some time interval, where data at the beginning of the interval corresponds to the attractor at a different time from the one derived from data at the end of the interval.Another effect is more subtle as it is related to the parameters used in the time delay embedding theorem [32,33].In particular, the time delay used ensures that the phase space dimensions are independent and the parameter is constant for the entire embedding.When the attractor is nonautonomous, and since the parameter is not adaptive, this means that the dependence of the dimensions may change in time, resulting in an attractor that appears more blurred and time dependent than it actually is.
Recent work has improved the embedding procedure for complex time-dependent systems such as the brain and cardiovascular system, revealing the direction of coupling between subsystems [34].However, the application of the ergodic theory to nonautonomous systems is still restricted.

IV. THE INVERSE APPROACH TO CHRONOTAXIC SYSTEMS
As a subclass of nonautonomous systems, it is not appropriate to analyze chronotaxic systems using many of the typical inverse approaches to dynamical systems.However, in contrast to the general case of nonautonomous systems, the characteristics provided in their definition can be taken advantage of, particularly those corresponding to the phase dynamics.By taking these defining properties into consideration a new framework of analysis is now presented.

A. Separation of phase and amplitude
One of the requirements for chronotaxic systems is that they possess a stationary limit cycle [1].The existence of a limit cycle means that the system can be transformed to polar coordinates where a phase and amplitude of the oscillation can be defined.These correspond to angular direction and radial distance from a reference point, respectively.
The amplitude dynamics of a chronotaxic system corresponds to the convergence of the system to the limit cycle, which is influenced only by a negative Lyapunov exponent and external perturbations.On the other hand, the phase dynamics corresponds to convergence to the time-dependent point attractor, which is characterized by a negative Lyapunov exponent and external perturbations, but also by the motion of the point attractor itself.The influence of this point attractor in the phase dynamics is the focus of investigation of chronotaxic systems in the inverse approach, which means the separation of the phase from the amplitude is an important first step.
In the time-frequency domain, the polar coordinates can be extracted from an oscillation if it has one continuous frequency in time."Continuous" is specified here because rapid jumps in frequency cannot be tracked due to the time-frequency resolution limit.This restriction also requires that only one oscillatory component is present in the time series or that multiple components can be extracted via decomposition, in which case each oscillation requires a separate transformation to polar coordinates.The exception to this is when the oscillations are in fact harmonics caused by nonlinearities in the system but corresponding to the same limit cycle oscillation [9].
Consider a time series f (t) of length L which contains a limit cycle oscillation.In an ideal situation the time series would only contain the dynamics of a sinusoidal-like oscillation, in which case the Hilbert transform could be applied to give an almost perfect extraction of the phase [35,36].However, when this is not the case it is fairly straightforward to calculate the frequency and amplitude of the oscillation for every point in time by using the continuous wavelet transform [37] where (s,t) is the mother wavelet, which is scaled according to the parameter s to change its frequency distribution and time shifted according to t.The oscillation can be traced in W T (s,t) using either a ridge-extraction method [38,39] or synchrosqueezed wavelet transform [11].While the amplitude dynamics are given by the amplitude of W T (s,t) at the positions of the oscillation in the s-t plane, extraction of the exact phase of the oscillation requires more careful consideration.
The extraction methods [11,38,39] are used to estimate the instantaneous frequencies of the oscillatory components in a time series, allowing identification of harmonics which can be used to determine the intracycle dynamics.The phase can then be found by integrating over the instantaneous frequency in time.This can be done regardless of the wavelet basis used, but the phase can also be found another way by using complex basis such as the Morlet wavelet [40], where f 0 is a parameter known as central frequency which defines the time and frequency resolution [37].Here, the corresponding frequency of the wavelet is given by 1/s.The phase of the observed system is then arg[W T (s,t)], where s and t denote the positions of the oscillation in the s-t plane.

B. Extracting the phase of the point attractor
Once the phase has been found the goal is to investigate the phase component of the oscillations, α x , to see if it is attracted to a time-dependent point attractor, α A x .Considering the case where the phases are from a chronotaxic system, if α x is left unperturbed then it would eventually converge to α A x .This means that one can determine the dynamics of α A x by finding α x when the system is unperturbed.The problem is that when given only the time series of the perturbed α x it is not always possible to see what part of dynamics comes from the perturbations and what is due to α A x as both are freely defined functions.In other words, the effect of a change in α A x on α x can be alternatively represented as deviations caused by the (noise) function ξ (t), Here f (α x ,α A x ,t) is the coupling function and ξ (t) is the function or "noise" describing the perturbations.In the case where ξ (t) = 0, α A x is the point attractor for the phase α x .In order to distinguish between the effects of α A x and ξ (t) we assume that the dynamics of α A x is confined to time scales larger than a single cycle.We also assume that the noise is either weak or comparable in magnitude to f (α x ,α A x ,t).These assumptions are valid in many real-world systems such as the cardiovascular system, where the low-frequency dynamics plays an important role [41].
For chronotaxic systems, an estimate of α A x , denoted α A * x , can therefore be found by filtering out high-frequency components in the dynamics of α x .At the same time, such a filter should not smooth over the dynamics of α A x .The possibility of phase slips also means that α A * x cannot be obtained by simply performing a moving average over the phase.Instead, smoothing over the frequency extracted from the wavelet transform provides the estimated angular velocity αA * x , which can in turn be integrated over to give α A * x .The procedure of extracting from a single time series both the estimate for the phase of the system, α * x , and the phase of the attractor, α A * x , is demonstrated in Fig. 1.In cases where phase x found by ridge extraction in the wavelet transform and then smoothing using a 2 s moving average (solid line) and the actual value of αA x from the analytic solution (dashed line).(d) The phase α A * x found by integrating over αA * x (solid line) and the actual dynamics of α A x from the analytic solution (dashed line).The offset of the phases is due to the estimated phase starting at a later time as αA x cannot be defined for the interval where the wavelets run over the ends of the time series.
slips are uncommon both α * x and α A * x can be extracted from the same wavelet transform.In some cases, however, when phase slips occur with almost every cycle, a high time resolution is required to track the phase during phase slips while a highfrequency resolution is needed to follow the driving frequency of the oscillation.Consequently, the central frequency of the wavelet transform needs to be small to extract α * x and large to extract α A * x .Considering now that the system is not chronotaxic, while α * x will still give an estimate of the phase of the system, α A * x will simply be a trend in the dynamics of α x and not correspond to a point attractor.The methods presented in the next sections show how such a case can be identified and, hence, determine when the phases are from a chronotaxic system.

C. Bayesian-based inference
One approach to determine whether or not an observed system is chronotaxic is to analyze the coupling function between α x and α A x , as shown in (7).Subsequent analysis on the inferred function can reveal the dynamics of the system in the absence of noise and provide evidence for chronotaxic behavior.
However, calculating the coupling function is a difficult task since time dependence plays an important role in chronotaxic systems.This can arise from changes in the motion of α A x but also in changes to the parameters describing the coupling.One solution to this problem is to use an approach based on Bayesian inference.
The Bayesian-based approach is optimized in a similar way to how the wavelet transform is optimized to track changes in frequency, but instead of frequency it is the dynamical equations of the system that are tracked in time [13,15,16,[42][43][44].This method assumes that the dynamics of two coupled phases is described by where ω i is the natural frequency of the oscillation, f i (α i ) are the self-dynamics of the phase, g i (α i ,α j ) are the cross couplings, and ξ (t) is a two-dimensional white Gaussian noise, ξ i (t)ξ j (τ ) = δ(t − τ )E ij .By assuming that the coupling terms are periodic functions, f i (α i ) and g i (α i ,α j ) are modeled using the Fourier bases and where k,r,s = 0.In practice, these bases are truncated to a finite number of functions, denoted A i,k (α i ,α j ).The corresponding parameters from ãi and bi then form the parameter vector c (i) k .In the method the inference of these parameters makes use of Bayes' theorem, where p X (M|X ) is the posterior probability distribution and (X |M) is the likelihood function for the values of the model parameters M given the data X , and p prior (M) is a prior distribution.The method assumes that probability functions of the inferred parameters take the form of a multivariate Gaussian distribution and that the noise, given by the matrix E ij , is white and Gaussian.An estimate can be found from the stationary point in the minus log-likelihood function where h is the time between samples in the time series and is assumed small enough to integrate the model given in (8) and summation over the repeated indices k,l,i,j is implicit.
In the case of the Euler midpoint scheme the time derivative of the phases is defined as α = (α n+1 − α n )/h, while α * = (α n+1 + α n )/2.Furthermore, the matrix A contains the base functions as shown in ( 9) and (10), while E accounts for the noise.The corresponding prior probability for the parameters is given by the multivariate Gaussian distribution where M is the number of parameters arranged in c (l) k and −1 is the covariance matrix of the parameters.The linearity of S means that if the prior probability distribution is assumed to be Gaussian, then the posterior distribution must also be Gaussian.Consequently, the form of the algorithm does not change when the posterior from a previous iteration is used to construct the prior for the next.To explicitly show the inference of the stationary point in S and the posterior distribution using Bayes' theorem, Eqs. ( 11)-( 13) are decomposed into the following: Here, (i,j ) prior exploits the information from the posterior matrix (i,j ) which is inferred from the previous time window (only for the initial window evaluation a flat prior (i,j ) prior = 0 and c prior = 0 is used).In this way, information is allowed to propagate between windows, facilitating the inference of parameters so that they become more accurate with time.However, to account for time variations in the values of the parameters the prior takes the form of a convolution between the posterior of the previous window and a user-defined diffusion matrix describing the change in c (i)  k or (i,j ) prior [15].The standard deviation corresponding to the diffusion of the parameters from (i,j ) prior is assumed to be a known fraction of the parameters themselves, σ (i)  k = pc (i) k , where p is a free parameter known as the propagation constant.This modification allows the method to track the change in the couplings over time.
In chronotaxic systems, α A x and α x are synchronized when α x is not influenced by noise.Therefore, determining if this synchronization occurs in the observed system provides an indication that it is chronotaxic.In particular, the synchronization between α A * x and α * x can be detected using the mapping procedure proposed in Refs.[15,16].In brief, two coordinates, x , are defined and bounded within the interval [0,2π ] so that the coordinates become "wrapped" on a circle.The system is integrated with the inferred parameters using the standard fourth-order Runge-Kutta algorithm from the initial condition ζ = 0, where it is assumed dζ (t)/dt| ζ =0 > 0 so that all trajectories travel in the same direction.A Poincaré section at ζ = 0 is then defined to detect the evolution of the trajectories after a complete cycle in ζ .The position of the trajectories after are designated by the map M, i.e., ψ n+1 = M(ψ n ), where ψ n are the initial conditions at ζ = 0.The existence of points where M(ψ n ) − ψ n = 0 indicate periodic orbits.In particular, the system is found to be chronotaxic in the case where two such fixed points exist; one stable, the other unstable.A chronotaxic index I chrono can then be defined so that I chrono = 1 when this condition is fulfilled, otherwise I chrono = 0.
The above method determines whether the inferred parameters of the coupling functions result in synchronization and, hence, provide evidence that the system is chronotaxic.However, rather than originating from the stability in the coupling function of α * x , this could also imply that the synchronization is a result of stability in the coupling function of α A * x .Therefore, in order to understand which of the coupling functions is responsible, it is necessary to calculate the direction of coupling.This is defined in Ref. [45] by where are the Euclidian norms of the parameters.Here, the odd parameters correspond to the coupling terms inferred for α 1 in the direction 2 → 1, while the even parameters correspond to the coupling terms inferred for α 2 in the direction 1 → 2. Note that the definition provided here for the direction of coupling is not unique.One can also calculate 12 and 21 by integrating over the coupling functions with respect to each phase [45].The direct use of 2  12 and 2 21 without calculating D has also been previously adopted [46].In the same paper the use of interval estimates, evaluated from both the parameters and their variances, were used to improve the inference of the couplings.This can easily be applied to the method shown here by using all of the elements in the covariance matrix −1 so that the inference also includes the variances of the parameters [47].
In summary, after α * x and α A * x have been extracted from the time series, Bayesian-based inference is applied to find the values of the coupling parameters for each phase in time.These parameters are then used to numerically map trajectories of the system to identify fixed points and calculate the index I chrono to determine if the system is chronotaxic or not.To check if the fixed points are really due to a coupling from α A * x to α * x , the directionality index D is calculated from the inferred parameters.If D points in the direction from α A x to α x and I chrono = 1, then this suggests that α x is a stable phase and is, hence, chronotaxic.

D. Phase fluctuation analysis
When the phases α x and α A x are almost perfectly synchronized the information flow between the two time series is significantly reduced, which can cause spurious inference in the above method but also in general phase dynamics approaches [36,48].For example, consider the system where the phase α x is neutrally stable with a zero-valued Lyapunov exponent.To distinguish between this system and a genuine chronotaxic system, an alternative method is needed.
The Bayesian approach searches directly for the "causes" of a point attractor in the dynamical equations of the system.However, another way to determine the existence of a point attractor is by using an indirect method that looks only at the "effects."By obtaining α A * x from the wavelet transform one aims to have the trajectory of the phase α x in the absence of perturbations.Calculating the difference x therefore provides information about the divergence from the point attractor due to these perturbations.
The restriction here is that α x only gives the exact divergence from the point attractor if α x remains within the same cycle.The discrepancy between α x and the actual divergence therefore changes by 2π with each phase slip.The wavelet transform is also blind to the direction of these phase slips so the true distance from the original "unperturbed" trajectory is also unknown.For a chronotaxic system, assuming the perturbations are due to an uncorrelated Gaussian process, two important effects can be observed in α x : (i) Plateaus where α * x runs more or less parallel to α A * x .These occur when the perturbations are small enough so that the system remains within the same cycle of the attractor.The influence of the attractor causes perturbations to decay, which means the divergence from the attractor is similar to the original Gaussian process.In contrast, when the phase is neutrally stable the perturbations are integrated over and therefore highly correlated in time, resulting in a random walk (i.e., Brownian noise).
(ii) Phase slips where large perturbations cause the system to either lag behind the attractor or move forward enough so that it becomes influenced by the attractor in the adjacent cycle.For a neutrally stable phase, phase slips do not occur-large perturbations can still cause α x to change by 2π but they are part of a continuous probability distribution.
A test for chronotaxicity can be performed by using a moving average (or similar time domain smoothing function) on α x over a time scale τ to define the trend, denoted α τ x .The distance between points in the delayed and original versions of this trend d α τ x (t) = α τ x (t + τ ) − α τ x (t) then provides information about the perturbations to the system over that time scale.However, because the direction of the perturbations that cause phase slips is essentially unknown, in practice it is best to consider only the magnitude |d α τ x |.If the system is nonchronotaxic and α x is a random walk, then the distribution of points in |d α τ x | should produce a one-sided Gaussian for all τ .On the other hand, for a chronotaxic system the time series will appear Gaussian only for time scales below those corresponding to phase slips.As soon as the averaging is performed over time scales close to the time taken for the phase to slip, the distribution will begin to shift and obtain a new peak.
If phase slips are rare then the extra peak observed in |d α τ x | might not be statistically significant.However, this also means that the system remains close to the attractor, in which case effect (i) from above dominates.This means that the perturbations should resemble noise with little correlation in time, as opposed to integrated noise with long-range correlation.In order to test this, detrended fluctuation analysis (DFA) can be performed on the time series [19,49].This technique explores the fractal self-similarity of fluctuations at different time scales in α x .The way in which these fluctuations scale is determined by the self-similarity parameter α, where fluctuations at time scales equal to t/a can be made similar to those at the larger time scale t simply by multiplying with the factor a α .
In order to estimate the self-similarity parameter α the time series is integrated in time and divided into nonoverlapping sections of length n.For each section the local trend is removed by subtracting a fitted polynomial.While the order of the polynomial is arbitrary, the convention is to use first-order fits [19,49].The root mean square fluctuation for the scale equal to n is then given by where Y (t) is the integrated and detrended time series and N is its length.The fluctuation F (n) provides a measure of the amplitude at the corresponding scale, which follows a scaling law if the time series is fractal.By plotting log F (n) against log n, the self-similarity parameter is then given by the gradient of the line.For completely uncorrelated white Gaussian noise (expected in chronotaxic systems) the parameter for α x returns a value of 0.5, while integrated white Gaussian noise (expected in nonchronotaxic systems) returns a value of 1.5.Hence, the self-similarity parameter for the fluctuations in α x provides a gauge of chronotaxicity.
To summarize, the procedure uses the difference α x between the estimated phase of the system and its attractor which are extracted from the time series.To identify the existence of phase slips, α x is averaged over different time scales τ and the change in the difference over these time scales |d α τ x | is calculated.For nonchronotaxic systems the Gaussian shape in the distribution of |d α τ x | remains the same for different values of τ .On the other hand, in chronotaxic systems containing phase slips the shape of the distribution is the same for small τ but shifts and obtains a new peak at longer time scales.Alternatively, in the absence of phase slips DFA can be applied to α x to obtain its fractal properties.The self-similarity parameter can then be used to determine whether or not the system is chronotaxic, where typically α < 1 for chronotaxic systems.

V. APPLICATIONS A. Phase oscillators
The above methodology is now applied to time series data of a chronotaxic and nonchronotaxic system.In this case the systems are phase oscillators, where the oscillatory time series x is generated from the phase α x using x = sin(α x ).
Consider the system αp = ω 0 (t), ( where ε is the coupling strength.The function ξ (t) is white Gaussian noise with standard deviation η = √ 2E, where ξ (t) = 0, ξ (t)ξ (τ ) = δ(t − τ )E.The frequency of α p is given as  19) with η = 0.8 using Bayesian-based inference.The parameter ε was varied so that the system switched between chronotaxic and nonchronotaxic behavior in time.In this particular example, the nonchronotaxic periods cause the system to continuously phase slip, which meant that αA * x was found from the Morlet wavelet transform with f 0 = 10 (high-frequency resolution), shown in (a).The frequency was then smoothed using a 100 s moving average, the result of which is shown as the dashed line.Conversely, the phase α * x was extracted using f 0 = 0.5 (high time resolution) which corresponds to frequency traced by the dashed line in (b).The coupling functions for the two phases were inferred, where (c) shows the directionality index from α A * x to α * x .In (d) the chronotaxic index based on the inferred coupling functions (black line) tracks the change between chronotaxic and nonchronotaxic dynamics due to the variation of the parameter ε (gray line).
with ω 1 = 0.002 rad s −1 ω 2 = 0.001π rad s −1 .The system is chronotaxic when |ε| 1, where for ω 0 (t) > 0 the attracting deterministic phase can be found analytically as [1] α where k is an arbitrary integer number.Now consider the following system which undergoes similar dynamics but is not chronotaxic: The parameters ω 0 (t) and ξ (t) are taken to be the same as before.While the system does not remain at a stable point in the phase, the oscillations retain a frequency around that defined by α p in the previous system because the natural frequency is given by the same function.
Figure 2 shows the analysis of the system given by ( 19) for the case where ε varies in time.After the estimated phases α * x and α A * x are extracted the Bayesian-based inference method is applied.The method correctly finds that for ε > 1 (i.e., when the system is chronotaxic) the inferred parameters result in synchronization and I chrono = 1, while for ε < 1 (i.e., when the system is nonchronotaxic) there is no synchronization and I chrono = 0.Meanwhile, the method finds a coupling in the direction from α A * x to α * x which implies that the synchronization is the result of a point attractor in the phase α * x .
It is also worth noting in Fig. 2(c) that D decreases or goes close to zero during the chronotaxic epochs.This effect occurs when the phases become strongly coherent, resulting in not enough available information for inference of the direction of coupling.In the case of moderate noise when phase slips exist, the directionality inference will be correct.If the noise becomes extremely large, introducing too many phase slips, the precision of the inference will be reduced.
The method based on DFA is demonstrated in Fig. 3, where the system in ( 22) is used for the nonchronotaxic example.This is a more difficult case for the Bayesian method because α x and α A x are close to being synchronized even with the effect of perturbations.In both the time and the time-frequency domains the two systems are also more or less indistinguishable.The difference can only be seen after calculating α x , which in this case was detrended below 200 s to remove the influence of the slow drift that resulted from imperfect phase extraction.The corresponding self-similarity parameter estimated from DFA shows the nonchronotaxic system to have a divergence almost the same as would be expected for integrated white Gaussian noise.On the other hand, the divergence in the chronotaxic system is much closer to the exponent for the perturbations.The value is not the same as for the noise because perturbations from the attractor do not decay instantaneously, meaning that some integration of the noise still occurs.It is also worth noting that even when the system is chronotaxic, if phase slips are present then the estimated self-similarity parameter is still close to that expected for a random walk.However, this is actually expected since the stability in the phase does not extend beyond a single cycle, which means the fluctuations from noise-induced phase slips are still integrated over.19) with ε = 1.7 and the non-chronotaxic system (22).The systems were perturbed with white Gaussian noise with strength η = 0.7 that was small enough to not cause phase slips in the chronotaxic system.The time series sin(α x ), with the first 20 s shown in (a), were analysed using the Morlet wavelet transform where |W T (s,t)| is shown in (b).An estimate of the phase α x was extracted from the wavelet transform with f 0 = 1 while an estimate of α A x was obtained by smoothing the frequency of the oscillation with a 100 s moving average and integrating over time.The estimated direction of coupling from α A * x to α * x calculated using the Bayesian-based method is shown in (c).The difference α x between the two phases was calculated and then detrended with a 200 s moving average, with the resulting time series shown in (d).DFA was applied to this time series as shown in (e), where the solid line gives the scaling of the fluctuations and the dotted line is a least-squares linear fit.In the case of the chronotaxic system the gray line corresponds to η = 1.4 where phase slips are present.The number shown next to the line is the gradient of the fit, which provides an estimate of the self-similarity parameter and gives a gauge of the chronotaxicity.
In the last example shown in Fig. 4 the noise is stronger, resulting in phase slips in the chronotaxic system.This strong noise also leads to large perturbations in the nonchronotaxic system, which causes the estimated direction of coupling to be spurious.However, the difference between phase slips and the large perturbations in the nonchronotaxic systems is difficult to quantify in the time and time-frequency domains.This difference is apparent in α x and becomes obvious after calculating the distribution of the perturbations for various time scales.While the shape of the distribution in the nonchronotaxic system remains the same at all time scales, the distributions of the chronotaxic system are the same shape for small values of τ but become shifted at longer time scales where the effect of the phase slips is included.The method is therefore able to detect the discontinuity in the effect of the perturbations which results from the presence of the point attractor in chronotaxic systems.

B. Cardiac system
Candidates for chronotaxic systems in the real world include both living and nonliving systems.In biological systems, the importance of "inner variables" analogous to p in (1) has been highlighted in Ref. [50].Similar driven steady states are also studied in quantum critical systems and solid state physics [51,52].
To illustrate how the approach presented here can influence physical models, the above methodology is now applied to data from one of these real systems.This comes from the experiment described in Ref. [1], where the breathing rate of a healthy young subject was paced quasiperiodically while the cardiac function was monitored using an electrocardiogram (ECG).x .The direction of coupling from α A * x to α * x calculated using the Bayesian-based method is shown in (c).The difference α x , shown in (d), was detrended using a 200 s moving average.The DFA of α x is shown in (e), where the number next to the line is the DFA exponent as estimated using a linear fit (dashed line).In (f) the estimated probability distributions of the perturbations are shown, where the numbers in seconds correspond to the time scales τ .
The phase dynamics of the cardiac function is comprised of deterministic interactions with the other oscillatory systems involved in cardiovascular regulation.The oscillatory processes that modulate cardiac activity are known or hypothesized to correspond to the breathing (0.145-0.6 Hz), smooth muscle (0.052-0.145Hz), sympathetic nerves (0.021-0.052Hz) and the vascular endothelium (0.005-0.021Hz) [41,53].Since this range of oscillations extends close to the oscillating frequency of the heart itself (0.6-2 Hz), the observed frequency must be averaged over short time scales in order to obtain α A * x .This essentially means that α * x and α A * x share very similar dynamics so that over short time windows the effect of the perturbations in α * x is not statistically significant and the two phases appear synchronized even when they are not.Consequently, the phase fluctuation methods are more appropriate than the Bayesianbased method in this situation.
Figures 5 and 6 show the results from the analysis of the data.In the case of the phase extracted from the ECG signal, Fig. 5, both methods indicate that the phase is not stable.The self-similarity parameter suggests that the divergence from the attracting phase is very close to a random walk.The distribution of the perturbations also retains the same shape for all time scales, implying that phase slips do not occur.
The analysis for the phase of the heart rate variability (HRV) is shown in Fig. 6.The most noticeable oscillation in the wavelet transform, but not the only one in the HRV, is known from physiology as respiratory sinus arrhythmia [54,55].Other studies have shown that this results from the heart being driven by the respiratory oscillation [15,19,56], which in this case had a frequency that varied in the range 0.06-0.25 Hz as it was paced.This respiration variability was obtained directly from the HRV during the extraction of αA * x .The analysis of the extracted phases gave a slightly smaller estimate of the self-similarity parameter but it is still close to the value expected for either a nonchronotaxic system or a chronotaxic system with many phase slips.However, the distribution of perturbations can be seen to flatten and shift away from 0 in a similar fashion to the distributions in Fig. 4. The nature of the fluctuations which are likely to be responsible for this shift is shown in the Supplemental Material [57], where the HRV oscillation is viewed relative to the respiration.Previously, these effects were accounted for by the fact that the perturbations to the system are inhomogeneous.However, chronotaxic systems now provide an alternative hypothesis, which is that a coupling from the respiration to the HRV results in an oscillation with a stable phase.The directionality index in Fig. 6 shows that such a coupling exists from α A * x to α x .This suggests that the respiration takes the form of p in (1).
In describing the chronotaxic properties of the healthy heart under quasiperiodically paced breathing, these results motivate the following addition to the model presented in Ref. [20].In this paper the coupling functions between the heart and respiration were extracted using Bayesian-based inference to produce a physical model of the interaction.In the model, the couplings were divided between the self-interaction s h of the heart's phase on its own dynamics, the direct interaction d h with the phase of the respiration, and the indirect interaction i h , which includes coupling terms dependent on both phases.The full model is presented in the paper as φh = 2πω h (t) + s h (φ h ,t) + d h (φ r ,t) where φ h is the phase of the heart, ω h (t) is its natural frequency, φ r is the phase of the respiration, and ξ h (t) is x .The direction of coupling from α A * x to α * x calculated using the Bayesian-based method is shown in (c).The difference α x , shown in (d), was detrended using a 200 s moving average.The DFA α x is shown in (e), where the number next to the line DFA exponent as estimated using a linear fit, which suggests that the system is either nonchronotaxic or is chronotaxic with phase slips.In (f) the estimated probability distributions of the perturbations are shown, where the numbers in seconds correspond to the time scales τ .The change in the shape of the distribution is consistent with a chronotaxic system with phase slips.assumed to be white Gaussian noise.The oscillation related to respiratory sinus arrhythmia originates primarily from the direct interaction, which is modeled as (1)  h (t) sin φ r + ϕ (1)  h (t) + (2)  h (t) sin 2φ r + ϕ (2)  h (t) , (24) where the parameters (1)  h (t), (2)  h (t), ϕ (1)  h (t), and ϕ (2)  h (t) allow the model to be time varying.However, in this form the subsequent oscillation in the HRV cannot be chronotaxic.This instead requires an intermediate phase to be present in the coupling, with the updated model taking the form φx = f x (φ x ,φ r ), d h = (1)  h (t) sin φ x + ϕ (1)  h (t) + (2)  h (t) sin 2φ x + ϕ (2)  h (t) , (25) where φ x is the phase of the oscillation observed in the HRV and f x (φ x ,φ r ) is a coupling function between this phase and the respiratory phase that results in a point attractor φ A x .The function f x (φ x ,φ r ) is likely to be dependent on time as well, which would mean the strength of the coupling from φ A x to φ x , and consequently the resistance to perturbations, is able to vary.
The model shown in (23) is simplified as the other oscillations present in the HRV are mainly accounted for by the time-varying parameters.In earlier studies, these oscillations are included in either the noise perturbations [58] or explicitly modeled as separate oscillators [59].Here we now go further by opening up the possibility that these are also chronotaxic oscillations with their own attracting phases.The nature of interaction of these oscillations through the HRV is in fact analogous to the parametric coupling model of the cardiovascular system proposed in Ref. [60].However, it should also be noted that the model presented here is based on data from a single experiment involving paced respiration, which means that further experiments are needed to confirm the applicability of these modifications in general.

VI. SUMMARY AND CONCLUSION
Chronotaxic systems have opened up a brand new area of time-dependent dynamical systems.The most important contribution of this subset of nonautonomous dynamical systems is that they can explain why oscillatory systems, such as the function of a healthy heart, can have a time-varying frequency and generate complex dynamics.However, detecting this stability using the inverse approach is a nontrivial task and requires analysis beyond time-dependent representations used in general for nonautonomous systems.
The framework presented here is able to distinguish chronotaxic and nonchronotaxic systems which appear similar in the main domains of analysis, including the time-frequency domain.The emphasis on single-variable time series means that these methods can be applied to almost any data from an oscillatory system to test for chronotaxicity.We apply a recently developed Bayesian-based approach because this method is able to accurately track nonautonomous dynamics when the system is perturbed by noise, as is the case in chronotaxic systems.In particular, transitions between chronotaxic and non chronotaxic behavior can be detected in the same time series.A second set of methods, named phase fluctuation analysis, uses another approach to measure the statistics of the observed perturbations to the system and can detect chronotaxic dynamics even when the analytical phases are close to synchronization.
The current framework should therefore find applicability to real data from many different sources.Moreover, the ability to distinguish between chronotaxic and nonchronotaxic systems provides a measure of when a system switches between corresponding physical states.In the case of the chronotaxic oscillation measured in the HRV, this could be used in either the diagnosis or prediction of disease.While the work here establishes this possibility, other methods may also be developed in the future to track chronotaxic dynamics and, hence, distinguish between physical states across a wider range of systems.

FIG. 1 .
FIG. 1. (Color online) Extraction of α A x from a time series.(a) An oscillation with a varying frequency and noise-induced phase slips is shown in the time domain.(b) Morlet wavelet transform of the time series where the colors indicate |W T (s,t)|.(c) The estimated angular velocity αA *x found by ridge extraction in the wavelet transform and then smoothing using a 2 s moving average (solid line) and the actual value of αAx from the analytic solution (dashed line).(d) The phase α A * x found by integrating over αA * x (solid line) and the actual dynamics of α A x from the analytic solution (dashed line).The offset of the phases is due to the estimated phase starting at a later time as αAx cannot be defined for the interval where the wavelets run over the ends of the time series.

FIG. 2 .
FIG. 2. (Color online) Analysis of the system shown in(19) with η = 0.8 using Bayesian-based inference.The parameter ε was varied so that the system switched between chronotaxic and nonchronotaxic behavior in time.In this particular example, the nonchronotaxic periods cause the system to continuously phase slip, which meant that αA *x was found from the Morlet wavelet transform with f 0 = 10 (high-frequency resolution), shown in (a).The frequency was then smoothed using a 100 s moving average, the result of which is shown as the dashed line.Conversely, the phase α *x was extracted using f 0 = 0.5 (high time resolution) which corresponds to frequency traced by the dashed line in (b).The coupling functions for the two phases were inferred, where (c) shows the directionality index from α A * x to α * x .In (d) the chronotaxic index based on the inferred coupling functions (black line) tracks the change between chronotaxic and nonchronotaxic dynamics due to the variation of the parameter ε (gray line).

FIG. 3 .
FIG.3.(Color online) Application of the phase fluctuation analysis methods to the chronotaxic system(19) with ε = 1.7 and the non-chronotaxic system(22).The systems were perturbed with white Gaussian noise with strength η = 0.7 that was small enough to not cause phase slips in the chronotaxic system.The time series sin(α x ), with the first 20 s shown in (a), were analysed using the Morlet wavelet transform where |W T (s,t)| is shown in (b).An estimate of the phase α x was extracted from the wavelet transform with f 0 = 1 while an estimate of α Ax was obtained by smoothing the frequency of the oscillation with a 100 s moving average and integrating over time.The estimated direction of coupling from α A * x to α * x calculated using the Bayesian-based method is shown in (c).The difference α x between the two phases was calculated and then detrended with a 200 s moving average, with the resulting time series shown in (d).DFA was applied to this time series as shown in (e), where the solid line gives the scaling of the fluctuations and the dotted line is a least-squares linear fit.In the case of the chronotaxic system the gray line corresponds to η = 1.4 where phase slips are present.The number shown next to the line is the gradient of the fit, which provides an estimate of the self-similarity parameter and gives a gauge of the chronotaxicity.

FIG. 4 .
FIG. 4. (Color online) Application of the phase fluctuation analysis methods to the chronotaxic system(19) with ε = 1.7 and the nonchronotaxic system(22).The systems were perturbed with white Gaussian noise with strength η = 1.4 that was large enough to cause phase slips in the chronotaxic system.The time series sin(α x ), with the first 20 s shown in (a), were analyzed using the Morlet wavelet transform where |W T (s,t)| is shown in (b).An estimate of the phase α x was extracted from the wavelet transform with f 0 = 1, while an estimate of α Ax was obtained by smoothing the frequency of the oscillation with a 100 s moving average and integrating over time.The estimated direction of coupling from α A * x to α * x calculated using the Bayesian-based method is shown in (c).The difference α x between these two phases was calculated, which is shown after detrending with a 200 s moving average in (d) for comparison.The distributions of |d α τx | are shown in (e) where the range shown was divided into 100 bins and the numbers in seconds indicate the value of τ .

FIG. 5 .
FIG. 5. (Color online) Analysis of the cardiac dynamics using phase fluctuation analysis.The oscillatory activity was extracted from an ECG, with a 20 s sample shown in (a).The phase and frequency of the dominant component in the ECG was extracted using the Morlet wavelet transform (b) with f 0 = 1.The frequency was smoothed using a 1 s moving average and integrated in time to find α A *x .The direction of coupling from α A * x to α * x calculated using the Bayesian-based method is shown in (c).The difference α x , shown in (d), was detrended using a 200 s moving average.The DFA of α x is shown in (e), where the number next to the line is the DFA exponent as estimated using a linear fit (dashed line).In (f) the estimated probability distributions of the perturbations are shown, where the numbers in seconds correspond to the time scales τ .

FIG. 6 .
FIG.6.(Color online) Analysis of theHRV using phase fluctuation analysis.The HRV was extracted from the Morlet wavelet transform of the ECG, with a 100 s sample shown in (a).The phase and frequency of the dominant component in the HRV was extracted using the wavelet transform (b) with f 0 = 0.5.The frequency was smoothed using a 6 s moving average and integrated in time to find α A *x .The direction of coupling from α A * x to α * x calculated using the Bayesian-based method is shown in (c).The difference α x , shown in (d), was detrended using a 200 s moving average.The DFA α x is shown in (e), where the number next to the line DFA exponent as estimated using a linear fit, which suggests that the system is either nonchronotaxic or is chronotaxic with phase slips.In (f) the estimated probability distributions of the perturbations are shown, where the numbers in seconds correspond to the time scales τ .The change in the shape of the distribution is consistent with a chronotaxic system with phase slips.