Dynamical Heart Beat Correlations during Running
Matti Molkkari, Giorgio Angelotti, Thorsten Emig, Esa Räsänen
DDynamical Heart Beat Correlations during Running
Matti Molkkari, ∗ Giorgio Angelotti, † Thorsten Emig,
2, 3, ‡ and Esa Räsänen
1, 2, § Computational Physics Laboratory, Tampere University, FI-33720 Tampere, Finland Laboratoire de Physique Théorique et Modèles Statistiques, CNRS UMR 8626,Université Paris-Sud, Université Paris-Saclay, 91405 Orsay cedex, France Massachusetts Institute of Technology, MultiScale Materials Science for Energy andEnvironment, Joint MIT-CNRS Laboratory (UMI 3466), Cambridge, Massachusetts 02139, USA (Dated: 2020-07-03)Fluctuations of the human heart beat constitute a complex system that has been studied mostlyunder resting conditions using conventional time series analysis methods. During physical exercise,the variability of the fluctuations is reduced, and the time series of beat-to-beat RR intervals (RRIs)become highly non-stationary. Here we develop a dynamical approach to analyze the time evolutionof RRI correlations in running across various training and racing events under real-world conditions.In particular, we introduce dynamical detrended fluctuation analysis and dynamical partial autocor-relation functions, which are able to detect real-time changes in the scaling and correlations of theRRIs as functions of the scale and the lag. We relate these changes to the exercise intensity quantifiedby the heart rate (HR). Beyond subject-specific HR thresholds the RRIs show multiscale anticor-relations with both universal and individual scale-dependent structure that is potentially affectedby the stride frequency. These preliminary results are encouraging for future applications of thedynamical statistical analysis in exercise physiology and cardiology, and the presented methodologyis also applicable across various disciplines.
INTRODUCTION
The increasing popularity and accuracy of wearabledevices and sensors present new opportunities to studyhuman physiology in a continuous, non-invasive man-ner for a huge number of subjects under real-world con-ditions. These devices enable the measurement of aplethora of physiological and mechanical signals such asthe heart rate, beat-to-beat (RR) intervals, overall mo-tion via GPS, motion of specific body locations via ac-celerations, and skin temperature. These data can berecorded in real time, often at one second intervals, anduploaded to web services. To date, most recorded dataare not analyzed in scientific rigour due to a lack of suit-able models for the dynamics of physiological signals un-der various intensities of exercise load, and also due torestricted availability of the data (property of industryand users). This limits opportunities for a better under-standing of complex physiological processes, diagnosticsand monitoring for patients in rehabilitation, and the op-timal training of athletes.However, it has been long known that a variety ofphysiological conditions and cardiac diseases affect heartrate variability (HRV) and the correlations in RR inter-vals [1]. In exercise physiology, HRV is often used atrest to evaluate recovery, fatigue and overtraining. Itis known that during exercise the overall variability ofthe RR intervals (RRI) is strongly suppressed. Regard-less, the RRI correlations contain valuable information ∗ matti.molkkari@tuni.fi † [email protected] ‡ [email protected] § esa.rasanen@tuni.fi even during exercise [2–4]. For example, the possibil-ity to determine certain physiological thresholds, such asthe anaerobic threshold, from the frequency spectrum ofHRV has been examined [5, 6]. Often the relative im-portance of low-frequency (LF: 0.04–0.15 Hz) and high-frequency (HF: 0.15–0.4 Hz) spectral power is studiedduring exercise. Using this concept as a measure of therelative sympathetic (SNS) and parasympathetic nervoussystem (PNS) activity, it has been shown that the PNSactivity decreases dramatically during exercise [7]. Incontrast, the SNS activity remains unchanged past thefirst ventilatory threshold before increasing abruptly [7].However, the use of the HF/LF ratio to measure cardiacsympatho-vagal balance has been criticized [8]. More-over, it is known that Fourier decomposition of dynamicsignals is often hampered by non-stationarity.To overcome the complications of Fourier methods andnon-stationarities, we base our analysis on detrendedfluctuation analysis[9] (DFA), which was developed tomeasure correlations in non-stationary time series by uti-lizing systematic detrending [9–11]. Furthermore, we areinterested in analyzing real world exercises recorded withreadily available commercial sports watches. Hence, westudy real-time correlations of RRIs during marathonraces (group M) and freeform training runs (group T).Such uncontrolled data may be plagued by severe non-stationary conditions, and the conventional division intoshort- and long-scale DFA exponents[10, 12] is likelyto be insufficient. To this end, we introduce dynamicDFA (DDFA) for the accurate determination of time-and scale-dependent scaling exponents α ( t, s ) with hightemporal resolution. To check the consistency of ourmethodology, we also apply similar dynamic approachto partial autocorrelation functions (PACFs) to obtaintheir dynamic counterpart (DPACF). a r X i v : . [ phy s i c s . d a t a - a n ] J u l RESULTS
Our main result is the discovery of scale-dependentanticorrelations ( α < . ) in the RRIs during runningthat vary with the heart rate. The anticorrelations ap-pear after the HR exceeds a subject-specific threshold.Their magnitude and the scale with the most dominantanticorrelations changes with exercise intensity. We findthat the DDFA method can reliably determine the dy-namic, scale-dependent scaling exponent α ( t, s ) (pleasesee the Supplementary Information for its numerical val-idation). Hence, it provides a powerful method for mea-suring multiscale correlations of non-stationary physio-logical signals. The results from the DDFA and DPACFmethods are found to be mutually consistent. Marathon Races
Figure 1 demonstrates our methods applied to a singlemarathon run (subject M1) of group M. The color-codedvalue of the scale-dependent exponent α ( s ) is shownin the first row as a function of the binned heart rate(HR) [Fig. 1(a)] and also as a function of running time t [Fig. 1(c)]. Over the studied scales s from 5 to 5000heart beats, the scaling exponent α ( s ) exhibits complexbehavior that could not be adequately described by theconventional division into short- and long-range scalingexponents. We consider the HR-dependent shift to anti-correlated RRIs at the shortest scales s < ∼ s > ∼ the RRIs become mostly non-stationary ( α > , fractional-Brownian-motion-like be-havior). In contrast, a typical 24-hour RR-tachogram ofa healthy subject at rest usually displays /f or pinknoise on long time scales (or low frequencies < ∼ . Hz),corresponding to α = 1 , and larger values for α at theshortest scales or higher frequencies [1].The black curve in Fig. 1(a) corresponds to the con-ventional DFA exponent α over the scales from 4 to16 heart beats. It shows almost linear decrease to val-ues around / and below, when the DDFA exponent α ( s ) displays short-scale anticorrelated behavior. How-ever, this simpler estimate is not sufficient for uniquelydistinguishing the presence of anticorrelations from theirshift to slightly longer scales. To explore the scale depen-dence of the anticorrelations in more detail, we show inFig. 1(e) the probability density for the values of α ( s ) forsix different scales s from 5 to 20 heart beats as a func-tion of the binned HR. On all six scales, the probabilityis maximum for α < / with a HR dependent modula-tion and an absence of anticorrelations on the two largestscales s = 15 , for lower beat rates. In order to estimate the relevant time scales of thephysiological processes behind the observed anticorre-lated beat intervals, we have also performed a DPACFanalysis. The result is shown for lags between 1 and 20heart beats in Figs. 1(b) and 1(d). The PACF revealsdirect anticorrelations (negative values) after a time lagof 1 and 2 beats, starting at low exercise intensities, andadditional anticorrelations up to about 10 beats beyondHR of about 175 BPM, being consistent with the DDFAresults. The probability density of the DPACF valuesfor lags between 1 and 6 beats, shown in Fig. 1(f), con-firm dominant direct anticorrelations on the shortest timescales of 1 to 2 beats, and 4 beats for high exercise in-tensities (here HR > ∼ ).Subject M1 as the chosen example has the most promi-nent anticorrelations and particularly simple, almost lin-ear, trend in the HR over the whole marathon. The indi-vidual DDFA and DPACF results, similarly to Fig. 1, forall the subjects of group M are shown in SupplementaryFig. S5. The results share qualitative similarities acrossthe subjects, as they all exhibit short-scale suppressionof correlations and the appearance of anticorrelations asa function of the HR. However, some differences are alsoapparent, as only three subjects (M1, M3, and M7) showthe shift of the anticorrelations to elevated scales at thehighest exercise intensities. Some short-scale anticorre-lations, particularly for subject M6 and to some extentfor M4 and M5, also appear at elevated scales, but thesehappen at lower intensities and are likely different in ori-gin. Regardless, additional research is required to deter-mine the effect of individual strains relative to standardphysiological thresholds on the results.To further study the consistency of the results betweenthe different subjects, the aggregated DDFA (top) andDPACF (bottom) results for all members of group M areshown in Fig. 2 as a function of both the absolute (left)and relative (right) HR. The most notable features arethe high-intensity elevated-scale anticorrelations startingto appear at 87% and 95% relative HR, or at the abso-lute HR of 175 BPM (this congruence on the absolutescale is likely to be coincidental). In these ranges alsothe conventional DFA exponent α (black curve) dropsslightly below / , but its limitations are apparent, asit is based on linear regression over the scales of 4 to 16beats. The anticorrelations at lower intensities (approx-imately 155–175 BPM) appear to be more condensed onthe relative scale (roughly 78–87%, with a more concen-trated maximum between 80–85%), which is apparentboth on the DDFA results and in the more pronounceddip of the short-scale exponent α . There is also a band ofshort-scale suppressed correlations with a trend towardslonger scales at even lower relative HR ( ≈ τ = 1 , whereas on the absolute scale thetransition is more gradual. These lag τ = 1 anticorrela-tions appear consistently for all the subjects and become t , s ) 1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00Dynamic PACF ( t , )10 S ca l e s ( RR I) (a) (c) H ea r t r a t e H R ( B P M )
160 170 180Heart rate HR (BPM)15101520 L a g ( RR I) (b) t (s) (d) D e n s it y D e n s it y ( t , s ) (e) s = 5 s = 6 s = 8 s = 10 s = 15 s = 20
160 170 1801.00.50.00.51.0 ( t , ) (f) = 1
160 170 180 = 2
160 170 180 = 3
160 170 180 = 4
160 170 180 = 5
160 170 180 = 6 C onv e n ti on a l Heart rate HR (BPM)
FIG. 1. Beat-to-beat (RR) interval correlations for the Marathon race of subject M1. Note that the upper-left and upper-rightcolor bars refer to (a,c) and (b,d), respectively. (a)
Color-coded dynamic (DDFA-1) scaling exponent α ( t, s ) on different scales s (y-axis) as a function of binned HR (x-axis). Here α ( t, s ) is averaged over those dynamic segments, whose average HR fallswithin . BPM wide bins. The values for empty bins are linearly interpolated if the gap does not exceed . BPM. The blacksolid line shows the mean together with the standard deviation (thin lines) and the the standard error of the mean (thicklines, barely visible) of the conventional short-scale (4–16 RRIs) scaling exponent α . The exponent is computed in movingwindows of 50 RRIs in HR bins of BPM. (b)
Color-coded partial autocorrelation functions (DPACF-0) C ( t, τ ) with differentlags τ (y-axis) as a function of the binned HR. (c) Similar to (a) but as a function of time during the marathon race. Theinstantaneous heart rate is overlaid on the data. (d)
Similar to (b) but as a function of time. The values that do not pass thenon-zero significance test as described in the text are shown in white. (e)
Probability density histogram for α ( t, s ) for differentscales s as a function of the HR. (f ) Probability density histogram of the DPACF-0 for different lags τ as a function of the HR.The histograms in (e-f) consist of 31-by-31 bins, and the probability densities are separately normalized for each HR bin, sothat they better depict the distributions as a function of the HR instead of measuring the prevalence of different HR regions.Furthermore, the color bar is capped at the 99.5th percentile to avoid outliers dominating the color scale. stronger with increasing HR, and also appear at longerlags at the regions where the DDFA anticorrelations shiftto larger scales. Naturally, the aggregated results shouldbe interpreted with care as they represent an average re-sult over all the samples. Secondly, there is uncertaintyin the maximum HR values of the subjects. Neverthe-less, the results suggest that neither the absolute nor therelative scale is universal for different individuals. Freeform Training Runs
In order to study the correlations of RRIs over a widerange of exercise durations and intensities, we performthe same analysis for subjects in group T. It is instruc-tive to consider first a single exercise of one subject whichis shown in Fig. 3. It consists of six intervals of high-intensity running, each interval lasting about 160 secondswith the subsequent intervals reaching higher and higherintensities. As a function of exercise time t the DDFA ex-ponent α ( s ) [Fig. 3(c)] and the PACF [Fig. 3(d)] consis- t , s ) 1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00Dynamic PACF ( t , )10 S ca l e s ( RR I) (a)
150 160 170 180Heart rate HR (BPM)1510 L a g ( RR I) (b) (c) max (d) C onv e n ti on a l FIG. 2. Aggregate beat-to-beat (RR) interval correlations as a function of heart rate for all subjects of group M. (a)
Averagevalues for α ( t, s ) for each scale s (y-axis) and HR bin (x-axis). The solid line depicts the conventional short-range α . (b) Average values for C ( t, τ ) for each lag τ (y-axis) and HR bin (x-axis). In (a-b), the data is processed as in Fig. 1. (c-d) Similarto (a-b) but as a function of the relative HR. In (c-d), the data is processed as in Fig. 1, but with the distinction that therelative HR bin width is . , the interpolation threshold is . , and the bin width for the conventional short-scale exponent α in (c) is . . tently reveal strong anticorrelations of RR intervals thatdevelop rapidly after the start of the intense interval. Theshortest-scale anticorrelations span to longer and longerscales with increasing HR in the latter intervals. Theearlier lower-intensity intervals exhibit anticorrelations atelevated scales, separated by a band of correlations fromthe shortest scale anticorrelations. This behavior was al-ready suggested by some of the marathon data (M4, M5,and M6), and in a following analysis we will relate theseto a distinct band of anticorrelations appearing at mod-erate exercise intensity. At rest between the intervalsthe anticorrelations rapidly vanish. The DPACF showsstrong lag τ = 1 , anticorrelations, whose magnitudesare in accordance with the short-scale DDFA anticorrela-tions as observed in group M. The existence of patches ofanticorrelations over time lags up to 10 beats is also con-sistently observed with the elevated-scale DDFA anticor-relations. As a function of HR, the anticorrelated behav-ior develops rapidly after an intensity threshold ( ≈ α ( s ) and DPACF clearly show two distinct bands of anticor-related RRIs [see Fig. 4(a) and (b)]. The first band atmoderate exercise intensity ( ≈ > ∼ BPM)does not show a clear trend as a function of the HR,although there is tendency towards spanning to longerscales with increasing intensity. Notably, the anticorre-lations remain present even at the shortest scales andlag. This is in contrast to some of the marathon data,where the highest-intensity anticorrelations appear at el-evated scales. This could be due to the different natureof the exercises, as in the marathon races these anticorre-lations appear after prolonged exercise at high intensity,and changes in, e.g., body temperature or electrolyte bal-ance may influence the results. On the other hand, in thediscussion section we make an argument that this couldbe due to interactions with the stride frequency. Theconventional α indicates the suppression of correlationsthat is consistent with the DDFA anticorrelations whentaking into account its limitation to the scales of 4–16beats. The α is clearly insufficient to capture anticor-related behavior concentrated on thin bands of scales.Figures 4(c) and 4(d) show the probability density plotsof α ( s ) and PACF values for different scales s and lags τ , respectively. The existence of two regions with anti-correlated RRIs is clearly visible. They are separated bya region with positive correlations (or α > / ).The aggregated data for all the subjects of group T are t , s ) 1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00Dynamic PACF ( t , )10 S ca l e s ( RR I) (a) (c) H ea r t r a t e H R ( B P M )
130 140 150 160 170 180Heart rate HR (BPM)1510 L a g ( RR I) (b)
200 400 600 800 1,000 1,200Time t (s) (d) C onv e n ti on a l FIG. 3. Beat-to-beat RR interval correlations for one interval exercise from group T. (a)
Dynamic scaling exponents (DDFA-1) α ( t, s ) (colors) and the conventional short-range α (solid line) as a function of the binned HR. (b) DPACF-0 correlations C ( t, τ ) as a function of the binned HR. (c-d) As in (a-b) but as a function of time, and in (c) the HR value is overlaid on thedata. For details on the data processing, see the caption of Fig. 1. shown in Supplementary Fig. S4. Most subjects showcommon qualitative similarities in the form of two an-ticorrelated bands as described for T07. However, forsome subjects the split into the two anticorrelated re-gions is not that clear; particularly there is the lack ofcorrelated shortest scale behavior separating these tworegions. In the absence of the correlated bands the be-havior is remarkably simple; the higher the intensity, themore prominent the anticorrelations are in both magni-tude and scales covered. In addition to individual in-trinsic cardiac variability, a possible explanation couldbe different training practices and external conditions,as for example T05 shows behavior that is most simi-lar to the marathon data. Another explanation could behighly regular running motion, which could promote cor-relations induced by, e.g., muscle contractions and bloodpressure variations, which is an argument set forth in thenext section. It is also worth noting that T11 reportedproblems with the chest strap, and as a result his datahas unusually high amount of missed beats (up to 50%).Despite of this, two regions of suppressed correlationsare present that are consistent with the other subjects,highlighting the robustness of the methodology.We assessed the suitability of higher order detrendingfor our analysis and decided to employ DDFA-1 due tothe following reasons: (i) The qualitative behavior re-mains the same at the shortest scales, which is the mostinteresting region for dynamic exercise intensity analysis,(ii) the short-scale bias in DFA is larger and crossoverscales are shifted with higher order methods, (iii) higherorders of DDFA appear to require longer dynamic seg-ments for similar statistical accuracy and have increasedcomputational cost. Finally, we point out that it is important to check thereliability of our DDFA and DPACF methods with re-spect to trends. Hence, we have filtered the data of sub-ject T07 according to the condition that the standard de-viation of the HR within the dynamic segments is smallerthan the values for certain quantiles. We find that theobservation of the bands with anticorrelations is robustand independent of the choice of the quantile filter. Infact, the anticorrelations appear stronger when limitingto dynamic segments with less HR variation, as the av-eraging is not performed over segments with transientchanges in the HR that could lead to spurious correla-tions. The exact results of this analysis for six differentchoices of quantiles are shown in Supplementary Fig. S6.
DISCUSSION
It is important to understand the physiological mech-anism causing the observed anticorrelations. Due to thelack of time series for other physiological variables, wecan present below only simple arguments that we con-sider to be potentially relevant for explaining the ob-served dynamic correlations. First, we point out thatthere are three physiologically relevant time scales thatfall into the range over which the anticorrelations oc-cur: (i) the stride frequency which is typically around 85strides per leg and per minute [13], (ii) the respirationcycle which is typically three to five heart beats long,and (iii) arterial blood pressure fluctuations, i.e., the so-called Mayer waves, which result from an oscillation ofsympathetic vasomotor tone and is of the order of tenseconds [14]. t , s ) 1.0 0.5 0.0 0.5 1.0Dynamic PACF ( t , )110 120 130 140 150 160 170 180 190Heart rate HR (BPM)10 S ca l e s ( RR I) (a)
110 120 130 140 150 160 170 180 190Heart rate HR (BPM)15101520 L a g ( RR I) (b) D e n s it y D e n s it y ( t , s ) (c) s = 5 s = 6 s = 8 s = 10 s = 15 s = 20
120 140 160 1801.00.50.00.51.0 ( t , ) (d) = 1
120 140 160 180 = 2
120 140 160 180 = 3
120 140 160 180 = 4
120 140 160 180 = 5
120 140 160 180 = 6 C onv e n ti on a l Heart rate HR (BPM)
FIG. 4. Aggregate beat-to-beat (RR) interval correlations for all the exercises of one subject (T07) in group T. (a)
Averagevalues for α ( t, s ) for each scale s (y-axis) and binned HR (x-axis). The solid line shows the conventional short-range α . (b) Average values for C ( t, τ ) for each lag τ (y-axis) and binned HR (x-axis). (c) Probability density histogram for α ( t, s ) fordifferent scales s as a function of the HR. (d) Probability density histogram for C ( t, τ ) for different lags τ as a function of theHR. Note that the probability densities are separately normalized for each HR bin. For details on data processing, please seethe caption of Fig. 1. All three processes are cyclic and hence can induceperiodic modulations to the heart rate through hemody-namics. Such periodicities could result in anticorrelatedRRIs when observed at scales similar to the period mea-sured in heart beats. Furthermore, the overall heart ratevariability is reduced under exercise due to withdrawalof cardiac vagal tone and parasympathetic control[2–4],i.e., the local short-term RRIs are more regular with-out /f - or Brownian-like diffusion for extended peri-ods of time. Therefore, subtler patterns should becomemore discernible, as they are not masked by the complexfluctuations of a healthy heart under resting conditions.Similarly, the relative magnitudes of the modulating sig-nals could affect the scale-dependence of the anticorre-lations. If some of the effects is much stronger, it willmask higher frequency periodicities as they will appearcorrelated when superimposed on the stronger lower fre-quency oscillations. Additionally, when a periodic signalis sampled at discrete intervals, the result is a new signalwhose period depends on the sampling frequency. Thiseffect is manifested, e.g., when the influence of the blood pressure variations due to the stride frequency is sampledat each heart beat.This latter phenomenon could explain some of thequalitative differences in the dynamic correlations be-tween the subjects. For some subjects (particularly T08,but also T02, T03, T10, and T12) the RRIs show clearlydefined behavior under exercise, becoming short-scale an-ticorrelated at moderate intensity, with the magnitudeand the scale of the anticorrelations increasing in con-junction with intensity. In contrast other subjects exhibitmore complex RRI-correlations where the simple anticor-relations are interrupted by bands of decreased or alteredcorrelations at shorter scales (please see SupplementaryFigs. S4 and S5 for the individual RRI-correlation plotsas a function of the HR). These more correlated bands ap-pear at heart rates corresponding to sampling frequencieswhere typical stride frequencies would look correlated atthe shortest scales. If these bands arise from the stridefrequency, that could also explain the better congruenceon the absolute HR scale in Fig. 2, as there is generallyless variance in the stride frequencies than in the maxi-mum heart rates.These considerations would imply that the(anti)correlations arise from underlying universalcardiolocomotor mechanisms, but detailed response toexercise may depend on individual physiology, biome-chanics of running and training status [6]. Furthermore,the onset of the anticorrelations, their strength, andscales of appearance show individual variability. Study-ing the relationship of these variables to standardizedthresholds and markers in exercise physiology couldallow utilizing the dynamic correlations for monitoringthe exercise intensity in real-time without the knowledgeof parameters such as the maximal oxygen uptake(VO2max) or maximum heart rate. We are aware ofthe possibility that the universal emergence of anticor-relations at elevated heart rates is most likely affectedby other physiological factors beyond the ones discussedhere. Clearly, further research is required, but theapproach herein provides a promising avenue forward. CONCLUSIONS
Our main result is the discovery of multiscale anticor-relations in RR intervals during running exercises underreal-world conditions. The anticorrelations have a dy-namical structure that depends on the exercise intensityas measured by the heart rate. The characteristics of thedynamical structure are revealed by our methodology,in particular the dynamic detrended fluctuation analysisand dynamic partial autocorrelation functions, which weanticipate becoming useful tools in data analysis acrossvarious disciplines. While we have demonstrated the ca-pability to study the dynamical RRI correlations dur-ing varying real-world circumstances, a more systematicevaluation of the methodology is required to control forexercise conditions.The observed anticorrelations appear on short scales(a few beats) at low to moderate exercise intensities. Asthe intensity is increased, the anticorrelations increasein magnitude and span to longer scales (up to 20–30beats). This simplified picture is complicated by correla-tions arising potentially from interactions with the regu-lar running motion when the stride frequency is appropri-ately proportional to the heart beat. These correlationsmask the anticorrelated behavior on bands of increas-ing or decreasing scales at moderate and high exerciseintensities, respectively. At rest, e.g., between runningintervals, the anticorrelations rapidly vanish, and appearimmediately when the intensity is increased again. Thesechanges happen before the HR saturates at the level nec-essary to maintain the ongoing exercise intensity. Hence,our findings allude the possibility of quantifying the rel-ative exercise intensity by measuring the dynamic corre-lation exponent α ( t, s ) in real time during exercise.This report of our initial findings serves as a preludefor highlighting the potential of the dynamic correlationanalysis so that further advances could be pursued. It is highly desirable to develop a theoretical model for thecomplex dynamics of the cardiovascular feedback loopsduring high-intensity exercise load that can explain theobserved time scales for the anticorrelated RR intervals.Clearly, a more systematic study with subjects perform-ing specific exercise protocol should be performed to ver-ify our observations. Besides, a thorough validation andcalibration of our results with data collected during run-ning exercise in a physiology laboratory is a natural nextstep for our study to relate the changes in the dynamiccorrelations to standard exercise physiology models. Theinclusion of accelerometer data, from which the stridefrequencies could be derived, would facilitate the verifi-cation of the running modalities as a possible cause forthe bands of correlations present for some subjects.Such controlled and systematic studies are not onlynecessary to elucidate the speculative nature of the re-sults herein, but they are further motivated by poten-tially enabling the application of this methodology in ex-ercise physiology. We expect that the reported RR in-terval correlations are suitable to represent a dynamical“fingerprint” of the exercise-induced cardiovascular load.Hence, our methodology – which could be integrated withthe present devices on the market – has a potential to be-come a new tool in real-time exercise monitoring withoutprevious knowledge of maximal thresholds such as themaximum hearth rate and lactate or ventilatory thresh-olds. MATERIALS AND METHODSHeart Rate Data during Exercise
We study real-time correlations of RRIs during exer-cises of various intensities. All heart rate data for thisstudy have been collected during regular running train-ing and racing under real-world conditions, i.e., outsidethe laboratory. Two groups of data were used for our the-oretical analysis. The first group of data was recorded byhuman volunteers during their regular running trainingwith freely chosen intensity and volume (group T). Thestudy period was at least 4 weeks, and some subjectsprovided data over a longer period of time. We obtainedinstitutional approval and informed consent (COUHESexemption for the employed protocol has been grantedunder protocol no. 1711132002). The research was per-formed in accordance with the rules and regulations setby the participating universities. This group involved 12volunteers (5 female, 7 male, with an age span from 27to 65 years). Their performances span a wide range fromtop national level to recreational runners: the personalbests in 10 km range from 29 min 31 sec to 44 min 57sec, in marathon from 2 hours 43 min 20 sec to 4 hours26 min 3 sec.During exercises, heart rate (HR), RR intervals (RRI),running velocity and distance were recorded using aGarmin heart rate monitor HRM4-Run and a GPS watch(Forerunner 935, Garmin Inc., Olathe, KS, USA). A pre-vious study has investigated and validated the accuracyof this HRM [15]. The data were recorded by the GPSwatch in the Flexible and Interoperable Data Transfer(FIT) format [16] and subsequently uploaded by the sub-jects to a web service that we had launched for this study.The total number of exercise files analyzed per subject(samples) varied between 18 and 261, with total covereddistances from 150 to 1889 km.The second group of data was obtained by selectingrandomly the marathon races of 7 subjects from datauploaded to the Polar Flow web service[17] (group M).Within registration to Polar Flow, the subjects havegiven their consent for the use of their data for researchpurposes. The metadata were provided by the users ofthis web service (all male, with an age span from 28 to 53years, and marathon finishing times between 3 hours 30min and 4 hours 17 min). HR and RRIs were recordedfor this group of subjects with a Polar heart rate monitorH10 and a Pro Strap (Polar Electro Oy, Kempele, Fin-land). Recently, the RR signal quality of this HRM hasbeen shown to be excellent from low- to high-intensityactivities in comparison to a ECG Holter device [18].In both groups T and M, the subjects provided theirmaximum and resting heart rates. Summaries of all themetadata for the two groups are shown in SupplementaryTable S1.As ECG data is not available, we do not attempt toremove ectopic beats or other artifacts based on physi-ological criteria. Therefore we merely remove technicalartifacts, such as missed beats, that can be isolated withreasonable certainty. The details for this data prepro-cessing are provided in Supplementary Appendix A.
Conventional Methods
For comparison we apply ordinary detrended fluctua-tion analysis [9, 11] to the RRI time series. By comput-ing the root-mean-squared fluctuations F ( s ) around localtrends at multiple scales s , the method assesses power lawscaling relations F ( s ) ∝ s α characterized by the scalingexponent α . In the context of HRV, typically two expo-nents are determined, for short- ( α ) and long-scale ( α )correlations, respectively [10, 12]. We extract the conven-tional short-scale (4–16 RRIs) scaling exponents α [10]in segments consisting of 50 RRIs. We compute the fluc-tuation functions in maximally overlapping windows forenhanced statistical properties [19]. A summary of theDFA method is provided in Supplementary Appendix B.We also provide a helpful summary of partial autocor-relation functions in Supplementary Appendix C beforeintroducing their dynamic counterparts here. Dynamic Segmentation
The dynamic behavior of the time series can be studiedby performing the analysis in moving temporal segments.However, to guarantee sufficient statistical accuracy, thelength of these segments is dictated by the largest scale s (DFA) or the lag τ (PACF), resulting in diminished tem-poral resolution for small scales. Therefore, we proposea dynamic segmentation procedure, where the segmentlength is varied as a function of the scale or the lag:1. Choose a function for determining the segmentlengths (cid:96) ( s ) as a function of the scale s . Here weadopt a simple linear relationship (cid:96) ( s ) = as where a is a constant. Smaller values increase the tem-poral resolution but also the statistical noise. Thedynamic length factor a itself may also be variedfor different scales.2. For each scale divide the time series into segmentsof length (cid:96) ( s ) . The segments themselves may beoverlapping if desired for smoother results. Identifythe segments S s,t by their temporal indices t , whichmay be, e.g., the mean time within the segment orany other suitable quantity. Dynamic Detrended Fluctuation Analysis (DDFA)
The dynamic segmentation together with the maxi-mally overlapping windows in the DFA scheme enablesthe following procedure for dynamic DFA (DDFA):1. Perform the dynamic segmentation for each scale s . The value of a = 5 was found to be an accept-able value for the dynamic length factor, which isemployed in all of our DDFA calculations.2. Utilizing overlapping windows, compute the fluc-tuation function in each segment S s,t at scales { s − , s, s + 1 } . Denote the logarithmic fluctua-tion function at these scales by ˜ F t ( s − , ˜ F t ( s ) and ˜ F t ( s + 1) , respectively.3. In each segment, compute the dynamic scaling ex-ponent α ( t, s ) by the finite difference approxima-tion [20] α ( t, s ) ≈ (cid:2) h − ˜ F t ( s + 1) + (cid:0) h − h − (cid:1) ˜ F t ( s ) − h ˜ F t ( s − (cid:3) / (cid:2) h − h + ( h + + h − ) (cid:3) , (1)where h − = log( s ) − log( s − and h + =log( s + 1) − log( s ) are the logarithmic backwardand forward differences. Fluctuation functionscomputed with maximally overlapping windows areempirically found to be smooth enough to per-mit the direct application of the finite differencescheme.The performance of the method is numerically vali-dated by applying it to simulated time series with knownproperties. Supplementary Appendix D explains the de-tails for analytically obtaining the theoretically expectedscale-dependency of DFA scaling exponents for differentprocesses. In Supplementary Appendix E these theoreti-cal results are utilized for confirming the acceptable per-formance of the DDFA method. Dynamic Partial Autocorrelation Function (DPACF)
In order to obtain a local estimate of the partial au-tocorrelation function C ( τ ) we compute it using an ap-proach similar to that of the DDFA algorithm. The stepsof this approach can be summarized as follows:1. Perform dynamic segmentation for each lag τ . Thevalue of a = 10 was found to be an acceptable valuefor the dynamic length factor, which is utilized inall of our DPACF calculations.2. In each segment S τ,t , perform polynomial detrend-ing of order m .3. For each segment, compute C ( τ ) by, for exam-ple, solving the Yule-Walkers equations with theLevinson-Durbin recursive scheme [21]. Chooseeach time for the maximum lag the parameter forwhich we are estimating the partial autocorrelationfunction. Denote this dynamic PACF by C ( t, τ ) .Resorting to the central limit theorem, it is a knownresult that the partial autocorrelation function is approx-imately non-zero at 5% significance level if |C ( t, τ ) | < . / (cid:112) (cid:96) ( τ ) . The evaluation of this significance band isstatistically valid only if (cid:96) ( τ ) > ∼ and therefore if τ > for a = 10 .Notice that in DPACF the detrending is applied to theoriginal time series in contrast to the integrated seriesin DDFA. Therefore, the results would be expected tobe qualitatively similar when the DPACF detrending or-der is one smaller than the DDFA detrending order n .This explanation is complicated by, e.g., the removal oflinear correlations in PACF. However, the relationship m ≈ n − is supported by empirical observations.While both DDFA and DPACF measure dynamic cor-relations, it is important to realize the qualitative dif-ference between them. DDFA describes the collective behavior of all beats over the scale s , whereas DPACFconsiders the average behavior of individual beats sepa-rated by the lag τ (with the linear dependence from thepreceding lags removed) ACKNOWLEDGMENTS
We acknowledge insightful discussions with Ary Gold-berger. We also thank Jyrki Schroderus, Laura Kar-avirta, Jussi Peltonen and Arto Hautala for all the support and suggestions. The authors acknowledgePolar Electro Oy (Oulu, Finland) for providing exer-cise data and CSC – IT Center for Science, Finland,for computational resources. Additional support wasprovided by an EMERGENCE grant of CNRS, INP,and the ICoME2 Labex (ANR-11-LABX-0053) and theA*MIDEX projects (ANR-11-IDEX-0001-02) funded bythe French program “Investissements d’Avenir” managedby ANR.
Appendix A: Dataset details and preprocessing
We study two datasets of running exercises and racesunder real-world conditions. These consist of voluntaryrunning exercises at freely chosen intensities and sched-ules (group T) and marathon races (group M). Metadatafor these datasets are presented in Supplementary TableS1.Artifacts in the data are removed prior to the analysisby utilizing the following scheme:1. RR intervals below or above the threshold values RR min and RR max are removed.2. Compute the local median of the RR intervals RR med ( t ) in a moving window of length l med .3. Remove RR intervals that fall outside the range (1 ± c med )RR med ( t ) , where the threshold c med is aconstant.The values for the filtering parameters are listed in Sup-plementary Table S2. Acceptable performance of the fil-ter is manually inspected for all the samples in group Mand a representative subset of 17 samples from the groupT. The filter is particularly adept at removing too longintervals arising from missed beats, which is the mostcommon error in exercise data [22]. As ECG data is notavailable, we do not attempt to filter the data based onphysiological criteria, and merely remove technical arti-facts that can be isolated with reasonable certainty. Appendix B: Detrended Fluctuation Analysis (DFA)
Ever since its introduction in the study of correla-tions in DNA sequences [9], DFA for time series hasbeen widely employed across multiple disciplines such asphysics [23, 24], medicine [1, 10, 25], finance [26], andeven music [27, 28]. The DFA method has been exten-sively studied [11, 29–36] and it has been expanded toaccount for effects such as multifractality [37] and cross-correlations [38].We briefly summarize the conventional DFA algorithmwhich has been developed to detect correlations in non-stationary time series [9, 10]. First, for a time series X ( j ) Supplementary Table S1. Summary of the metadata for the two groups of subjects: Training runs of various durations andintensities (T) and marathon races (M). For group T shown are the age, gender, resting heart rate HR rest and maximum heartrate HR max (as reported by the subject), the personal best (PB) time for 10 km and marathon (time in hh:mm:ss) within 3years before this study, the number of exercise samples analyzed, the average heart rate HR avg of all samples, and the totaldistance and duration of all samples. For group M all subjects were male, and shown are their age, resting heart rate HR rest and maximum heart rate HR max (as reported by the subject), marathon finishing time, and average heart rate HR avg duringthe marathon. Group Tsubject age [y] gender HR rest HR max PB 10 km PB Marath. samples HR avg distance [km] duration [h]T01 48 m 40 192 00:36:49 02:50:49 186 151 1527 121.5T02 27 m 40 193 00:29:31 – 54 130 576 43.5T03 29 f 43 200 00:36:50 02:47:56 20 155 170 13.7T04 29 f 50 205 00:42:30 – 103 167 1076 92.8T05 39 f – 200 00:44:57 03:40:09 15 176 150 11.7T06 – m 40 192 00:35:22 02:43:25 26 153 666 49.3T07 33 m 42 214 00:33:56 02:43:20 261 154 1889 138.7T08 37 f 43 200 – 04:26:03 20 149 209 22.1T09 27 m 48 195 00:31:32 – 21 143 199 16.7T10 37 f 45 179 – – 18 154 215 18.1T11 65 m 47 170 00:43:10 03:13:29 26 134 287 28.2T12 47 m 48 182 – 03:09:49 53 129 1133 95.4Group Msubject age [y] HR rest HR max Finishing time [hh:mm:ss] HR avg
M1 53 56 185 03:40:09 172M2 50 46 174 03:36:47 149M3 28 58 198 04:17:04 171M4 34 50 195 04:12:36 155M5 43 55 200 03:49:03 162M6 39 55 194 04:19:36 155M7 50 50 200 03:30:33 174Supplementary Table S2. Data filtering parameters.Group RR min (ms) RR max (ms) l med (beats) c med M 250 600 15 . T 250 1000 11 . of length N a cumulative summation is performed, Y ( k ) = k (cid:88) j =1 ( X ( j ) − (cid:104) X (cid:105) ) , (B1)where the mean (cid:104) X (cid:105) of the time series is subtracted, butthat is not strictly necessary for DFA [11]. Convention-ally, the integrated time series of (B1) is divided intonon-overlapping windows of length s . In each window w ,a local trend is determined as the least-squares fit of alow order polynomial p s,w ( k ) to the data. (The methodis denoted by DFA- n if the degree of the detrending poly-nomial is n [11].) The fluctuations are measured as thevariance from the local trend p s,w ( k ) in each window: F s,w = s (cid:80) k ∈ w ( Y ( k ) − p s,w ( k )) . These squared fluctu-ations are averaged over the windows to yield the fluctu-ation function F ( s ) = (cid:10) F s,w (cid:11) / . (B2)Allowing the windows to overlap enhances the statistical properties of this estimate [19]. When this procedure isrepeated for different window sizes, or scales s , a power-law increase of the fluctuations with the window size maybe observed, i.e., F ( s ) ∼ s α . Here α is a scaling exponentthat can be considered as a generalization of the Hurstexponent H [39]. See Supplementary Appendix E belowfor more details.However, experimental time series rarely exhibit exactscaling over several scales. Many previous studies havefocused on finding a robust determination of the scalingregimes [40, 41], or on extracting a spectra of scaling ex-ponents α ( s ) [42–46]. These methods are based on thenotion that the spectra may be defined as the local slope of the logarithmic fluctuation function, α ( s ) = d [log F ( s )]d [log s ] . (B3)In the context of HRV, these methods generalize and ex-pand the conventional division into short (4–16 beats)and long-range (16–64 beats) scaling exponents. In prac-tice, the behavior may also change over time, either dueto external influences, or the process itself may compriseseveral distinct intrinsic modes. This paper develops amethodology that takes these temporal variations intoaccount in a consistent manner.Depending on the value of the exponent α , differentdegrees of correlations of the time series or its increments1 Supplementary Table S3. Meaning of values of the DFA scal-ing exponent α . The qualitative interpretation remains thesame for even higher exponents that become discernible withhigher-order DFA: For each integral interval the lower andupper halves correspond to originally anticorrelated or corre-lated increments, respectively.Scaling exponent Interpretation Stationarity < α < / anti-correlated stationary α = / white noise / < α < correlated α = 1 1 /f (pink) noise < α < / anti-correlated increments non-stationary,stationaryincrements α = 1 / Brownian noise / < α < correlated increments can be identified. The meaning of the different ranges for α are summarized in Supplementary Table S3.The scaling exponent α is related to other scaling ex-ponents in time series analysis. Scale invariance is alsoobserved in the Fourier domain as a function of the fre-quency f with a power spectral density that scales in thelow frequency limit as P ( f ) ∼ f − β . The exponent β isrelated to the DFA exponent by the scaling relation β =2 α − [31, 32]. In exercise physiology, the power spectrumof heart rate time series is a frequently employed tool toquantify the cardiological response to exercise. However,analyses in the frequency domain are potentially plaguedby non-stationarity. For stationary signals, the autocor-relation function C ( τ ) = (cid:104) X ( τ ) X ( τ + τ ) (cid:105) decays forlong lags τ with a power law ∼ τ − γ . Then the scalingrelation γ = 2 − α holds [33]. For more details on theDFA method and its relation to correlation functions, seeSupplementary Appendix D. Appendix C: Partial Autocorrelation Function(PACF)
It is instructive to supplement DFA analyses by a di-rect study of correlations at different time scales by com-puting the autocorrelation function C ( τ ) and the par-tial autocorrelation function C ( τ ) at lag τ . The latterhas been successfully used to identify the best autore-gressive (AR) process to fit a time series, using the factthat C ( τ ) = 0 for all τ > p for an AR model of order p [47]. The autocorrelation function C ( τ ) is dominatedby trends in the data, suggesting apparent correlations.On the contrary, C ( τ ) is less affected by trends due tothe subtraction of the linear dependence on intermediatelags from the autocorrelation function. If the time seriescontains oscillations, C ( τ ) shows a periodic pattern witha frequency that modulates the data. On the contrary, C ( τ ) shows anticorrelations. The lags τ for which C ( τ ) assumes negative values can be related to the periodic-ity of the oscillations. More specifically, for a time series X ( τ ) the partial autocorrelation function is given by C ( τ ) = (cid:68) [ X ( τ ) − ˆ X τ τ ( τ )][ X ( τ + τ ) − ˆ X τ τ ( τ + τ )] (cid:69) (C1)where ˆ X τ τ ( τ (cid:48) ) is the best linear predictor, determinedby ˆ X τ τ ( τ (cid:48) ) = c + τ − (cid:88) i =1 c i X ( τ + i ) , (C2)where the coefficients c i are determined by the conditions c + τ − (cid:88) i =1 c i (cid:104) X ( τ + i ) (cid:105) = (cid:104) X ( τ (cid:48) ) (cid:105) , (C3) c (cid:104) X ( j ) (cid:105) + τ − (cid:88) i =1 c i (cid:104) X ( τ + i ) X ( j ) (cid:105) = (cid:104) X ( τ (cid:48) ) X ( j ) (cid:105) (C4)for j = τ + 1 , . . . , τ + τ − . The function C ( τ ) canbe computed practically from the Yule-Walker equations[47]. The relation between the AR fits and C ( τ ) is usefulin the performed analysis. Indeed, it has been shownthat a signal with non-trivial periodic behavior can be,for short time scales, successfully fitted by an AR process,and its dominant frequency of oscillation can be extractedfrom the estimated coefficients. These findings have beenalso confirmed by DFA [48]. Appendix D: Additional remarks on DFA for thevalidation of DDFA
In this section we provide some known theoretical re-sults for the conventional DFA algorithm to supportthe validation procedure presented in the following sec-tion. In DFA, the range of detectable exponents is de-termined by the degree of detrending, n , and is givenby ≤ α ≤ n + 1 [32]. While the existence of values α > may be criticized as a failure of the detrend-ing procedure [49], they may also be understood as anadvantage of the method for allowing meaningful quan-tification of non-stationary processes [35, 36]. The de-trending may be considered successful if it achieves thestatistical equivalence over the DFA windows, so thatthe fluctuation function F ( s ) does not depend on thewindow [35, 36]. This condition is fulfilled for DFA- n with time series exhibiting polynomial trends of de-gree n − . In general, for two uncorrelated signals X A ( t ) , X B ( t ) (random processes or trends), a superposi-tion principle holds, stating that the squared fluctuationfunction of the sum X A + B ( t ) = X A ( t ) + X B ( t ) is givenby F A + B ( s ) = F A ( s ) + F B ( s ) [29].For stationary processes and for non-stationary pro-cesses with stationary increments the fluctuation func-tion F s does not depend on the window and may be an-alytically computed [33, 34]. Its squared value is deter-mined as the weighted sum of the autocovariance function2 ˆ C ( τ ) = (cid:104) X ( τ ) X ( τ + τ ) (cid:105) − (cid:104) X (cid:105) in the former case, andthat of the variogram S ( τ ) = (cid:68) [ X ( τ + τ ) − X ( τ )] (cid:69) inthe latter case, F s = s − (cid:88) j = − s +1 G ( j, s ) ˆ C ( j ) (D1) F s = − s − (cid:88) j =1 G ( j, s ) S ( j ) (D2)with the weight function G ( j, s ) given by [50] G ( j, s ) = 1 s s −| j | (cid:88) k =1 a k,k + | j | , (D3)where a k,k (cid:48) are the elements of the matrix A = D (cid:62) (cid:104) I − B (cid:62) (cid:0) BB (cid:62) (cid:1) − B (cid:105) D , (D4)where the elements d i,j of the matrix D are unity if i ≥ j and zero otherwise [33, 34]. The effect of detrending isincorporated into the so-called design matrix B of leastsquares regression, which for DFA-1 is given by [34] B = (cid:20) · · ·
11 2 · · · s (cid:21) . (D5)This matrix A describes an operator for constructing thesquared fluctuations F s = X (cid:62) s,w AX s,w from the valuesof the time series X s,w in window w at the scale s [34].The autocovariance function for fractional Gaussiannoise (fGn) and the variogram for fractional Brownianmotion (fBm) with Hurst parameter H are known to be ˆ C fGn H ( τ ) = σ (cid:16) | τ + 1 | H − | τ | H + | τ − | H (cid:17) , (D6) S fBm H ( τ ) = σ | τ | H , (D7)where σ corresponds to the variance of ordinary Gaus-sian noise [51]. From these correlations and by usingEqs. D1 and D2, the theoretical fluctuation function maybe computed for these processes. The spectrum of thescaling exponent α ( s ) is then obtained from (B3). Thesetheoretical results may be utilized for studying the be-havior of the DFA method as a function of the scale s .For example, the deviation between the DFA estimatefor α ( s ) and the asymptotic large scale exponent α isvisualized in Supplementary Fig. S1 for fGn and fBm.The well-known overestimation of the scaling exponentat the shortest scales is clearly visible, and it is mostpronounced in the anticorrelated region with α < / .Around the asymptotic value of α = 1 there is an abruptqualitative change as the scaling exponent is suddenlyunderestimated for an extended range of scales. This hasbeen observed previously [34].In general, the short scale behavior depends on thedetails of the underlying process, and hence can be dif-ferent for other processes such as an autoregressive model Scale s A s y m p t o ti c o f f G n a nd f B m -0.5-0.2-0.1-0.05-0.02-0.010.010.020.050.10.20.5 D e v i a ti on fr o m a s y m p t o ti c Supplementary Figure S1. Theoretical deviation of the DFAestimate α ( s ) from asymptotic scaling exponent α for fGnand fBm as a function of the scale s . The deviation is definedas the theoretical scale-dependent exponent α ( s ) minus theasymptotic scaling exponent α . Note the quasi-logarithmicscale for the deviation. AR(p). For more details on the scale dependence of thedeviation between asymptotic α and α ( s ) , please see alsoRef. [52]. Appendix E: Numerical Validation of DDFA
Fractional Brownian motion (fBm) and its increments,fractional Gaussian noise (fGn), are commonly utilizedfor benchmarking DFA (see, e.g., Ref. [53] and referencestherein). These processes are characterized by the Hurstparameter ≤ H < , and exhibit long-range correla-tions with the asymptotic (large-scale) scaling exponents α = H for fGn and α = H +1 for fBm [34, 51]. Deviationsfrom the asymptotic behavior occur at shorter scales dueto the finite length of the samples and the intrinsic biasin DFA due to the detrending. It is, however, possibleto compute the exact theoretical scale-dependent scalingexponent α ( s ) for these processes, as described in Sup-plementary Appendix D.We validate the dynamic DFA (DDFA) method by ap-plying it to simulated fGn and fBm, and comparing theresults to the theoretically expected values. We utilizethe Davies–Harte method, which is an efficient methodfor simulating these processes with their exact covariancestructure [54][55]. We generate samples of fGn andfBm of length for each value of the Hurst parame-ter H . From these simulated time series, we computethe dynamic scaling exponent α ( t, s ) in non-overlappingdynamic segments with various dynamic segment lengthfactors a . The mean difference between the DDFA expo-nent α ( t, s ) and the theoretically expected DFA exponent α ( s ) is illustrated in Supplementary Fig. S2. The generaltrend is that the limited sample size results in underesti-mation of the scaling exponents: the shorter the dynamicsegments, the greater the bias. Similarly, processes with3 -0.125-0.1-0.075-0.05-0.0250.0250.05 B i a s r e l a ti v e t o t h e o r e ti ca l ( s )
105 20 400.00.51.01.52.0 A s y m p t o ti c o f f G n a nd f B m (a)
105 20 40 (b)
105 20 40 (c)
105 20 40 (d)
Scale s Supplementary Figure S2. Bias of the dynamic detrendedfluctuation analysis (DDFA) estimate of α ( t, s ) relative toits theoretically expected value for fractional Gaussian noise(fGn) and fractional Brownian motion (fBm). Characterizedby the Hurst parameter ≤ H < , the asymptotic α forthese processes is α = H and α = H + 1 respectively. Detailsfor computing the theoretical scale-dependent DFA exponent α ( s ) are provided in Supplementary Appendix D. The biasis defined as the observed α ( t, s ) minus its theoretically ex-pected value α ( s ) . The dynamic segment length factors a are4, 5, 7, and 10 in (a-d), respectively. larger asymptotic α suffer from larger bias, with a weakdiscontinuity at α = 1 (when the process changes fromfGn to fBm). This tendency can be understood as aris-ing from the increased abundance and length of streaksin more correlated time series. This effect is not fullycaptured by the relatively short segments. However, forthe shortest scales (5 and 6), particularly in the anti-correlated region, the exponent is slightly overestimatedinstead. We also observe that the contour lines for thebias in α ( t, s ) have nearly converged to a constant valueof the asymptotic α already at the scale s = 40 .The standard deviation of the DDFA estimation of theexponent α is shown in Supplementary Fig. S3 for seg-ment lengths 4, 5, 7, and 10 in (a-d), respectively. Thedeviation consistently decreases with the increasing seg-ment length. On the other hand, the deviation reducedas the DFA-1 exponents approach the limits α = 0 and α = 2 , since at these boundaries there is less room forvariations. In particular, as α is reduced from / to-wards zero (corresponding to the anticorrelated regime)the deviations are strongly reduced. The local reductionin the deviation just above α = 1 is due to the most anti-correlated increments in this regime (see SupplementaryTable S3).Generally, the bias and the standard deviation of theDDFA method are found to be acceptable for our pur-poses, especially in view of the fact that we have par-ticular interest in the anticorrelated region as underlinedbelow in the results. All our DDFA computations are per-formed with the dynamic segment length factor a = 5 .This was found to be a good compromise between the S t a nd a r d d e v i a ti on o f ( t , s )
105 20 400.00.51.01.52.0 A s y m p t o ti c o f f G n a nd f B m (a)
105 20 40 (b)
105 20 40 (c)
105 20 40 (d)
Scale s Supplementary Figure S3. Standard deviation of the dy-namic detrended fluctuation analysis (DDFA) estimate of α ( t, s ) for fGn and fBm as a function of the scale s . Charac-terized by the Hurst parameter ≤ H < , the asymptotic α for these processes is α = H and α = H + 1 respectively. Thedynamic segment length factors a are 4, 5, 7, and 10 in (a-d),respectively. accuracy of the DDFA method and the dynamical reso-lution requiring a sufficiently small segment size. Appendix F: Additional Heartbeat Correlation Plots
Here we present beat-to-beat (RR) interval (RRI) cor-relations for all the subjects in the study. In Supplemen-tary Fig. S4 we illustrate the average RRI correlationresults as a function of the heart rate aggregated over allthe runs for each subject of Group T. The relative heartrate is utilized to better facilitate the comparison be-tween different individuals. Similar correlation plots forthe marathons of Group M are shown in SupplementaryFig. S5, along with the correlation landscapes as a func-tion of time during the marathon races. Additionally, weestablish the consistency of the anticorrelated bands inthe presence of possible trends in Supplementary Fig. S6.We demonstrate this by limiting the analysis to subsets ofdata where the heart rate within the dynamic segmentsexhibits subsequently lower and lower standard devia-tion. The analysis is performed for subject T07, who hasthe most data.Each correlation plot consists of pairs of color-codedDDFA (upper panels) and DPACF (lower panels) results.Plots as a function of the heart rate (HR) are based ondata that is averaged of dynamic segments whose averageHR falls within bins with widths of . BPM or . forthe absolute and relative HR, respectively. The valuesfor empty bins are linearly interpolated if the gap doesnot exceed . BPM (absolute) or . (relative). TheDDFA plots (as a function of the HR) also display theconventional short-scale (4–16 RRIs) scaling exponents α by a semi-transparent black line with error bars (thinbars: standard deviation, thick bars: standard error of4the mean, barely visible). The exponent is computed inmoving windows of 50 RRIs in HR bins of BPM (abso-lute) or . (relative). The DDFA plots displaying thecorrelation landscapes for single runs show the instanta-neous HR with the semi-transparent black line instead,and in the corresponding single run DPACF plots thevalues that do not pass the non-zero significance test areshown in white.5 t , s ) 1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00Dynamic PACF ( t , )10 S ca l e s ( RR I) T01 L a g ( RR I) T02 C onv e n ti on a l S ca l e s ( RR I) T03 L a g ( RR I) T04 C onv e n ti on a l S ca l e s ( RR I) T05 L a g ( RR I) T06 C onv e n ti on a l S ca l e s ( RR I) T07 L a g ( RR I) T08 C onv e n ti on a l S ca l e s ( RR I) T09 L a g ( RR I) T10 C onv e n ti on a l S ca l e s ( RR I) T11
90 100 110 120 130 140 150 160 170 180 190Heart rate HR (BPM)1510 L a g ( RR I) T12
90 100 110 120 130 140 150 160 170 180 190Heart rate HR (BPM) 0.00.51.01.52.0 C onv e n ti on a l Supplementary Figure S4. Aggregate RRI correlation results for each subject of Group T. For each subject the average DDFA-1scaling exponents α ( t, s ) (upper panels) and DPACF-0 correlations C ( t, τ ) (lower panels) as a function of binned relative heartrate. For a detailed explanation about how the data is computed, please see Supplementary Appendix F. t , s ) 1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00Dynamic PACF ( t , )10 S ca l e s ( RR I) M01 H ea r t r a t e H R ( B P M ) L a g ( RR I) C onv e n ti on a l S ca l e s ( RR I) M02 H ea r t r a t e H R ( B P M ) L a g ( RR I) C onv e n ti on a l S ca l e s ( RR I) M03 H ea r t r a t e H R ( B P M ) L a g ( RR I) C onv e n ti on a l S ca l e s ( RR I) M04 H ea r t r a t e H R ( B P M ) L a g ( RR I) C onv e n ti on a l S ca l e s ( RR I) M05 H ea r t r a t e H R ( B P M ) L a g ( RR I) C onv e n ti on a l S ca l e s ( RR I) M06 H ea r t r a t e H R ( B P M ) L a g ( RR I) C onv e n ti on a l S ca l e s ( RR I) M07 H ea r t r a t e H R ( B P M )
140 150 160 170 180 190Heart rate HR (BPM)1510 L a g ( RR I) t (s)0.00.51.01.52.0 C onv e n ti on a l Supplementary Figure S5. Overview of all the marathons in Group M. Left: Average RRI correlation results as a function ofbinned relative heart rate for each subject of Group M. Right: RRI correlation landscapes of the marathon runs. For eachsubject the DDFA-1 scaling exponents α ( t, s ) (upper panels) and DPACF-0 correlations C ( t, τ ) (lower panels) are shown. Fora detailed explanation about how the data is computed, please see Supplementary Appendix F. t , s ) 1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00Dynamic PACF ( t , )120140160180 H ea r t r a t e , m ea n DFA
DPACF Scale s DDFA Q : 1/2 1/4 1/8 1/16 1/32 2.55.07.510.012.515.0 H ea r t r a t e , S . D . DFA
DPACF Scale s DDFA S ca l e s ( RR I) Q = 1 L a g ( RR I) Q = 1/2 C onv e n ti on a l S ca l e s ( RR I) Q = 1/4 L a g ( RR I) Q = 1/8 C onv e n ti on a l S ca l e s ( RR I) Q = 1/16
90 100 110 120 130 140 150 160 170 180 190Heart rate HR (BPM)1510 L a g ( RR I) Q = 1/32 C onv e n ti on a l
90 100 110 120 130 140 150 160 170 180 190Heart rate HR (BPM)
Supplementary Figure S6. Consistency of the results with respect to trends. Top row: different quantiles Q of the mean andstandard deviation of the heart rate within the dynamic segments for all the data of subject T07. For DPACF-0 and DDFA-1the dynamic segment length factor a has values of 10 and 5, respectively, and for conventional DFA the short-range (4–16 beats)scaling exponent is computed in moving windows consisting of 50 RRIs. Lower panels: the average RRI correlation results asa function of the heart rate when the data is limited to dynamic segments with the heart rate standard deviation less thanthe value for the specified quantiles Q . For a detailed explanation about how the data is computed, please see SupplementaryAppendix F. [1] A. L. Goldberger, L. A. N. Amaral, J. M. Hausdorff, P. C.Ivanov, C.-K. Peng, and H. E. Stanley, Fractal dynamicsin physiology: Alterations with disease and aging, Pro-ceedings of the National Academy of Sciences , 2466(2002).[2] R. Karasik, N. Sapir, Y. Ashkenazy, P. C. Ivanov, I. Dvir,P. Lavie, and S. Havlin, Correlation differences in heart-beat fluctuations during rest and exercise, Phys. Rev. E , 062902 (2002).[3] A. J. Hautala, T. H. Mäkikallio, T. Seppänen,H. V. Huikuri, and M. P. Tulppo, Short-term cor-relation properties of R–R interval dynamics atdifferent exercise intensity levels, Clinical Phys-iology and Functional Imaging , 215 (2003),https://onlinelibrary.wiley.com/doi/pdf/10.1046/j.1475-097X.2003.00499.x.[4] T. Gronwald and O. Hoos, Correlation proper-ties of heart rate variability during enduranceexercise: A systematic review, Annals of Non-invasive Electrocardiology , e12697 (2020),https://onlinelibrary.wiley.com/doi/pdf/10.1111/anec.12697.[5] M. Buchheit, S. R., and G. P. Millet, Heart-rate deflec-tion point and the second heart-rate variability thresholdduring running exercise in trained boys, Pediatric Exer-cise Science , 192 (2007).[6] R. Di Michelle, G. Gatta, A. Di Leo, M. Cortesi, F. An-dina, E. Tam, M. Da Boit, and F. Merni, Estimationof the anaerobic threshold from heart rate variability inan incremental swimming test, Journal of Strength andConditioning Research , 3059 (2012).[7] Y. Yamamoto, R. L. Hughson, and J. C. Peterson, Au-tonomic control of heart rate during exercise studiedby heart rate variability spectral analysis, Journal ofApplied Physiology , 1136 (1991), pMID: 1757310,https://doi.org/10.1152/jappl.1991.71.3.1136.[8] G. Billman, The LF/HF ratio does not accurately mea-sure cardiac sympatho-vagal balance, Frontiers in Phys-iology , 26 (2013).[9] C.-K. Peng, S. V. Buldyrev, S. Havlin, M. Simons, H. E.Stanley, and A. L. Goldberger, Mosaic organization ofDNA nucleotides, Physical Review E , 1685 (1994).[10] C. Peng, S. Havlin, H. E. Stanley, and A. L. Goldberger,Quantification of scaling exponents and crossover phe-nomena in nonstationary heartbeat time series, Chaos:An Interdisciplinary Journal of Nonlinear Science , 82(1995).[11] J. W. Kantelhardt, E. Koscielny-Bunde, H. H. Rego,S. Havlin, and A. Bunde, Detecting long-range correla-tions with detrended fluctuation analysis, Physica A: Sta-tistical Mechanics and its Applications , 441 (2001).[12] F. Shaffer and J. P. Ginsberg, An overview of heart ratevariability metrics and norms, Frontiers in Public Health , 258 (2017).[13] D. E. Lieberman, A. G. Warrener, J. Wang, and E. R.Castillo, Effects of stride frequency and foor position atlanding on braking force, hip torque, impact peak forceand the metabolic cost of running in humans, J. Exp.Biol. , 3406 (2015).[14] C. Julien, The enigma of Mayer waves: Facts and models,Cardiovascular Research , 12 (2006).[15] J. Cassirame, R. Vanhaesebrouck, S. Chevrolat, and L. Mourot, Accuracy of the Garmin 920 XT HRM toperform HRV analysis, Austr. Phys. Eng. Sci. in Med. , 831 (2017).[16] Garmin Canada, The flexible and interoperable datatransfer (FIT) protocol software development kit, (2020).[17] Polar Electro, Polar Flow web service, https://flow.polar.com/ (2020).[18] R. Gilgen-Ammann, T. Schweizer, and T. Wyss, RR in-terval signal quality of a heart rate monitor and an ECGHolter at rest and during exercise, Euro. J. App. Phys. , 1525 (2019).[19] K. Kiyono and Y. Tsujimoto, Nonlinear filtering proper-ties of detrended fluctuation analysis, Physica A: Statis-tical Mechanics and its Applications , 807 (2016).[20] B. Fornberg, Generation of finite difference formulas onarbitrarily spaced grids, Mathematics of Computation , 699 (1988).[21] J. Durbin, The fitting of time-series models, Revue del’Institut International de Statistique / Review of the In-ternational Statistical Institute , 233 (1960).[22] D. A. Giles and N. Draper, Heart rate variability duringexercise: A comparison of artefact correction methods,The Journal of Strength & Conditioning Research ,10.1519/JSC.0000000000001800 (2018).[23] M. S. Santhanam, J. N. Bandyopadhyay, and D. Angom,Quantum spectrum as a time series: Fluctuation mea-sures, Phys. Rev. E , 015201(R) (2006).[24] V. Kotimäki, E. Räsänen, H. Hennig, and E. J. Heller,Fractal dynamics in chaotic quantum transport, Phys.Rev. E , 022913 (2013).[25] J. Kim, D. Shah, I. Potapov, J. Latukka, K. Aalto-Setälä,and E. Räsänen, Scaling and correlation properties of RRand QT intervals at the cellular level, Scientific Reports , 3651 (2019).[26] N. Vandewalle and M. Ausloos, Coherent and randomsequences in financial fluctuations, Physica A: StatisticalMechanics and its Applications , 454 (1997).[27] H. Hennig, R. Fleischmann, A. Fredebohm, Y. Hag-mayer, J. Nagler, A. Witt, F. Theis, and T. Geisel, Thenature and perception of fluctuations in human musicalrhythms, PloS one , e26457 (2011).[28] E. Räsänen, O. Pulkkinen, T. Virtanen, M. Zollner, andH. Hennig, Fluctuations of hi-hat timing and dynamicsin a virtuoso drum track of a popular music recording,PLOS ONE , 1 (2015).[29] K. Hu, P. C. Ivanov, Z. Chen, P. Carpena, and H. E.Stanley, Effect of trends on detrended fluctuation analy-sis, Phys. Rev. E , 011114 (2001).[30] Z. Chen, P. C. Ivanov, K. Hu, and H. E. Stanley, Ef-fect of nonstationarities on detrended fluctuation analy-sis, Phys. Rev. E , 041107 (2002).[31] C. Heneghan and G. McDarby, Establishing the relationbetween detrended fluctuation analysis and power spec-tral density analysis for stochastic processes, Phys. Rev.E , 6103 (2000).[32] K. Kiyono, Establishing a direct connection between de-trended fluctuation analysis and fourier analysis, Phys.Rev. E , 042925 (2015).[33] M. Höll and H. Kantz, The relationship between thedetrendend fluctuation analysis and the autocorrelation function of a signal, The European Physical Journal B , 327 (2015).[34] O. Løvsletten, Consistency of detrended fluctuation anal-ysis, Phys. Rev. E , 012141 (2017).[35] M. Höll, H. Kantz, and Y. Zhou, Detrended fluctuationanalysis and the difference between external drifts andintrinsic diffusionlike nonstationarity, Phys. Rev. E ,042201 (2016).[36] M. Höll, K. Kiyono, and H. Kantz, Theoretical founda-tion of detrending methods for fluctuation analysis suchas detrended fluctuation analysis and detrending movingaverage, Phys. Rev. E , 033305 (2019).[37] J. W. Kantelhardt, S. A. Zschiegner, E. Koscielny-Bunde,S. Havlin, A. Bunde, and H. E. Stanley, Multifractal de-trended fluctuation analysis of nonstationary time series,Physica A: Statistical Mechanics and its Applications , 87 (2002).[38] B. Podobnik and H. E. Stanley, Detrended cross-correlation analysis: A new method for analyzing twononstationary time series, Phys. Rev. Lett. , 084102(2008).[39] Originally it was the exponent in Hurst’s R/S analysis[56], and Mandelbrot and van Ness used it as a parame-ter for defining fGn and fBm [51]. For these processes itcan be related to the (asymptotic) DFA scaling exponent(fGn: α = H , fBm: α = H + 1 ). In the < α < rangethe relationship (between R/S and DFA exponents) isapproximate, and for fGn/fBm holds asymptotically forlarge scales. For α > the interpretation α = H + 1 isessentially due to the definition of fBm. It also could beargued to hold, at least approximately, for processes thatare defined as the cumulative sum of another process,similarly as the (discrete) fGn/fBm.[40] E. Ge and Y. Leung, Detection of crossover time scalesin multifractal detrended fluctuation analysis, Journal ofGeographical Systems , 115 (2013).[41] A. Habib, J. P. Sorensen, J. P. Bloomfield, K. Muchan,A. J. Newell, and A. P. Butler, Temporal scaling phenom-ena in groundwater-floodplain systems using robust de-trended fluctuation analysis, Journal of Hydrology ,715 (2017).[42] G. M. Viswanathan, C.-K. Peng, H. E. Stanley, and A. L.Goldberger, Deviations from uniform power law scalingin nonstationary time series, Phys. Rev. E , 845 (1997).[43] J. C. Echeverría, M. S. Woolfson, J. A. Crowe, B. R.Hayes-Gill, G. D. H. Croaker, and H. Vyas, Interpreta-tion of heart rate variability via detrended fluctuation analysis and αβ filter, Chaos: An Interdisciplinary Jour-nal of Nonlinear Science , 467 (2003).[44] P. Castiglioni, G. Parati, A. Civijian, L. Quintin, andM. D. Rienzo, Local scale exponents of blood pressureand heart rate variability by detrended fluctuation anal-ysis: Effects of posture, exercise, and aging, IEEE Trans-actions on Biomedical Engineering , 675 (2009).[45] J. Xia, P. Shang, and J. Wang, Estimation of local scaleexponents for heartbeat time series based on DFA, Non-linear Dynamics , 1183 (2013).[46] M. Molkkari and E. Räsänen, Robust estimation of thescaling exponent in detrended fluctuation analysis of beatrate variability, in Computing in Cardiology (2018).[47] G. E. Box, G. M. Jenkins, and G. C. Reinsel,
Time SeriesAnalysis, Forecasting and Control , 4th ed., Wiley Seriesin Probability and Statistics (Wiley, 2008).[48] P. G. Meyer and H. Kantz, Inferring characteristictimescales from the effect of autoregressive dynamics ondetrended fluctuation analysis, New Journal of Physics , 033022 (2019).[49] R. M. Bryce and K. B. Sprague, Revisiting detrendedfluctuation analysis, Scientific reports , 315 (2012).[50] We adapt the notation from Ref. [34] with the followingchanges: The factor s − is included in the weight func-tion G ( j, s ) that is symmetrically extended for negative j by G ( − j, s ) = G ( j, s ) . This allows for a more succinctnotation for Eqs. D1 and D2.[51] B. Mandelbrot and J. Van Ness, Fractional brownian mo-tions, fractional noises and applications, SIAM Review , 422 (1968).[52] M. Molkkari, Advanced Methods in Detrended Fluctua-tion Analysis with Applications in Computational Cardi-ology , Master’s thesis, Tampere University (2019).[53] Y.-H. Shao, G.-F. Gu, Z.-Q. Jiang, W.-X. Zhou, andD. Sornette, Comparing the performance of FA, DFA andDMA using different synthetic long-range correlated timeseries, Scientific Reports , 835 (2012).[54] R. B. Davies and D. S. Harte, Tests for Hurst effect,Biometrika , 95 (1987).[55] The accuracy of the simulated time series is establishedby comparing their joint fluctuation functions (the meanin (B2) is taken over the squared fluctuations F s,w of allrealizations of the simulated time series) to the theoreti-cal fluctuation functions of fGn and fBm.[56] H. E. Hurst, The problem of long-term storage in reser-voirs, International Association of Scientific Hydrology.Bulletin1