Correlations in magnitude series to assess nonlinearities: application to multifractal models and heartbeat fluctuations

TL;DR

This paper analytically links magnitude series correlations to nonlinearities in time series, demonstrating their application in distinguishing heart rate variability during rest and exercise, revealing increased non-linearity during rest.

Contribution

It provides analytical expressions for magnitude autocorrelation in linear Gaussian noise and introduces a non-linearity index applicable to short records without scaling.

Findings

Higher non-linearity in heartbeat signals during rest compared to exercise

Method applicable to short time series without requiring scaling

Supports link between non-linearity and autonomic nervous system activity

Abstract

The correlation properties of the magnitudes of a time series (sometimes called volatility) are associated with nonlinear and multifractal properties and have been applied in a great variety of fields. Here, we have obtained analytically the expression of the autocorrelation of the magnitude series of a linear Gaussian noise as a function of its correlation as well as several analytical relations involving them. For both, models and natural signals, the deviation from these equations can be used as an index of non-linearity that can be applied to relatively short records and that does not require the presence of scaling in the time series under study. We apply this approach to show that the heart-beat records during rest show higher non-linearities than the records of the same subject during moderate exercise. This behavior is also achieved on average for the analyzed set of 10 semiprofessional soccer players. This result agrees with the fact that other measures of complexity are dramatically reduced during exercise and can shed light on its relationship with the withdrawal of parasympathetic tone and/or the activation of sympathetic activity during physical activity.

Equations

$$|x_i|=|y_{i+1}-y_i|.$$

$$F_d(\ell )\propto {\ell }^{\alpha }.$$

$$C_x(\ell )=\frac{⟨x_i\cdot x_{i+\ell }⟩-⟨x_i⟩⟨x_{i+\ell }⟩}{{\sigma }_x^2}=⟨x_i\cdot x_{i+\ell }⟩,$$

$$K_x(\ell )\equiv ⟨x_i\cdot x_{i+\ell }⟩-⟨x_i⟩⟨x_{i+\ell }⟩=C_x(\ell ).$$

$${\sigma }_{|x|}^2$$

$$C_{|x|}(\ell )$$

$$K_{|x|}=\frac{2}{\pi }\left[\sqrt{1-K_x^2}+K_x\arcsin K_x-1\right],$$

$$C_{|x|}=\frac{2\left[C_x\arcsin C_x-1+\sqrt{1-C_x^2}\right]}{\pi -2}$$

$$C_{|x|}=\frac{1}{\pi -2}{C}_x^2+\mathcal{O}(C_x^4)$$

$$C_x(\ell )\simeq \frac{(1-\gamma )(2-\gamma )}{2{\ell }^{\gamma }}\propto \frac{\mathrm{sgn}(1-\gamma )}{{\ell }^{\gamma }}$$

$$H=\frac{\beta +1}{2}=\frac{2-\gamma }{2}$$

$$C_{|x|}\propto \frac{1}{{\ell }^{2\gamma }}$$

$$C_{x^2}(\ell )=\frac{⟨x_i^2\cdot x_{i+\ell }^2⟩-⟨x_i^2⟩⟨x_{i+\ell }^2⟩}{{\sigma }_{x^2}^2}$$

$$K_{x^2}(\ell )\equiv ⟨x_i^2\cdot x_{i+\ell }^2⟩-⟨x_i^2⟩⟨x_{i+\ell }^2⟩=2K_x(\ell {)}^2.$$

$$C_{x^2}=C_x^2$$

$${\rm sgn}(x)=\left\{\begin{array}[]{rll}-1&{\rm if}&x<0\\ 0&{\rm if}&x=0\\ 1&{\rm if}&x>0\end{array}\right.$$

$$C_x=\sin \left(\frac{\pi }{2}{C}_s\right)$$

$$C_s=\frac{2}{\pi }{C}_x+\mathcal{O}(C_x^3)$$

$$C_{|x|}=\frac{2}{\pi -2}\left[\frac{\pi }{2}{C}_s\sin \left(\frac{\pi }{2}{C}_s\right)+\cos \left(\frac{\pi }{2}{C}_s\right)-1\right]$$

$$x_{\mathrm{comp}}(i)=x_{\mathrm{mod}}(i)\cdot x_{\mathrm{sign}}(i)$$

$$C_x(\ell )=C_s(\ell )\frac{(\pi -2)C_{|x|}(\ell )+2}{\pi }$$

$$C_x(\ell )\simeq \frac{2H_2(2H_2-1)}{{\pi }^2{\ell }^{2-2H_2}}\left[\frac{H_1^2(2H_1-1{)}^2}{{\ell }^{4-4H_1}}+2\right].$$

$$\Delta t_{RR,i}=t_{RR,i+1}-t_{RR,i}$$

$$x^{\prime }={\Phi }^{-1}\left[F(x)\right]$$

$$\delta C(\ell )=C_{|x|}(\ell )-C_{|x|,\mathrm{linear}}(C_x(\ell ))$$

$$\Delta =\sum _{\ell =1}^{{\ell }_{\max }}\delta C(\ell {)}^2$$

$$\rho (x,y)$$

$$K_{\mathrm{mag}}$$

$${\int }_0^{\infty }$$

Full text

††preprint: APS/123-QED

I Introduction

In the field of time series analysis, the concept of nonlinearity can be interpreted in different ways Yosi_Geo . An intuitive definition is that nonlinear time series are those generated by nonlinear dynamic equations, i.e. the values of the series depend on time, or on other values of the series, according to nonlinear expressions — squares, logarithms, trigonometric functions, etc. But usually we do not have prior information about this dependence, in fact, in most of the cases the goal is nothing but finding such dynamic equations. Nonlinearity is also frequently defined in terms of the autocorrelation function: a time series is nonlinear when there is dependence between the values of the series at different positions even though its autocorrelation vanishes. Although a bit more complicated, the definition of of Schreiber and Schmitz Schreiber2000 is quite suitable for practical purposes. According to this definition, a time series is linear when its Fourier phases are random, i.e. the series of phases of its Fourier transform is a random number uniformly distributed in the interval $[-\pi,\pi]$. Thus, the presence of nonlinear correlations in a time series can be assessed by means of surrogate data tests: (i) Given a time series, compute its Fourier transform, randomize its Fourier phases and transform it back. The resulting surrogated series preserves the distribution of the data and the linear correlations because its power-spectrum remains unchanged Schreiber2000 . (ii) Some relevant statistics is evaluated in the original as well as in the surrogated signal and, if there is a statistically significant difference between both signals, it means that the original Fourier phases were not random and thus, the original signal was nonlinear, i.e. the null hypothesis of linearity can be rejected. Sometimes instead of accepting or rejecting the null hypothesis, the goal is simply to compare the degree of nonlinearity of two different time series (e.g. records obtained under different physiological conditions) and the value of the statistics is directly used as a measure of nonlinearity.

The autocorrelations in the magnitude series is also a good indicative of the presence of nonlinear correlations. For a given time series $\{y_{i}\}$, $i=1,...,N$, its magnitude series (sometimes also called volatility) is usually defined as the absolute value of the series increments:

[TABLE]

It is defined as the magnitude of the increments rather than the magnitude of the series itself because in most cases the series of increments is fairly stationary while the original series is not. Apart from its utility in revealing nonlinear properties, the magnitude series together with the sign series (magnitude-sign analysis Yosi_PRL ) provides complementary information about the original series: the magnitude measures how big the changes are and the sign indicates their direction.

Once obtained the magnitude series, the standard procedure to quantity its correlations is the use of the Detrended Fluctuation Analysis (DFA). In brief, the DFA method obtains the root mean square fluctuations of the series around the local trend $F_{d}(\ell)$ in all windows of a given size $\ell$ and repeats the procedure for different window sizes. Scaling is present when

[TABLE]

Typically, $\alpha$ is estimated as the slope of a linear fitting of $\log(F_{d})$ vs. $\log(\ell)$. The exponent $\alpha$ quantifies the strength of the correlations present in the time series and is also related to the power spectrum exponent $\beta$ and the autocorrelation function exponent $\gamma$ allegrini ; rangarajan . The scaling analysis of the magnitude series, was first introduced to study nonlinearities in heart-beat fluctuations Yosi_PRL but since then, examples of quantifying nonlinearity using the DFA exponent of the magnitude series can be found in many other fields such as Fluid Dynamics fluid , Geophysical geo01 ; geo02 ; Yosi_Geo and Economical time series economy .

The scaling exponent of the magnitude fluctuations is easy to compute and is also related to the width of the multifractal spectrum Yosi_volatility ; manolo , another quantity also frequently used to unveil the nonlinear properties of a signal Plamen1999 .

Nevertheless, this approach shows three main drawbacks:

• (i)

In order to properly define the scaling exponent $\alpha$, $F_{d}(\ell)$ vs. $\ell$ must show a good fit to a power-law, which is not the case in many natural series. Also, the interpretation of crossovers in $F_{d}(\ell)$ vs. $\ell$ as a signature of the existence of regions with different scaling has been recently challenged. In particular, it has been shown that the evaluation of short-range scaling exponent ($\alpha_{1}$), a quantity widely used in heart rate analysis Pena09 , could be affected by spurious results Holl15 and that $\alpha_{1}$ is strongly biased by the breathing frequency Pandelis09 . Without judging the validity of these criticisms, the truth is that some results obtained with $\alpha_{1}$ are contradictory Sandercock06 . These problems affect DFA in general as a technique to evaluate scaling exponents.

• (ii)

Furthermore, particularizing to the evaluation of the scaling exponent of the magnitude series, we have shown very recently that in some situations DFA does not properly detect the correlations and assigns uncorrelated behavior to correlated magnitude series CarpenaDFA .

• (iii)

It is assumed implicitly that the presence of correlations in the magnitude series is a signature of nonlinearity but, as we show later, even for a linear time series, $C_{|x|}>0$ when $C_{x}\neq 0$.

For these reasons we propose here a different approach: We consider a linear Gaussian noise $\{x_{i}\}$, i.e. a time series whose values follow a Gaussian distribution and its Fourier phases are random, with correlations given by $C_{x}$ and obtain analytically the correlations $C_{|x|}$ of the magnitude series $\{|x_{i}|\}$ which result to depend only on $C_{x}$. Note that $C_{|x|}$ are the magnitude correlations expected in purely linear noises. When analyzing an experimental time series $\{x_{{\rm exp}}(i)\}$, the deviation of the correlation of its magnitude $C_{|x_{\rm exp}|}$ with respect to the linear expectation $C_{|x|}$ is then a good signature of nonlinear correlations. Taking into account that natural data does not always follow Gaussian distributions, prior to the computation of the magnitude correlations, we transform the distribution of the data to a normal distribution with zero mean and unit standard deviation, ${\cal N}(0,1)$.

This article is organized as follows: Motivated by the fact that in most of the examples cited above, correlated non-stationary natural series are modeled as fractional Brownian motions (fBm), and thus their stationary increments as fractional Gaussian noises (fGn), in section II we obtain the analytical expression of the autocorrelation of magnitude series of a linear Gaussian noise as a function of its autocorrelation as well as a quadratic approximation. We also obtain the corresponding expression for the series of squares, $\{x_{i}^{2}\}$, which is sometimes used to study nonlinear correlations (Sec. II.1) and discuss about the autocorrelation of the sign series and its relation to the autocorrelation of the magnitude series (Sec. II.2). In section III we explore the nonlinear properties of artificial series generated with a model that produces Gaussian noises with multifractal properties and in section IV, as an example of their utility, we apply the relations derived here to the study of heart beat time series during rest and moderate exercise. Section V concludes the paper.

II Autocorrelation of magnitude series

Given a time series $\{y_{i}\}$, with its corresponding series of increments $\{x_{i}\}$, our aim is to obtain the autocorrelation $C_{|x|}(\ell)$ of its magnitude series $\{|x_{i}|\}$ (1) as a function of the autocorrelation of the series of increments $C_{x}(\ell)$ provided that $\{x_{i}\}$ is a linear Gaussian noise, i.e. all $x_{i}\sim{\cal N}(0,1)$ and that only linear correlations are present in the series. Thus, the autocorrelation function at distance $\ell$ is given by:

[TABLE]

where $\langle\cdot\rangle$ denotes average over the series and $\sigma_{x}^{2}$ is the variance of the series.

Under the assumption of $x_{i}\sim{\cal N}(0,1)$ the autocorrelation coincides with the autocovariance

[TABLE]

On the other hand, for the magnitude series we have:

[TABLE]

and we can write for its autocorrelation:

[TABLE]

Taking into account that $x_{i}$ and $x_{i+\ell}$ are two linearly correlated Gaussian random variables, the autocovariance of the magnitude series, $K_{|x|}(\ell)$, can be expressed as a function of $K_{x}(\ell)$ according to Eq. (36) in appendix A:

[TABLE]

and replacing in (6):

[TABLE]

It is easy to check that (8) is an even and positive function which implies that the magnitude of a linear Gaussian noise cannot be anticorrelated (Fig. 1).

If we consider small values of $C_{x}$, Eq. (8) can be approximated by a Taylor expansion and obtain:

[TABLE]

Thus, for small values, $C_{|x|}$ behaves essentially as the square of $C_{x}$. In fact, the error of (9) is around $2\%$ for $C_{x}=0.5$ which makes this approximation virtually correct for most real data. In figure 1 we plot Eq. (8), its quadratic approximation Eq. (9) as well as several examples of artificial series created with Gaussian linear models.

This result is especially interesting when studying the scaling behavior of series with power-law correlations that have been found in a great variety of complex systems. We can characterize these series by their power spectral exponent $\beta$ because most methods of generating power-law correlated Gaussian noises consist in the generation of series with $1/f^{\beta}$ decay in their power spectrum with $-1<\beta<1$ (e.g. multifractal ; Makse ). In particular, these methods are widely used to generate approximate fractional Gaussian noises (fGn) which are indeed linear Gaussian noises whose autocorrelation function decays asymptotically as a power law Beran :

[TABLE]

where $\gamma=1-\beta$. It is also quite common to characterize the fGns by their Hurst exponent ($H$) which is related to both $\beta$ and $\gamma$ by:

[TABLE]

For stationary time series ($0<H<1$), the Hurst exponent also coincides with the DFA exponent $\alpha$ (2). Note that the term ${\rm sgn}(1-\gamma)$ in the numerator of (10) vanishes for $\beta=0$ ($H=0.5$, white noise) thus leading to an uncorrelated random noise. For $\beta>0$ ($H>0.5$) the series is long-range correlated, also known as having “long memory” rangarajan ; Beran , in the sense that its autocorrelation decays very slow with exponent $\gamma<1$. In fact $\sum_{\ell=1}^{L}C_{x}(\ell)$ diverges as $L\rightarrow\infty$. Likewise for $\beta<0$ ($H<0.5$), i.e. $\gamma>1$, $C_{x}$ is negative and the series is anticorrelated. In this situation, although the autocorrelation also decays as a power law, we cannot properly speak about long-range anti-correlations because they decay relatively fast, in the sense that now the autocorrelation function is summable. Another conclusion drawn from (10) is that we cannot obtain linear Gaussian noises with positive autocorrelation functions decaying faster than $1/\ell$.

We obtain from (9) that the autocorrelation of the magnitude series of a fGn also decays as a power law with exponent $2\gamma$ and is always positive, even for $H<0.5$ when the fGn is anticorrelated:

[TABLE]

Nevertheless, we must distinguish two different situations:

• (i)

$H>0.75$. Here $2\gamma<1$ and $C_{|x|}$ decays slower than $1/\ell$ thus leading to long-range power-law correlations in the magnitude series.

• (ii)

$H<0.75$. Now $2\gamma>1$ and $C_{|x|}$, although still being positive and following a power-law, decays very fast. For example, in Fig. 2 we can see that for $H=0.7$, $C_{|x|}$ reaches the background noise level for relatively short scales ($\ell<100$) even for a time series as long as $2^{24}\simeq 1.6\times 10^{7}$.

Indeed, the methods quantifying correlations by means of the study of fluctuations fail to detect the power-law correlations present in magnitude series for $H<0.75$. For example, two widely used techniques like Fluctuation Analysis (FA) or Detrended Fluctuation Analysis (DFA) Peng94 wrongly classify as “white noise” the magnitude series of Gaussian noises with $H<0.75$ manolo ; Yosi_volatility despite being true only for $H=0.5$ CarpenaDFA .

II.1 Relation with the autocorrelation of square series

For simplicity, sometimes the autocorrelation of square series, $\{x_{i}^{2}\}$, is studied instead of the magnitude series, i.e.:

[TABLE]

Indeed, it has been shown numerically that the scaling properties of the correlations of both series are quite similar Yosi_volatility . Below we justify analytically this similarity.

As we did for the magnitude series, first we obtain the autocovariance of the square series, $K_{x^{2}}$, as a function of the autocovariance of the series, $K_{x}$ (Appendix B):

[TABLE]

Taking into account that $x_{i}\sim{\cal N}(0,1)$ and thus $\sigma_{x^{2}}=\sqrt{2}$, we obtain:

[TABLE]

Which obviously implies that the squares of a linear Gaussian noise, just as the magnitude, cannot be anticorrelated.

Eq. (15) also justifies the fact that for power-law correlated series, $C_{|x|}$ and $C_{x^{2}}$ scale asymptotically with the same exponent: for long enough values of $\ell$ we have $C_{x}\ll 1$ and thus the approximation (9) is valid, leading to $C_{|x|}\propto C_{x^{2}}$.

II.2 Relation with the autocorrelation of the sign series

Apart from its relevance in the study of nonlinear correlations, the magnitude series together with the sign series (defined below) provide complementary information about the original series $\{y_{i}\}$: while the magnitude measures how big the changes are, the sign indicates their direction. Sign series are also relevant for the study of first-passage time in correlated processes Conchita12 . Below we obtain a relation between $C_{|x|}$ and the autocorrelation of the sign series, $C_{s}$.

Given a time series $\{x_{i}\}$, the series $\{{\rm sgn}(x_{i})\}$ is defined by:

[TABLE]

If the series of increments $\{x_{i}\}$ is a linear Gaussian noise, Apostolov et al. signum_corr have shown that the autocorrelation of the sign series $C_{s}(\ell)$ can be expressed in terms of the autocorrelation $C_{x}(\ell)$ by:

[TABLE]

Again, we can also obtain an approximation for small values of $C_{x}(\ell)$:

[TABLE]

which implies that, if $C_{x}$ is a power law, $C_{s}$ scales asymptotically with the same exponent as $C_{x}$. In particular, this result holds for fGns (Fig. 3).

In addition and, taking into account that $-1\leq C_{s}\leq 1$, from here it is clear that $C_{x}$ and $C_{s}$ always have the same sign and thus, the sign series will be correlated where $\{x_{i}\}$ is correlated and anticorrelated where $\{x_{i}\}$ is anticorrelated (Fig. 3).

Equation (17), together with (8) allows us to express the autocorrelation of the magnitude series as a function of the autocorrelation of the corresponding sign series:

[TABLE]

III Example of a nonlinear model

Up to now we have only shown examples of linear Gaussian signals for which the derived relations among $C_{x}$, $C_{|x|}$ and $C_{s}$ (Eqs. (8), (17) and (19)) must hold. Nevertheless, if we consider nonlinear Gaussian signals, i.e. signals that, despite having a Gaussian distribution have nonrandom Fourier phases, these relations are no longer valid and the deviation from these equations can be used as a signature of nonlinearity. For example, if $C_{|x|}(\ell)\neq 0$ and $C_{x}(\ell)=0$, i.e. eq. (8) does not hold, two values of the signal at distance $\ell$ are not linearly correlated ($C_{x}(\ell)=0$) but they are not independent because $C_{|x|}(\ell)\neq 0$ and thus, the signal is nonlinear according to one of the definitions given in the introduction.

Here it is important to stress that these equations are valid for each individual value of the autocorrelation function and the possible deviations from nonlinearity can be observed without the assumption of any kind of scaling or power-law behavior in the signal.

We concentrate here on equation (8) because the correlations in the series of magnitudes have been related to the presence of nonlinear correlations and multifractal structure Yosi_PRL ; Yosi_volatility ; manolo . To show the effect of nonlinearities we generate artificial series using a simple method proposed by Kalisky et al. Yosi_volatility which is able to generate multifractal Gaussian noises just by multiplying the sign and the magnitude of two independent linear Gaussian noises. Despite its simplicity, this method is able to independently control both the linear correlations of the signal and its multifractal spectrum width — see also manolo for a systematic exploration of the method.

In brief this procedure, composition method from now on, works as follows:

• (i)

Obtain the magnitude series of a ${\cal N}(0,1)$ fGn $\{x_{{\rm mag}}(i)\}$, with Hurst exponent $H_{1}$ and the sign series of another ${\cal N}(0,1)$ fGn $\{x_{{\rm sign}}(i)\}$, with Hurst exponent $H_{2}$, where $i=1,...,N$, being $N$ the size of the series.

• (ii)

Obtain the composed series as the product of the magnitude and sign series:

[TABLE]

for $i=1,...,N$.

The resulting series $\{x_{\rm comp}(i)\}$ is Gaussian by construction but it presents nonlinear correlations and Eq. (8) is not fulfilled. Instead, it can be shown that its autocorrelation function is given by manolo :

[TABLE]

where obviously $C_{|x|}$ and $C_{s}$ coincide with the autocorrelation functions of the magnitude of the fGn with $H_{1}$ and the sign of the fGn with $H_{2}$ respectively. Note that, although $C_{x}$ is not exactly a power law, it decays asymptotically as $1/\ell^{2-2H_{2}}$, i.e. the autocorrelation of the composed series decays asymptotically with the same exponent as the autocorrelation of the fGn used to obtain the sign series. Indeed, just take into account approximations (9) for $C_{|x|}$ and (18) for $C_{s}$ and the asymptotic expression for the autocorrelation of a fGn (10) to obtain:

[TABLE]

For $0<H_{1},H_{2}<1$ the second summand is the leading one and we get asymptotically $C_{x}(\ell)\propto 1/\ell^{2-2H_{2}}$. In that sense we say that the linear correlations of the composed series are controlled by the sign manolo .

In Fig. 4 we show $C_{|x|}$ vs. $C_{x}$ for several examples of nonlinear series generated by means of the composition method. For all the series shown $H_{2}=0.85$ and thus, all of them have the same scaling behavior for the linear correlations, nevertheless the different values of $H_{1}$ lead to different degrees of nonlinearity according to the deviation of $C_{|x|}$ from the linear expectation (dashed line in Fig. 4). Note that, no matter the value of $H_{1}$, in all cases we observe a deviation from linearity. For smaller values of $H_{1}$ this deviation is more evident at small $C_{x}$ (longer scales) while for large $H_{1}$ it appears mainly at great $C_{x}$ (small scales).

This means that the uncoupling of magnitude and sign (i.e. the magnitude of the changes is independent of its direction) always leads to a nonlinear behavior or, conversely, in a linear Gaussian signal magnitude and sign are not independent but coupled in a specific way that leads to the behavior described by Eq. (19). In general, for natural signals where magnitude and sign are neither independent nor Gaussianly coupled, plots of $C_{|x|}$ vs. $C_{x}$ can be of great utility to shed light about the way in which the magnitude of the changes is related to its direction, i.e. the magnitude and sign coupling.

III.1 Nonlinearity and multifractality

Multifractality and nonlinearity are two concepts that usually go together. Indeed, the width of the multifractal spectrum is considered to be linked to the degree of nonlinearity of the signal Parisi85 ; Badin16 and the finding of multifractal properties is usually associated with complex nonlinear interactions in the systems under study. Nevertheless, although related concepts, multifractality and nonlinearity describe the properties of the signal from different points of view Tang2015 .

The nonlinear signals generated by means of the composition method described above are a good example to show that nonlinearity is not always related to multifractality. This method was originally developed Yosi_volatility to generate Gaussian signals with multifractal properties, in fact, it has been shown that the width of the multifractal spectrum grows linearly with the Hurst exponent $H_{1}$ of the signal used to obtain the magnitude series for $H_{1}>0.75$ and when $H_{1}<0.75$ the width of the multifractal spectrum almost vanishes manolo . Nevertheless, we show here (Fig. 4) that for all values of $H_{1}$ (including the white noise for which $H_{1}=0.5$) the composed signal is clearly nonlinear, despite having an almost zero mutifractal spectrum width. Here, it is important to point out that the region where multifractal detrending techniques give null multifractal width ($H_{1}<0.75)$ manolo coincides with the region where $C_{|x|}$ lies below the linear expectation, at least in the region where the power law fits are carried out ($\ell>4$). According to this, we can say that for this model multifractality is a signature of nonlinearity only when the autocorrelations in the magnitude are larger than expected in a linear model. In the opposite situation, the time series is indeed nonlinear but the multifractal analysis will not reveal it.

IV Example of natural signals: Heart rate during rest and exercise

Since the pioneering works Peng93 , much attention has been paid to the study of correlations in time series of interbeat intervals, i.e. series of times between consecutive heart beats $\{t_{RR,i}\}$, also known as $RR$ time series. In fact, the correlations in such series have been revealed as a powerful tool to evaluate alterations due to disease or aging Goldberger02 ; Plamen1996 , discriminate between physiological states PlamenSleep1999 and assess the state of fitness Aubert03 ; Dong16 . In most cases, studies are limited to linear correlations (power-spectrum, autocorrelation function, DFA, etc.) but nonlinear correlations are indeed present in $RR$ time series Baillie09 and are supposed to play an important role in heart dynamics as their reduction or absence has been related to aging and certain pathological conditions Plamen1999 ; Yosi_PRL . It is worth mentioning that frequently the correlations of the data are supposed to scale as power laws.

Regarding the heart rate during exercise, it is well known that heartbeat dynamics can change dramatically with physical activity. The most evident changes are the abrupt increase in the heart-rate (i.e. reduction of the mean $RR$ intervals) and the reduction of the heart rate variability (HRV), i.e. the variance of the $RR$ times series Sarmiento13 . In addition to these features that can be observed by direct inspection of raw $RR$ time series (Fig. 5.a), it has also been found that exercise modifies the distribution of the power spectrum by reducing the low frequency components Sarmiento13 ; Sandercock06 ; Anosov00 and introducing very high frequencies related to the respiration rate Lewis10 , decreases the sample entropy Platisa08 . Also, the linear correlations measured by the short scale DFA exponent ($\alpha_{1}$) are not only reduced with exercise Karasik02 ; Platisa08 but also can be correlated with the intensity of the exercise Hautala03 . Nevertheless it is fair to say that the opposite result can also be also found in the literature Tulppo2001 . In summary, despite this last contradiction, the general agreement is that in a wide sense the complexity of the $RR$ time series is reduced during exercise and that this effect is related to the breakdown of the equilibrium between the two branches of the autonomic nervous system due to the withdrawal of parasympathetic tone and/or the activation of sympathetic activity (see Sandercock06 ; Lewis10 for reviews).

Here we hypothesize that this reduction in complexity should also be reflected in the lost of nonlinearity in the heart dynamics during exercise. In particular, we focus ourselves on short scales because it has been reported that in this range ($\ell<11$ beats) linear correlations seem to be clearly affected by the intensity of the exercise and because, in practice, the typical length of the records at rest is rarely longer than 10-15 minutes (500-1000 beats) to avoid excessive interferences with the training sessions, thus preventing from accurate evaluation of autocorrelation functions at long distances.

We analyze records during rest and moderate exercise from 10 semi-professional soccer players all of them healthy males (age $23.8\pm 2.9$ yr) without any prior history of cardiovascular disease. Each record includes two stages: (i) 10 minutes of normal wake rest condition, laying in supine position on the soccer field (ii) followed by 20 minutes of moderate running, i.e. at typical warming-up pace (Fig. 5.a). Heart rate was monitored beat-by-beat using a Polar S810i $RR$ cardiotachometer (Polar Electro, Oy, Finland) validez_polar2 .

As $RR$ time series are typically non-stationary, especially during exercise (Figs. 5.b and c), it is a common practice to analyze the series of its increments:

[TABLE]

which are quite stationary, at least in weak-sense (Figs. 5.d and e). Following the notation introduced in Sec. II $\{y_{i}\}$ would be the series of $RR$ intervals while $\{x_{i}\}$ would be the series of interbeat intervals increments ($\Delta t_{RR}$).

The distributions of $\Delta t_{RR}$ are fairly symmetric, although they are not exactly Gaussian but Levy-stable distributions with tails decaying slower than in the Gaussian case Peng93 (Figs. 5.f and g). For this reason, prior to the analysis we convert the distribution of the data to a standard normal distribution by means of the transformation:

[TABLE]

where $F(\cdot)$ is the cumulative distribution of the original $\Delta t_{RR}$ data and $\Phi(\cdot)$ is the cumulative standard normal distribution ${\cal N}(0,1)$. We have observed that this transformation practically does not modify the linear correlations $C_{x}$ (not shown).

For each subject we compute the autocorrelation function of the series of increments $C_{x}(\ell)$ and of the magnitude series $C_{|x|}(\ell)$ for both rest and exercise records.

In Fig. 6 we show the results for one of the subjects for $\ell=1,...,20$. In general, we observe that $C_{x}$ reaches similar values during rest and exercise or even greater values for the latter (Fig. 6.b) but, on the other hand, $C_{|x|}$ is typically greater during rest (Fig. 6.c). In addition, if we inspect carefully Fig. 6.a it is clear that not only the values of $C_{|x|}$ are greater on average for rest than for exercise but also the exercise records are closer to the thick line (expectation for a Gaussian linear noise). For this reason, a good measure of nonlinearity is not simply the autocorrelation in the magnitude $C_{|x|}$ but its difference with the expectation for a linear Gaussian noise computed using Eq. (8) (thick line in Fig. 6.a):

[TABLE]

This quantity takes into account not only the value of $C_{|x|}$ but also its difference with the linear expectation. For example $C_{|x|}(\ell=1)$ reaches a relatively high value for both, rest an exercise (Fig. 6.c) but, once subtracted the linear expectation $\delta C(\ell=1)$ is much higher for rest than for exercise (Fig. 6.d).

In order to obtain a single number to quantify the nonlinearity of a signal we propose here the sum of the squares of the curve $\delta C(\ell)$:

[TABLE]

In particular, as we are interested in the short scale correlations, and following most of the authors in the bibliography, we adopt $\ell_{\rm max}=10$. We obtain that our nonlinearity index $\Delta$ is clearly higher during rest than during exercise (Fig. 7). For each individual subject $\Delta$ is higher for his record during rest than for his corresponding record during exercise and also the group averages are clearly different for rest and exercise ($p<10^{-7}$). Nevertheless, we have to take into account that when dealing with relatively short records, comparisons between series of different length can lead to spurious results due to finite size effects. Here, we have that the records during exercise are two times longer than those during exercise; in addition due to the fact that HR increases with physical activity, the records during exercise are 4-5 times longer in number of beats. For this reason, we check the validity of our findings by comparing our records during rest with records of the same number of beats during exercise: consider a subject with a $N_{r}$ beats record during rest and a corresponding record during exercise of length $N_{e}>N_{r}$ beats and let be $n=\lfloor N_{e}/N_{r}\rfloor$. We extract $n$ non-overlapping windows from the exercise record starting from left to right and another $n$ non-overlapping windows from right to left (in order to use all available data). For all these $2n$ sub-series we compute $\Delta$ and average for each subject. Results are shown in Fig. 7 (red triangles). Although now the differences between rest and exercise are a bit smaller, all the values of $\Delta$ for rest are above the corresponding values for exercise (including the error bars) and the difference between group averages is still statistically significant ($p=3\times 10^{-4}$ ).

V Conclusions

We have obtained analytically the expression of the autocorrelation of the magnitude series $C_{|x|}$ of a linear Gaussian noise as a function of its autocorrelation $C_{x}$ as well as several analytical relations involving $C_{|x|}$, $C_{x}$ and the autocorrelation of the sign series $C_{s}$. These expressions are useful to study the nonlinear properties of artificial series obtained by models as well as natural series with the great advantage that our approach does not make any prior assumption about the scaling or functional form of the autocorrelation functions. Indeed, the nonlinearity index proposed in section IV has the advantage that can be evaluated on relatively small samples and does not require scaling in the autocorrelation function.

In particular, we study the nonlinear properties of a Gaussian model designed to produce series with multifractal properties and show that this model generates nonlinear signals for all the values of the parameters even for those leading to monofractal behavior. This means that, although multifractality seems to imply nonlinearity, the reverse is not always true.

We also analyze natural time series. Specifically, we have shown that the heart-beat records during rest show higher nonlinearities than the records of the same subject during moderate exercise. This behavior is also achieved on average for the analyzed set of 10 semiprofessional soccer players. With this result we show that the nonlinear properties of the heart-beat dynamics is yet another feature supporting that the complexity of the heart-beat is reduced during exercise. It is also worth mentioning that our nonlinearity index is sensible to moderate exercise. This means that it could probably be applied to the study of nonlinear properties during exercise at different levels of intensity and thus, it could be of interest to study the changes in the balance between sympathetic and parasympathetic nervous systems during exercise.

VI Acknowledgments

This work is partially supported by grants: FQM-7964 and FQM-362 from the Spanish Junta de Andalucía. P.B. thanks Oleg Usatenko for helpful discussions.

Appendix A Autocovariance of the magnitudes of two linearly correlated Gaussian variables

Consider two random variables $\{X,Y\}$, both with zero mean and unit standard deviation and following the bivariate Gaussian distribution Multivariate :

[TABLE]

where $K=\langle xy\rangle$ is the covariance of variables $X$ and $Y$, which also coincides with their correlation taking into account that both of them have zero mean and unit standard deviation. Note that from (A) it follows that $K=0$ if and only if $\{X,Y\}$ are independent, i.e. they have only linear correlations.

The covariance of $|X|$ and $|Y|$ is given by:

[TABLE]

where we have used that $\langle|x|\rangle=\langle|y|\rangle=\sqrt{2/\pi}$.

Now, changing the integration variables $\xi=x/A_{k}$ and $\varphi=y/A_{k}$, where $A_{k}\equiv\sqrt{1-K^{2}}$, we obtain:

[TABLE]

The integral over $\varphi$ can be written as

[TABLE]

where we have used the identity Gradshteyn :

[TABLE]

with $a=K\xi$, $b=1/2$ and the fact that

[TABLE]

Replacing (A) in (A)

[TABLE]

and using the identity Ng

[TABLE]

with $a=\frac{K}{\sqrt{2}}$ and $b=\frac{A_{k}}{\sqrt{2}}$, we get:

[TABLE]

Finally, after some trigonometric manipulation:

[TABLE]

Appendix B Autocovariance of the squares of two linearly correlated Gaussian variables

Considering again, as in Appendix A, two Gaussian variables $\{X,Y\}$ following the bivariate Gaussian distribution (A), the autocovariance of their squares is given by:

[TABLE]

where we have used that $\langle x^{2}\rangle=\langle y^{2}\rangle=1$. Now, changing the integration variables $\xi=x/A_{k}$ and $\varphi=y/A_{k}$ we obtain:

[TABLE]

Taking into account that $A_{k}^{2}=1-K^{2}$:

[TABLE]

and finally:

[TABLE]