Detecting chronotaxic systems from single-variable time series with separable amplitude and phase

The recent introduction of chronotaxic systems provides the means to describe 1 nonautonomous systems with stable yet time-varying frequencies which are resistant to 2 continuous external perturbations. This approach facilitates realistic characterization of the 3 oscillations observed in living systems, including the observation of transitions in dynamics 4 which were not considered previously. The novelty of this approach necessitated the 5 development of a new set of methods for the inference of the dynamics and interactions 6 present in chronotaxic systems. These methods, based on Bayesian inference and detrended 7 fluctuation analysis, can identify chronotaxicity in phase dynamics extracted from a single 8 time series. Here, they are applied to numerical examples and real experimental EEG data. 9 We also review the current methods, including their assumptions and limitations, elaborate 10 on their implementation, and discuss future perspectives. 11


Introduction
The theory of nonautonomous dynamical systems has increasingly been recognised as a necessity in the treatment of the inherent time-variability of biological systems [1].Closer inspection of the dynamics observed in nature suggests that previous approaches to the characterization of temporal fluctuations in these observations may be insufficient.At first glance, biological fluctuations may appear random, leading to their description by stochastic models [2].The complexity observed in biological systems has also led to attempts to treat them with chaos theory [3], however this does not allow for the apparent stability of these systems, irrespective of their initial conditions.Such characteristics of biological oscillators suggests underlying determinism or control of both their amplitudes and frequencies, even with continuous perturbations.This phenomenon of biological systems resisting a natural tendency to disorder was discussed in terms of free energy minimization [4] and separation of internal and external states, but this approach is still based on random dynamics.A closely related, yet more natural, approach is to consider them as nonautonomous systems, which are explicitly time-dependent.Approaches based on reformulation of nonautonomous systems as higher dimensional autonomous systems introduce unnecessary complexity, whilst failing to accurately describe dynamics arising from nature due to the fact that these are open systems, subject to continuous variable external perturbations.Many living systems may be considered as nonautonomous oscillatory systems, with such time-varying dynamics observed in individual mitochondria [5], the cardio-respiratory system [6,7], the brain [8], and blood flow [9].
Although stability of the amplitude dynamics of an oscillator can be achieved with autonomous self-sustained limit cycle oscillators, the frequency of this oscillation could be easily changed by weak external perturbations [10].To account for a case where this frequency of oscillation is also robust to perturbations, yet time-dependent, a completely new approach is required.Thus, nonautonomous systems with stable, yet time-varying frequencies were recently addressed, and formulated as chronotaxic systems [10][11][12].Chronotaxic systems possess a time-dependent point attractor provided by an external drive system.This allows the frequency of oscillations to be prescribed externally through this driver and response system, giving rise to determinism, even when faced with strong perturbations.
Once these properties of the underlying system have been recognised, the next problem is how to infer these dynamics and interactions from direct observations, i.e., via the inverse approach.In a chronotaxic system, particularly one found in nature, whilst the underlying dynamics are defined by the external driver, the system will likely still be affected by other influences and noise.These may mask the chronotaxic dynamics if the correct analytical approach is not applied.For example, the inherent time-variability of the frequency of the dynamics arising from a chronotaxic system means that it cannot be accurately characterized by any method based on averaging.This novel class of systems requires new inverse approach methods, with the focus on the extraction and identification of the dynamics of the drive system, and its influence on the response system.Here, we review the current state of inverse approach methods for the identification of chronotaxicity from a single time series of the response system in which the phase and amplitude dynamics are separable.We then apply these methods to numerically simulated and real experimental data.Section 2 presents the mathematical formulation of chronotaxic systems, Section 3 describes current inverse approach methods and their application to the detection of chronotaxicity.In Section 4, numerical examples are presented to demonstrate the methods, and their assumptions and limitations are discussed.The inverse approach methods are also applied to real experimental data.Finally, Section 5 discusses future directions of the work.

Chronotaxic Systems
The crucial concept in the theory of chronotaxic systems is the ability to resist continuous external perturbations.In autonomous systems, such an ability is provided by a stationary steady state, allowing the system to always return to the vicinity of this steady state when continuously externally perturbed.However, only in nonautonomous (thermodynamically open) systems can the position of this steady state change in time, i.e., be outside of equilibrium.In such a case, not only the stationary state of a system, but also its time-dependent dynamics will be able to resist continuous external perturbations.These oscillatory nonautonomous dynamical systems with time-dependent steady states were introduced by Suprunenko et al. [10] and named chronotaxic systems, emphasizing that their dynamics is ordered in time (chronos -time, taxis -order).
Mathematically, nonautonomous dynamical systems and, consequently, chronotaxic systems, are defined by the following system of equations: where p ∈ R n , x ∈ R m , f : R n → R n , g : R m × R n → R m , in which n and m can be any positive integers.Importantly, the solution x(t, t 0 , x 0 ) of Equation (1) depends on the actual time t as well as on the initial conditions (t 0 , x 0 ), whereas the solution p(t − t 0 , p 0 ) depends only on initial condition p 0 and on the time of evolution t − t 0 .The subsystem x is nonautonomous in the sense that it can be described by an equation which depends on time explicitly, e.g., ẋ = g(x, p(t)).A chronotaxic system is described by x which is assumed to be observable, and p which may be inaccessible for observation, as often occurs when studying real systems.Rather than assuming or approximating the dynamics of p, we focus on the dynamics of x and use only the following simple assumption: system p is assumed to be such that it creates a time-dependent steady state in the dynamics of x, which is schematically shown in Figure 1a.Therefore, the whole external environment with respect to x is divided into two parts.The first part is given by p which is only that part which makes the system x chronotaxic (defined below), i.e., an unperturbed chronotaxic system.The second part contains the rest of the environment, and is therefore considered as external perturbations.
Firstly, we provide a mathematical formulation of unperturbed chronotaxic systems.The defining component of an unperturbed chronotaxic system is a time-dependent steady state, also called a point attractor, and denoted x A (t).
Usually, a steady state is defined using a so-called forward limit, when forward time approaches infinity.Assuming that the whole space R m of x is a basin of attraction, i.e., that for any initial condition x 0 at time t 0 the solution of the system asymptotically approaches the time-dependent steady state x A , a condition of forward attraction for x A is the following, ( This condition can only be satisfied when the chronotaxic system is not perturbed.However, taking into account the time-dependence of x A (t), this condition is not satisfactory in terms of defining the time-dependent point attractor.Any solution x(t, t 0 , x 0 ) which satisfies Equation (2) with given x A (t) can also be considered as a time-dependent point attractor.Moreover, when dealing with living systems, it is crucially important to describe stability at the current time t and not in the infinite future.This problem is resolved by employing a condition of pullback attraction, which should also be satisfied by x A (t) in a chronotaxic system, One can see, that this condition defines a time-dependent point attractor at current time t.
Considering the condition in Equation ( 3) at all times t > −∞, it follows that the time-dependent point attractor should also satisfy the invariance condition, i.e., the condition that x A is a solution of the system Equation (1), Equations ( 2) and (3) determine asymptotic convergence in the infinite future or starting from the infinite past.Asymptotic convergence allows the dynamics of x(t, t 0 , x 0 ) to deviate from x A during a certain finite time interval.Thus, during this time-interval the ability to resist continuous external perturbations will be absent.Therefore, in order to characterize the ability of living systems to sustain their time-dependent dynamics at finite time intervals, a chronotaxic system should satisfy the condition of contraction, or equivalently the attraction at all times.This means that in the phase space R m , x ∈ R m , there should be a contraction region C(t) such that for any two trajectories x 1 , x 2 of a system inside the contraction region x i (t, t 0 , x 0i ) ∈ C(t), i = 1, 2, the distance between them can only decrease, However, in general the contraction region C(t) can be finite, and different trajectories can eventually leave this region.Therefore, in a chronotaxic system the contraction region should contain a finite area A , A ⊂ C, such that solutions of the system starting in A never leave it, In such a case, fulfillment of these conditions guarantees that the time-dependent point attractor x A is located inside the area A inside the contraction region C.
Alternatively, the trajectory x A (t) can be viewed as a linearly attracting uniformly hyperbolic trajectory [13], so that the distance between a neighboring trajectory and x A (t) can only decrease in an unperturbed chronotaxic system.For more details and for relations between chronotaxic and other dynamical systems see Reference [12].A simple example of an unperturbed chronotaxic system is given by unidirectionally coupled phase oscillators with unwrapped phase ϕ x ∈ (−∞, ∞) driven by a phase where φp (t) = ω(t).
As an example, for a particular choice ε(t) = ω(t) > 0 and ω 0 (t) = 0, the equation can be integrated, and the limit t 0 → −∞ can be calculated, leading to the explicit expression for a time-dependent point attractor of an unperturbed chronotaxic system, ϕ A x (t) = ϕ p (t) − π/2 + 2πk, and k is any integer.It should be noted that the dynamics of p can be complex, stochastic, or chaotic, provided that the above conditions are met.Nevertheless, the dynamics of x will be determined by the dynamics of p, and therefore it will be deterministic, at least in unperturbed chronotaxic systems.When considering perturbed chronotaxic systems, for simplicity it is sufficient to consider perturbations only to the x component, as any perturbations to p can be included in its dynamics assuming that x does not influence p.In perturbed chronotaxic systems, which model real life systems, the general external perturbations will create complex dynamics of x with a stochastic component.Such dynamics may look very complex: perturbations can push trajectories away from the contraction region, therefore they can temporarily deviate before they converge.Despite this, due to the existence of the contraction region, the system will resist continuous external perturbations.The time-dependent dynamics of a perturbed chronotaxic system will be very close to the dynamics of the unperturbed chronotaxic system provided the perturbations are weak enough.In the case of very strong continuous perturbations, such perturbations may override the driving system p, and become effectively a new driver, causing the initial point attractor of the chronotaxic system to disappear.However, it may be restored once the perturbations again become sufficiently weak.
Thus, when perturbations do not destroy the chronotaxic properties of a system, the stable deterministic component of its dynamics can be identified, as will be shown below.This reduces the complexity of the system, allowing us to filter out the stochastic component and focus on deterministic dynamics and interactions between system x and its driver p, [10,14].For such complex systems as living systems, it has the potential to extract properties of the system which were previously neglected.

Inverse Approaches to Nonautonomous Dynamical Systems
A wide range of observed properties in living systems can be explained by considering them as nonautonomous.Despite this, difficulties in their analysis as such have led to many unsuccessful attempts to apply methods more suited to autonomous systems.From a single time series arising from a dynamical system, inverse approach methods can be used to infer the underlying dynamics of this system, in terms of phase or amplitude.In deterministic systems, phase space analysis is usually the first point of call, i.e., reconstruction of the attractor in phase space.This can be achieved with only a single time series using embedding, in which the dimensions of the reconstructed attractor are composed of time-delayed versions of the data in the time series [15].This approach works well for autonomous systems, but does not consider the possibility of time-dependent attractors [1].Phase space methods are particularly suited to the treatment of the dynamics observed in chaotic systems.In contrast, nonautonomous systems appear very complex in phase space.To incorporate time-dependence into these systems, extra dimensions in phase space are required, introducing unnecessary complexity to the problem.
In the case that there is only one oscillation in the time series, the Hilbert transform can be applied to obtain the complex analytic signal, from which the instantaneous phase can be determined directly.If the time series contains more than one component of interest, for example different oscillatory modes, it can be decomposed into its constituent parts using a method such as empirical mode decomposition (EMD), in which peak/trough detection is used to create upper and lower envelopes.From these, a trend is defined and subtracted from the signal to produce a series of intrinsic mode functions, each one representing an oscillatory component of the time series [15].However, there are some limitations when applying EMD to nonstationary data.
Many signal analysis methods assume stationarity of the frequency distribution of the data, but in nonautonomous systems this assumption is not valid.Single variable time series, particularly those from living systems, must be treated as arising from nonautonomous dynamical systems, due to time-dependent influences of variables other than the one under study.Approaches based on windowing have been applied in order to attempt to treat time variability in data, but these potentially lose crucial information.For example, in phase space reconstruction the window may not be of a sufficient size to capture the whole of the attractor, or its variations in time.Application of the Fourier transform to nonstationary data will result in a blurred or misleading power spectra, severely limiting its usefulness.The windowing approach has been applied here with some success, in the form of the Short Time Fourier transform (STFT), but the use of windowing leads to limitations; the better the time resolution, the worse the frequency resolution (known as the Gabor limit [16]).In addition, the fixed time-frequency relationship at all scales in a windowed Fourier transform severely limits its usefulness for the analysis of low frequency oscillations.This problem can be addressed by using the continuous wavelet transform (CWT), which provides a logarithmic frequency scale (see Section 3.3).The CWT is based on wavelets rather than the sines and cosines used in Fourier based methods.The simultaneous observation of the time and frequency domains is extremely useful in the visualization of dynamical systems and their time evolution.As a result, development of wavelet-based methods specifically for the treatment of time-dependent dynamics is now a very active field of research [15], including wavelet phase coherence [17], the synchrosqueezed transform [18,19] and wavelet bispectrum.
In addition to determining the characteristics of the underlying dynamics of single nonautonomous oscillatory modes, inverse approach methods are also used to decompose their interactions.One of the most well known characteristics of interacting systems is synchronization between oscillatory components, i.e., a fixed relationship between their phases or amplitudes.Once the phases of oscillations have been extracted from a time series, a measure of phase synchronization can be calculated using synchronization indices or phase coherence [15].However, these methods do not account for time-varying synchronization.Dynamical Bayesian inference is able to detect time-varying synchronization in a system, whilst simultaneously inferring the direction of coupling and time-evolving coupling functions [20,21].In the time-frequency domain, wavelet phase coherence can be used to monitor phase relationships over time and frequency by utilising the phase information obtained from the continuous wavelet transform [15,17].In a similar way, couplings between oscillators can be detected and quantified using wavelet bispectrum [22].The ability of these methods to directly take into account the time-varying characteristics of data makes them ideal for application to nonautonomous systems.

Detecting Chronotaxicity
Here we present two distinct inverse approach methods which may be utilised in the detection of chronotaxicity: phase fluctuation analysis (PFA) and dynamical Bayesian inference.It should be noted that the current methods are only applicable to phase dynamics in the context of the detection of chronotaxicity, i.e., we focus on the ability of the time-varying frequency to resist continuous external perturbations.The two methods rely on a different inferring basis.Phase fluctuation analysis provides a measure of statistical effects observed in a signal, whilst the dynamical Bayesian inference method infers a model of differential equations and gives a measure of dynamical mechanisms, i.e., the evaluation of chronotaxicity relies on the inferred parameters of the model.PFA is said to infer a functional connectivity, while the dynamical Bayesian inference method infers effective connectivity [23].The optimal method to use depends on the characteristics of the data, as detailed below.
It is possible to detect whether a system is chronotaxic or not by observing the distribution of the fluctuations in the system relative to its unperturbed trajectory.This comes from the fact that if the original distribution of the perturbations is known, then the stability of the system relative to the unperturbed trajectory (which by definition follows the time-dependent point attractor in a chronotaxic system) can be determined from how these perturbations grow or decay over time.For example, take the non-chronotaxic phase oscillator [24] φx = ω 0 (t) + η(t), where ω 0 (t) is the time-dependent natural frequency and the observed phase ϕ x is perturbed by noise fluctuations η(t).Integrating we find, Assuming that ω 0 (t) > 0 and η(t) is an uncorrelated Gaussian process, this means that the dynamics of ϕ x will consist of a monotonically increasing phase perturbed by a random walk noise (Brownian motion).However, the situation is different for a chronotaxic phase oscillator, e.g., where ϕ p is an external phase and |ε| > 1.In this case the stability provided by the point attractor causes each noise perturbation to decay over time, preventing η(t) from being integrated over to the same extent.The perturbations still do not decay instantly, as the system takes time to return to the point attractor, meaning that some integration of the noise still takes place.However, the size of the observed perturbations over longer time scales is greatly reduced, causing a change in the overall distribution from that expected for Brownian motion.

Extracting the Perturbed and Unperturbed Phases
The first problem of generating a method based on the above principle is how to extract both the perturbed and unperturbed phases of the system from an observed time series.This usually requires the separation of the amplitude and phase for an oscillation in the time series, which is possible using time-frequency domain analysis [15].The analytic signal generated by the Hilbert transform can also be used, but this requires corrections for nonlinear oscillations and cannot be used when more than one oscillation is present in the time series [14,25,26].Additionally, the use of the Hilbert transform requires the use of protophase-to-phase conversion [26].
A time-frequency representation with an optimal frequency resolution of the time series f (t) of length L is provided by the continuous wavelet transform [27], 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 Morlet wavelet is a commonly used mother wavelet and is defined as [28], where the corresponding frequency is given by 1/s and f 0 is a parameter known as central frequency which defines the time/frequency resolution [27].
Oscillations can be traced in W T (s, t) using either a ridge-extraction method [29,30] or the synchrosqueezed wavelet transform (SWT) [18].These extraction methods can be 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 intra-cycle dynamics.The phase ϕ x 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.
With the estimated perturbed phase ϕ * x extracted, further work is needed to obtain the unperturbed phase ϕ A * x .In particular, it is difficult to separate the dynamics corresponding to ϕ A x from the effect of the noise perturbations η(t).This task is simplified by assuming that the dynamics of ϕ A x are confined to time scales larger than a single cycle and that the noise is either weak or comparable in magnitude to these dynamics.
With these assumptions, an estimate of ϕ A x can be found by filtering out high frequency components of ϕ * x .However, such a filter should not smooth over the dynamics of ϕ A x .An optimal way of removing these high frequency noise fluctuations without affecting the unperturbed dynamics is to instead smooth over the frequency extracted from the wavelet transform [15].This provides the estimated angular velocity φA x , which can in turn be integrated over time to give the estimated phase of the driver ϕ A * x .For further methodological details on phase extraction see Section 4.2.

Dynamical Bayesian Inference
One approach to the detection of chronotaxicity is the application of dynamical Bayesian inference to the extracted perturbed (ϕ * x ) and unperturbed (ϕ A * x ) phase estimates in order to model their interactions.
In dynamical systems Bayesian inference can simultaneously detect time-varying synchronization, directionality of coupling and time-evolving coupling functions [20,21].The characteristics of these coupling functions between ϕ * x and ϕ A * x may reveal the dynamic mechanisms of the system in terms of chronotaxicity.Bayesian inference is able to track time-dependent system parameters, meaning that it is particularly useful for the detection of chronotaxicity in systems which move in and out of a chronotaxic state.
Following extraction of phases from the continuous wavelet transform, as described in Section 3.3, we assume their dynamics is described by [14,20,21] 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 Based on the periodic nature of the system, the basis functions are modeled using the Fourier bases and where k, r, s = 0.In practice, it is reasonable to assume that the dynamics will be well described by a finite number of Fourier terms, denoted A i,k (ϕ i , ϕ j ).The corresponding parameters from ãi and bi then form the parameter vector c (i) k .The inference of these parameters utilises 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 negative log-likelihood function is with implicit summation over repeated indices k, l, i, j.The log-likelihood is a function of the Fourier coefficients of the phases [20].Assuming a multivariate normal distribution as the prior for parameters c k with means c and covariances Σ = Ξ −1 , the stationary point of S can be calculated recursively from These are calculated within a moving time window, with the current prior depending on information from the posterior of the previous window.The inferred parameters of the basis functions can be used to determine whether synchronization results.The presence of synchronization provides evidence that the system is chronotaxic, however it remains unclear from which coupling function the stability arises without calculating the direction of coupling [31], where are the Euclidean norms of the parameters.The odd parameters correspond to the coupling terms inferred for ϕ 1 in the direction 2 → 1, and the even parameters correspond to the coupling terms inferred for ϕ 2 in the direction 1 → 2. See [32] for further details and an in depth tutorial on dynamical Bayesian inference and its implementation.In summary, Bayesian inference is applied to ϕ * x and ϕ A * x , following their extraction from the time series (see Section 3.3).The time evolution of the coupling parameters for each phase is inferred and these are used to determine the synchronization state of the system, and the direction of coupling between the phases.In a chronotaxic system we require the driver and response systems to be almost or fully synchronized, and also that the direction of coupling is only from the driver ϕ A x to ϕ x .The basis of this method is the calculation of the synchronization and direction of coupling of the system in order to determine chronotaxicity.However, the more synchronized the driver is with the response system, the less information flow between the two.With less information from which to infer parameters, most directionality methods, including Bayesian inference, become less reliable, and whilst synchronization may still be accurately detected, the direction of coupling will become less accurate the closer the system gets to synchronization.With frequent external perturbations, intermittent transitions, and moderate dynamical noise, there is greater information flow, and thus the inference is more precise, but this cannot be assumed in chronotaxic systems.In real systems, the synchronization state is not known beforehand, thus a more robust method is required, which can identify chronotaxicity even in systems close to synchronization.

Phase Fluctuation Analysis
Phase fluctuation analysis is effective even when ϕ x and ϕ p are almost synchronized.Given the estimates of ϕ x and ϕ A x , the next step is to analyse x to find the distribution of fluctuations in the system relative to the unperturbed trajectory.
In order to quantify the distribution of fluctuations, detrended fluctuation analysis (DFA) can be performed on ∆ϕ x [6,33].Following from the observations of Section 3.2, this method tries to estimate the fractal self-similarity of fluctuations at different time scales in order to distinguish the random walk fluctuations of non-chronotaxic systems from the less-integrated fluctuations of chronotaxic systems.The scaling of these fluctuations 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 by multiplying with the factor a α .The DFA exponent of (f) incorrectly suggests the system is non-chronotaxic.To distinguish between a non-chronotaxic system and a chronotaxic system with phase slips, the delayed distributions were calculated (see Section 3.5) in the non-chronotaxic (g) and chronotaxic (h) case.
In order to calculate α, the time series ∆ϕ x is integrated in time and divided into sections of length n.For each section the local trend is removed by subtracting a fitted polynomial-usually a first order linear fit [6,33].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 amplitude F (n) follows a scaling law if the time series is fractal.By plotting log F (n) against log n, the value of α is simply the gradient of the line.For completely uncorrelated white Gaussian noise (the noise assumed to perturb the system) the parameter for α has a value of 0.5, while integrated white Gaussian noise (expected in non-chronotaxic systems) returns a value of 1.5.Note that this assumes that the noise does not cause phase slips in ϕ x .This would cause perturbations over large time scales (i.e., greater than one cycle) to not decay even if the system was chronotaxic.In these cases another approach should be used instead [14].If there are large perturbations which cause the system to move far enough forward or behind the current cycle to be attracted instead by an adjacent cycle, known as a phase slip, this will result in an increased DFA exponent.This can result from large jumps in the extracted phase fluctuations.To distinguish between the case of a chronotaxic system with phase slips, and a non-chronotaxic system, we consider the fact that in the latter, perturbations may cause ∆ϕ x to change by 2π, but these are part of a continuous probability distribution, in contrast to the chronotaxic case.Phase slips can be detected by calculating the distribution of the difference between the phase fluctuations ∆ϕ x (t) and these fluctuations delayed by a time scale τ .d∆ϕ τ x (t) = ∆ϕ x (t + τ ) − ∆ϕ x (t) therefore gives information about the perturbations of the system over that time scale.When phase slips are present, the distribution of |d∆ϕ τ x | changes with respect to τ [14].An example of this difference is shown in Figure 2g,h, and can also be seen in real biological systems, as previously demonstrated in the heart rate variability [14].

Numerical Simulations
The basis of the PFA method is the quantification of the fundamental difference between phase fluctuation distributions in oscillatory systems, depending on their chronotaxicity.Here, we illustrate this characteristic using the simplest realisation of a chronotaxic system, two unidirectionally coupled oscillators (see Figure 1b): where ϕ p and ϕ x are the instantaneous phases of the driving and the driven oscillators, respectively, ω p > 0 and ω x > 0 are the natural frequencies of the oscillators, ε > 0 is the strength of the coupling and η is white Gaussian noise with standard deviation σ = √ 2E where η(t) = 0, η(t)η(τ ) = δ(t − τ )E.Note that when ε = 0 the system is reduced to φx = ω x + η(t) and becomes non-chronotaxic; when η = 0 and ε > |ω x − ω p | the system becomes chronotaxic with ϕ A x (t) = ϕ p (t) − arcsin((ω p − ω x )/ε).The system was integrated using the Heun scheme [15], with an integration step of 0.001 and noise strength σ = 0.3.∆ϕ x , shown in Figure 2, was obtained by subtracting the unperturbed phase (ϕ A x (t) and ω x t in the chronotaxic and non-chronotaxic cases, respectively) from the perturbed phase ϕ x , as obtained numerically.DFA was then performed on ∆ϕ x , with exponents shown in Figure 2. The values of the exponents demonstrate the differences in the noise distributions between chronotaxic and non-chronotaxic systems.In the chronotaxic case, the noise is closer to white, whereas in the non-chronotaxic case it is closer to a random walk.It is this difference which is exploited in the PFA method.
In many systems, particularly those originating from nature, there will be more than one oscillation present in a signal, with different chronotaxicity characteristics.To test the PFA method in the case of multiple modes a signal containing two distinct oscillations was simulated, with dynamics described by Equation ( 21), with time varying angular frequencies, where f p and f x are the average frequencies of oscillation in Hz of the chronotaxic and non-chronotaxic case, respectively, and f m is the frequency of variation.Frequencies of oscillation were chosen to vary around 1 and 0.25 Hz in the non-chronotaxic and chronotaxic cases, respectively, with f m = 0.003.Both systems were perturbed with white Gaussian noise with strength σ = 0.5.The logarithmic frequency scale of the wavelet transform is very useful for identifying and separating the presence of oscillatory modes, which may otherwise appear as merged in other time-frequency representations, such as the windowed Fourier transform.Figure 3 shows the results of PFA on the signal.It correctly identifies mode A (around 0.25 Hz) as chronotaxic, and mode B (around 1 Hz) as non-chronotaxic.In single variable time series obtained from real dynamical systems, it is highly unlikely that the observed dynamics will result from a simple, unidirectional constant coupling as described above.Rather, the system may be influenced by continuous perturbations, couplings to other oscillators, and temporal fluctuations in chronotaxicity.Here, we demonstrate the applicability of the described inverse approach methods to these more complex cases.We model a system of two bidirectionally coupled oscillators: with drivers ϕ p1 and ϕ p2 , and ω p1 (t) = 2π − 0.5 sin(2π0.005(t))and ω p2 (t) = π − 0.5 cos(2π0.005(t)).First, we consider the case of strong influence of the driver ϕ p1 on the system, resulting in chronotaxicity of both oscillators.Phase fluctuation analysis was applied to the system, and successfully identified both ϕ x1 and ϕ x2 as chronotaxic (see Figure 4a).Second, we consider the case in which ϕ x1 is chronotaxic but ϕ x2 is not, and demonstrate that despite continuous influences from multiple drivers and other oscillators, single variable time series arising from the same system can be distinguished in terms of their chronotaxic dynamics.Again, PFA correctly distinguishes between the two oscillators.This could be of great importance when investigating composite parts of a larger dynamical system, and seeking to identify causal relationships between observed oscillations.For example, recent advances in cellular imaging are providing the means to observe the dynamics of individual cellular processes in different cellular compartments [34].Applying inverse approach methods for the detection of chronotaxicity to these dynamics could provide valuable information on the current state of the cell.x2 was extracted with f 0 = 2 and smoothed using a polynomial fit (red line), whilst ϕ * x2 was extracted from the wavelet transform with f 0 = 0.5 (grey line).Bayesian inference was applied, using a time window of 90 s.The inferred direction of coupling can be seen in (c).Positive values show coupling from the driver to the oscillator only.(d) I sync was calculated and shows excellent agreement with changes in ε 3 .I chrono was also calculated, and was slightly less accurate due to the direction of coupling becoming negative very briefly, because of reduced information flow between systems to accurately infer parameters during synchronization.
So far, we have only considered scenarios in which a system constantly remains as either chronotaxic or non-chronotaxic.Real dynamical systems may exhibit time variation in their coupling strengths, allowing the system to fluctuate between chronotaxic states.In these cases, it is possible to use dynamical Bayesian inference to track variations in chronotaxicity in time.To demonstrate this, ε 3 was allowed to vary in time in Equation ( 23), whilst ε 1 = ε 2 = 0.1 and ε 4 = 0, resulting in intermittent chronotaxicity of the oscillator ϕ x2 .ϕ A * x2 and ϕ * x2 were extracted from the synchrosqueezed wavelet transform of sin(ϕ x2 ).Results of the application of dynamical Bayesian inference are shown in Figure 5.This method is able to track the intermittent changes in chronotaxicity, through changes in synchronization and direction of coupling, demonstrating its usefulness for the detection of chronotaxicity in systems where the interactions between oscillators are time-varying.

Practical Considerations
Both presented methods, phase fluctuation analysis and dynamical Bayesian inference, rely on precise phase extraction of the estimated attractor ϕ A * x and the perturbed dynamics ϕ * x .Therefore, the parameters in the respective methods must be carefully selected depending on the characteristics of the given data.
The continuous wavelet transform provides an optimal compromise between time and frequency resolution.In the majority of examples used in this paper, f 0 = 1 has been used.However, the wavelet central frequency, f 0 , can be altered to suit specific needs.For example, in a case where there are many phase slips, it may be necessary to extract the estimate of the attractor, ϕ A * x , with a higher f 0 to obtain a better frequency resolution and smoother dynamics, whilst the perturbed phase ϕ * x is extracted using a lower f 0 , leading to an increased time resolution for the purpose of locating each phase slip.The parameter f 0 may also be increased to provide greater distinction between oscillatory modes, but this will be at the expense of time resolution.It should be noted that modes must be separable in time-frequency representations in order for these inverse approach methods to be applicable.
One fundamental assumption of chronotaxicity is that the system under consideration is oscillatory.Although the presented methods can be applied to any extracted phases, one should take great care to ensure that these phases correspond to a true oscillatory mode, otherwise all results will be meaningless.In the numerical simulations presented here, we predetermine the characteristics of the oscillations which are present, and ensure that they are not concealed by noise, allowing their successful extraction directly from the wavelet transform.These extracted phases can be verified using the specified parameters as a reference signal, and thus we can be confident with the final results.On the contrary, in real experimental data, the first question must be whether the signal contains any significant oscillations at all.To determine whether this is the case, the recently developed method of nonlinear mode decomposition (NMD) may be used.NMD is an adaptive, time-frequency representation based decomposition tool, which decomposes any given signal into a set of physically meaningful oscillations (if present) and residual noise.In the detection of chronotaxicity, the crucial advantage of NMD over other decomposition methods, such as EMD or bandpass filtering, is its use of surrogate data testing in order to distinguish between deterministic and random activity [35].The success of surrogate testing for the identification of nonlinear oscillatory modes in neural data has also been demonstrated previously [36], and more generally in [37].By verifying the presence of oscillations, and their underlying nature, e.g., whether they are nonlinear, these methods reliably inform the user which analysis approach to take.In this way, we can ensure that any oscillatory modes extracted from real experimental data are physically meaningful, and their characteristics, including their instantaneous phase, are accurately determined.Once a significant oscillatory mode has been located and extracted using NMD, its smoothed instantaneous frequency provides ϕ A * x for use in phase fluctuation analysis.ϕ * x can then be extracted from the wavelet transform as before, with the parameter f 0 chosen to give sufficient time resolution to follow the noise fluctuations which are removed by NMD.An example of the use of NMD in PFA is provided in Figure 6, and explained further in Section 4.3.x − ϕ A * x .The DFA exponent was calculated and was 1.57, suggesting that the system is not chronotaxic.Checking for phase slips in (g) shows no change in distribution.
The reliability of the presented inverse approach methods increases with data availability, i.e., a longer time series will give a more reliable result.However, it is not often feasible to collect hundreds of cycles of oscillation.When recording data from live subjects, for example blood flow recordings, the time of recording must be a compromise between long time series and subject comfort.In the case of cellular recordings, such as cell membrane recordings via the patch clamp technique, the health of the cell can rapidly deteriorate, and thus affect the reliability of results.Therefore, it is useful to determine the lowest possible number of recorded oscillations for which we may still reliably test for chronotaxicity.
In order to address this question, two unidirectionally coupled phase oscillators (Equation ( 21)) were simulated for 1000 cycles with frequencies 1 and 0.1 Hz, with h = 0.01 and σ = 0.07.With coupling ε = 2, the system is chronotaxic.The important parameters to consider in DFA are n min and n max , the lower and upper values for the range of the first order polynomial fits performed in order to calculate ∆ϕ x .The lower value, n min , is set to be 2 cycles of the slowest oscillation, to ensure observation of the dynamics over a longer range than one cycle.The smallest n max required to still obtain a reliable DFA exponent was observed to be n max = 3 cycles of oscillation (see Figure 7), provided that the time series is sufficiently long.The second test seeks to identify the required length of the whole time series when Based on these results, the time series should be at least 8 times n max , thus, there should be at least 24 oscillations in the time series.However, to ensure universal applicability, the length of the time series should be at least 10 times n max , the generally accepted value in DFA [38], resulting in the requirement of 30 cycles.using these values of n min and n max in DFA.The DFA exponent was calculated from varying lengths of the same noise signals, from 3 to 10 times n max , to identify the point where the result is no longer reliable.It was found that the time series should be at least 8 times n max in order to obtain a reliable result, therefore at least 24 cycles of the slowest oscillation are required to test for chronotaxicity.However, if possible, the time series should be at least 10 times n max [38], to reduce noise by providing more data windows.Overlapping within DFA is also possible, and will go some way toward reducing noise, and improve reliability.The results shown in Figure 7 were obtained with an overlap of 0.8.
Whilst we expect the value of the DFA exponent α to be around 0.5 in a chronotaxic system, and 1.5 in a non-chronotaxic system, it is unlikely to be so definitive in reality.In fact, the value of α will depend on a number of factors.The type of noise in a real system is not necessarily white, however the point of phase fluctuation analysis is to identify changes in its distribution.α will also vary depending on how strong the chronotaxicity is in the system, i.e., how strongly driven the observed oscillator is.In our models, this can be represented by varying the coupling strength ε; weaker coupling will result in a higher DFA exponent as the noise is partially integrated.The ratio of the natural frequency of the chronotaxic oscillator to the frequency of the external driver, or detuning, may also affect the value of α.

Application to Experimental Data
Chronotaxicity will manifest in nature as a result of a driving system which is strong enough that the oscillatory response system maintains stability in its frequency and amplitude, even when subject to continuous external perturbations.Chronotaxicity was previously demonstrated in the heart rate variability (HRV) [14], when influenced by paced breathing.It has been shown that the main direction of coupling between the cardiac and respiratory oscillators is the influence of the respiration on the heart rate, known as respiratory sinus arrhythmia (RSA), and this was clearly demonstrated.Here, we provide an example of the application of phase fluctuation analysis to real experimental data in the form of an electroencephalogram (EEG) recording from an anaesthetised human subject.
Distinct oscillations have long been observed in the brain since the invention of EEG by Hans Berger in 1924.Briefly, from lowest to highest frequency, there are at least five frequency bands which have been identified in approximately the following frequency intervals: delta (0.8-4 Hz), theta (4-7.5 Hz), alpha (7.5-14 Hz), beta (14-22 Hz) and gamma .Different frequencies of oscillation have been attributed to distinct states of the brain.For example, the alpha and theta bands have been shown to reflect cognitive and memory performance [39].One active area of research utilising the information provided by these oscillations is in attempts to quantify the depth of anaesthesia based on their temporal evolution in different states of consciousness.Despite the daily worldwide use of general anaesthesia (GA), the mechanisms leading to this state are still poorly understood in terms of how it truly affects the brain.Thus, brain-state monitoring is still not an accepted practice in GA, due to the lack of reliable markers [40].However, recent studies in which the spectral power of the oscillations in different frequency bands has been tracked both temporally and spatially during anaesthesia with propofol have shown promising results.For example it was shown that during consciousness, alpha oscillations are concentrated in occipital channels, whilst during propofol induced anaesthesia, these oscillations are concentrated in frontal channels [40].An increase in power in the frequency interval 0.1-1 Hz (delta) was also observed in this study during anaesthesia.Understanding the mechanisms underlying these changes in brain function could not only lead to new approaches to anaesthesia monitoring but may be widely applicable in many areas of neuroscience, including in the study of various neurological disorders.
It has been clearly demonstrated that phase interactions are highly important for healthy brain functioning, with by far the most widely reported observations revolving around phase synchronization, which can, as an example, be used to infer information about short and long range behaviours [41].Brain waves arise from networks of synchronized neurons, and the detected phase of these oscillations determines the degree of excitability of the neurons, and influences precise discharge times of the cells in the network, therefore affecting relative timing of action potential in different brain regions [42].
Before any conclusions may be drawn about the phase dynamics of a system, the phase must be accurately extracted from the time series.The problem of the extraction of phase from EEG data has been approached from many directions, some more physically meaningful than others.Early approaches to the investigation of phase interactions between brain waves used spectral coherence, but this does not separate phase and amplitude components, thus amplitude effects may influence coherence values when only phase locking information is required [43].A widely used phase extraction approach is the use of the Hilbert transform to obtain the analytic signal [44], usually preceded by band-pass filtering in the frequency interval of interest, highlighting the necessity of the separation of the oscillation of interest from background brain activity-either other oscillations or noise.Lachaux, et al. recognised the necessity of the separation of amplitude and phase when seeking to detect synchrony between brain waves, introducing phase-locking statistics (PLS) [43] to measure the phase covariance between two signals, verified by surrogate testing.This method also allows for non-stationarity in the signal.However, based on very narrow band-pass filtering, this method does not allow for time-variability in the frequency of oscillation, but did highlight the usefulness of complex wavelets in the extraction of phase dynamics.The Hilbert transform and wavelet convolution methods were compared in the analysis of neural synchrony, and found not to differ substantially [45], but both of these methods relied on narrow band-pass filtering beforehand.However, the use of band-pass filtering to extract an oscillatory EEG component with a time-varying frequency has limited usefulness.An instantaneous frequency defined from the analytic signal obtained from band-passing in a particular frequency range in a real signal containing multiple spectral components and noise may be ambiguous and meaningless [41].To address this problem, ridge extraction methods [29] applied to the complex wavelet transform were used to track the instantaneous frequency of a single oscillatory mode [41], providing a much higher precision of phase extraction, and importantly allowing the phase dynamics of nonautonomous systems to be accurately traced in time.Another rarely considered issue when tracing instantaneous frequencies in time is the presence of high harmonics in the signal.Narrow restriction of the frequency range will remove these harmonics, and thus remove valuable intra-cycle phase information.This issue has been addressed directly by the introduction of nonlinear mode decomposition [35].The inverse approach methods applied here take into account all these issues in order to accurately extract the instantaneous phase of brain oscillations.
In order to demonstrate the method and search for evidence of chronotaxicity in the phase dynamics of brain waves we applied phase fluctuation analysis to a real EEG signal.The EEG of an anaesthetised subject was recorded for 20 minutes at 1200 Hz (Figure 6a).The signal was resampled to 100 Hz by splitting the time series into windows, and setting their mid-point to their mean.As expected, strong oscillations were observed in the alpha and delta frequency bands.Nonlinear mode decomposition (see Section 4.2) extracted the oscillatory mode around 10 Hz in the alpha frequency band and identified it as physically meaningful through surrogate testing (Figure 6c).The instantaneous frequency of this mode was then smoothed using a moving average of 4 s.This value was chosen to provide the best match between the instantaneous phase of the extracted nonlinear mode ϕ x and its smoothed version ϕ A * x .As NMD by nature removes the noise from the modes which it extracts, ϕ * x must then be extracted from the continuous wavelet transform with a time resolution which will allow the noise fluctuations to be included in the extracted mode.Here, it is very important to check that the extracted phase corresponds to that extracted using NMD (see Figure 6e).Once the viability of the extracted fluctuations is confirmed, ∆ϕ x can be calculated as ϕ * x − ϕ A * x .The DFA exponent of ∆ϕ x was then calculated, and was 1.57.The distribution |d∆ϕ τ x | was calculated to check for phase slips in the extracted phase fluctuations, but the distribution did not change over any time scale τ .
The analysis suggests that the alpha oscillation as extracted is not chronotaxic.However, the current inverse approach methods are based on a single point attractor and a single response system.As discussed by Sheppard, et al. [46], the spectral peaks observed in the EEG, including those observed in the alpha band, result from frequency synchronization between thousands of neurons.In this sense, the observed phase is, in fact, only a statistical measure highlighting the preferred phase of the underlying ensemble of neurons.A method to quantify this was provided by the mean-field variability index, κ, which changes depending on the interactions in the observed network of oscillators [46].For a non-interacting network, with purely random phasors, κ will converge to 0.215, whereas in a state of complete phase synchronization, κ will tend to zero.Based on the current assumptions of the inverse approach methods, if the detection of chronotaxicity relied only on phase dynamics, we would expect the value of κ to tend to zero in a chronotaxic system.However, when applied to real EEG data, κ was actually greater than 0.215 in most cases, suggesting amplitude synchrony (possibly intermittent), intermittent phase coherence, or both.Therefore, it is apparent that in the case of brain dynamics, to truly characterise chronotaxicity, it must be reconsidered within a network of many oscillators, as known to be present in the brain.Here, the driving system may be a subnetwork of synchronized oscillators or the mean-field or mean-phase of ensembles of neurons, influencing other areas of the brain in complex ways, with both temporal and spatial dynamics to take into account.
The presented methods are restricted by the fact that they are currently only applicable to determining chronotaxicity in phase dynamics.Traditionally, in brain dynamics, it is the amplitude of the oscillations observed in the distinct frequency bands which receives the most attention, although phase dynamics is now gaining considerable recognition [47].In addition to the dynamics within individual frequency bands, there are also interactions between frequency bands [48], known as cross frequency coupling (CFC).The importance of phase information in oscillatory brain activity has been clearly demonstrated, for example phase synchronization between frequencies has been shown to be correlated with certain cognitive processes [49].Phase measures also provide the advantage of high temporal precision [49].However, the nature of CFCs have not only been observed as phase-phase interactions [50], but also as amplitude-phase [51] and amplitude-amplitude interactions [52].Whilst some efforts have been made to isolate phase information in neural oscillations [53], the importance of amplitude-phase interactions cannot be ignored, for example the observed modulation of gamma amplitude by the phase of theta oscillations has been identified as a code utilised in multi-item formation in the brain [54].Other functional roles of amplitude-phase coupling have also been highlighted [55], thus it is clear that both amplitude and phase must be considered simultaneously to accurately characterise brain dynamics.Indeed, phase-amplitude coupling has been demonstrated during anaesthesia [56], meaning that the current inverse approach methods may be insufficient to determine chronotaxicity in this system.

Discussion
The recent formulation of chronotaxic systems provides a completely novel approach to the characterisation of time-varying dynamics in real data.Crucially, they provide a framework in which systems may be time-varying, both in terms of their amplitude and phase dynamics, continuously perturbed, and yet still exhibit determinism.Whilst the apparent complexity of some real time-varying oscillatory systems previously led to their consideration as stochastic or chaotic, chronotaxicity facilitates a much more natural approach to the description of their dynamics.The introduction of this approach required the development of new inverse approach methods for the detection of chronotaxicity in time series arising from dynamical systems.Here, we reviewed the currently available methods for the identification of chronotaxicity from a single time series, and also expanded on various issues regarding their implementation, in order to facilitate the application of the methods to any data set containing at least one oscillatory component.This ability to characterise oscillations in terms of their chronotaxicity, i.e., to determine whether the observed dynamics arise as a result of influence from an external driver, provides the potential to unlock new information about dynamical systems and their interactions with their environment.
As they currently stand, the inverse approach methods for the detection of chronotaxicity are only applicable to systems in which the amplitude and phase dynamics are separable, as they are applied directly to the extracted phases of the system and all amplitude information is discarded.This assumption is valid if considering that the amplitude dynamics of a chronotaxic system corresponds to the convergence of the system to the limit cycle, influenced only by a negative Lyapunov exponent and external perturbations.In contrast, the phase dynamics corresponds to convergence to the time-dependent point attractor, which is also characterized by a negative Lyapunov exponent and external perturbations, but also the motion of the point attractor itself [14].As it is this point attractor in phase dynamics that we are interested in-separation of amplitude and phase follows naturally.Using this approach, an example of chronotaxic dynamics was succesfully demonstrated in a real system, in the case of heart rate variability [14].However, in generalized chronotaxic systems [12], the amplitude and phase are not required to be separable, providing even greater applicability to real systems, allowing amplitude-amplitude and amplitude-phase interactions in addition to the phase-phase dynamics considered in [10,11].Therefore, the incorporation of the ability to identify these new possibilities for chronotaxicity is crucial in the further development of these inverse approach methods.This will then provide the means to detect chronotaxicity in systems where amplitude and phase are not separable, as previously discussed in the case of brain dynamics (see Section 4.3).The current definition of chronotaxicity is based on a time-varying point attractor, exerting influence over a system such that it can remain stable despite continuous external perturbations.Numerical results presented here assume that this point attractor results from a single oscillatory drive system acting on a maximum of two coupled oscillators.However, as highlighted in the brain dynamics example, in reality we must consider that this point attractor could result from multiple interacting influences, for example a network of oscillators, perhaps acting as one synchronized drive system.
Regardless of the mechanisms of the underlying oscillations, if they manifest as a point attractor, characterisation of their chronotaxicity necessitates the application of methods which can extract both their phase and amplitude dynamics with utmost accuracy.Methods reliant on averaging will not provide the required precision.Both amplitude and phase information can be extracted from the continuous wavelet transform, a fact which may be utilised in the further development of inverse approach methods for the detection of chronotaxicity.Extending these methods to simultaneously take into account both phase and amplitude dynamics whilst incorporating the effects of their couplings, may lead to a method based on an optimal combination of time-frequency representations and effective connectivity methods such as dynamical Bayesian inference.This will then provide even wider applicability to real oscillatory systems such as those observed in brain dynamics.

3 .
Identifying chronotaxicity in signals with more than one oscillatory mode.(a) The first 250 s of a time series of a simulated signal containing two distinct oscillations, with coupling strengths ε = 2 for mode A (chronotaxic) and ε = 0 for mode B (non-chronotaxic).(b) The continuous wavelet transform of the signal in (a).(c) The instantaneous frequency (light grey) of both components is extracted from the wavelet transform, with central frequency f 0 = 0.5, and smoothed (red), using a polynomial fit.The smoothed frequency is then integrated in time to obtain an estimate of the unperturbed phase, ϕ A * x , which is then subtracted from the perturbed phase ϕ * x as extracted directly from the wavelet transform.(d) and (e) show ∆ϕ x = ϕ * x − ϕ A * x for each mode.(f) and (g) show the results of DFA on ∆ϕ x , with DFA exponents α correctly identifying mode A as chronotaxic and mode B as not chronotaxic.

Figure 5 .
Figure 5. Identifying intermittent chronotaxicity using dynamical Bayesian inference.Bayesian inference was performed on ϕ * x2 and ϕ A * x2 extracted from sin(ϕ x2 ) (see Equation (23)) with ε 3 varying as shown in (d).(a) the CWT of sin(ϕ x2 ).(b) Instantaneous frequencies extracted from the wavelet transform.ϕ A *x2 was extracted with f 0 = 2 and smoothed using a polynomial fit (red line), whilst ϕ * x2 was extracted from the wavelet transform with f 0 = 0.5 (grey line).Bayesian inference was applied, using a time window of 90 s.The inferred direction of coupling can be seen in (c).Positive values show coupling from the driver to the oscillator only.(d) I sync was calculated and shows excellent agreement with changes in ε 3 .I chrono was also calculated, and was slightly less accurate due to the direction of coupling becoming negative very briefly, because of reduced information flow between systems to accurately infer parameters during synchronization.

20 20 Figure 6 .
Figure 6.An example of the application of phase fluctuation analysis to an electroencephalogram (EEG) signal obtained from the forehead of an anaesthetised patient, shown in (a).(b) The continuous wavelet transform of the electroencephalogram (EEG) signal in (a).(c) Using NMD, a significant oscillatory mode in the alpha frequency band was identified and extracted (dark grey line).(d) The instantaneous frequency extracted using NMD (grey line), and smoothed using a moving average of 4 s (red line).(e) The extracted phases of the mode from NMD (grey), smoothed NMD (red), and from the CWT (black) with f 0 = 1.5.(f) ∆ϕ x was calculated as ϕ *x − ϕ A * x .The DFA exponent was calculated and was 1.57, suggesting that the system is not chronotaxic.Checking for phase slips in (g) shows no change in distribution.

Figure 7 .
Figure 7.In order to test the reliability of the DFA exponent when reducing n max , the maximum number of cycles of oscillation used in its calculation was varied.(a) Chronotaxic oscillation of 1 Hz.(b) Chronotaxic oscillation of 0.1 Hz.(c) Non-chronotaxic oscillation of 1 Hz.(d) Non-chronotaxic oscillation of 0.1 Hz.The same noise signals were then tested with n max = 3 for different lengths of the time series from 10 to 3 times n max .Based on these results, the time series should be at least 8 times n max , thus, there should be at least 24 oscillations in the time series.However, to ensure universal applicability, the length of the time series should be at least 10 times n max , the generally accepted value in DFA[38], resulting in the requirement of 30 cycles.