Heart Rate Dynamics During Acute Recovery From Maximal Aerobic Exercise in Young Adults

Introduction Resting heart rate (HRrest), heart rate variability (HRV), and HR recovery (HRR) from exercise provide valuable information about cardiac autonomic control. RR-intervals during acute recovery from exercise (RRrec) are commonly excluded from HRV analyses due to issues of non-stationarity. However, the variability and complexity within these trends may provide valuable information about changes in HR dynamics. Purpose Assess the complexity of RRrec and determine what physiologic and demographic information are associated with differences in these indices in young adults. Methods RR-intervals were collected throughout maximal treadmill exercise and recovery in young adults (n = 92). The first 5 min of RRrec were (1) analyzed with previously reported methods that use 3-interval lengths for comparison and (2) detrended using both differencing(diff) and polynomial regression(res). The standard deviation of the normal interval (SDNN), root mean square of successive differences (rMSSD), root mean square (RMS) of the residual of regression, and sample entropy (SampEn) were calculated. Repeated measures analysis of covariance (ANCOVA) tested for differences in these indices for each of the methodological approaches, controlling for race, body fat, peak oxygen uptake (VO2peak), and resting HR (HRrest). Statistical significance was set at p < 0.05. Results VO2peak and HRrest were significantly correlated with traditional measures of HRR and the variability surrounding RRrec. SampEndiff and SampEnres were correlated with VO2peak but not HRrest or HRR. The residual-method provided a significantly (p = 0.04) lower mean standard error (MSE) (0.064 ± 0.042) compared to the differencing-method (0.100 ± 0.033). Conclusions Complexity analysis of RRrec provides unique information about cardiac autonomic regulation immediately following the cessation of exercise when compared to traditional measures of HRR and both HRrest and VO2peak influence these results.


INTRODUCTION
Historically, the response of heart rate (HR) and heart rate recovery (HRR) to-or following-a perturbation, has been commonly utilized within research and clinical settings as a non-invasive physiological measure of cardiovascular regulation. Several studies have shown that impairment in autonomic function-represented by delayed HRR-indicates adverse cardiovascular outcomes in otherwise healthy individuals (Qiu et al., 2017;Lachman et al., 2018). While changes in HR and traditional indices of HRR can provide important information about differences in cardiac autonomic control, heart rate variability (HRV) has been shown to offer a more sensitive measure of cardiac autonomic regulation at rest (Task-Force, 1996) and during exercise (Karapetian et al., 2008) compared to a mean HR value. As measures of HR, HRR, and HRV have clear associations with cardiovascular health, a better understanding of RR intervals following various exercises and exercise intensities during the acute recovery phase (RR rec ) is needed within the literature. Doing so, provides implications for clinical, research, and performance-related exercise testing and monitoring. Michael et al. (2017) highlight and discuss the usefulness of assessing HRV at the onset, during, and immediately following the cessation of exercise to provide insight into autonomic stress reactivity. They, Michael et al. (2017), suggest that upon the onset of exercise, inputs from higherorder brain centers feed into the medullary control of cardiovascular function to reset the arterial baroreflex. This shift in the operating point for the arterial baroreflex requires a physiological shift in cardiac output that is met primarily by a withdrawal in parasympathetic input to increase HR. HRV provides context to changes in cardiac regulation beyond the assessment of mean HR values but these methods are limited to stationary time-series. The statistical estimation of a time-series (whether that statistic is assessing the variability or the complexity of the time-series) is dependent upon the underlying process producing the time-series and because each element of a non-stationary time-series can have a different entropy, these statistical calculations become consistently biased (Berry et al., 2020). Prior research has examined the relations between changes in HRV and the nonstationary (data with significant slopes in the HR response) Abbreviations: HR, Heart Rate; HRR, Heart Rate Recovery; HRV, Heart Rate Variability; RR rec, RR Intervals During Recovery; SDNN, Standard Deviation of RR-intervals; rMSSD, Root Mean Square of Successive RR Differences; ApEn, Approximate Entropy; SampEn, Sample Entropy; RMS, Root Mean Square Index; VO 2peak, Peak Oxygen Uptake; diff, Difference Calculation; res, Residual Calculation; MSE, Mean Square Error; ARS, Adaptive Regression Splines. trends in exercise recovery with conflicting HRV results (Arai et al., 1989;Bernardi et al., 1990;Perini et al., 1990;Breuer et al., 1993;Nakamura et al., 1993;Oida et al., 1997;Perini and Veicsteinas, 2003).
While a few studies have investigated methodological approaches for HRV indices during and/or after exercise, these studies have produced conflicting results (Arai et al., 1989;Bernardi et al., 1990;Perini et al., 1990;Breuer et al., 1993;Nakamura et al., 1993;Oida et al., 1997;Perini and Veicsteinas, 2003). As it pertains to exercise recovery, these methods are less commonly reported in the literature than traditional measures of HR and HRR. A recent review by Peçanha et al. (2017) examined the strengths and weaknesses of a variety of methods that have previously been utilized to assess post-exercise autonomic recovery via HRR and HRV indices. As discussed in this review (Peçanha et al., 2017), Goldberger et al. (2006) first introduced the idea of separating-out the 5 min period following the cessation of exercise into twenty 15s, ten 30s, and five 60s segments, which could be assessed by different variability indices to evaluate changes in cardiac autonomic reactivation. The piecewise analysis of these first 5 min following the cessation of exercise utilized HRV indices such as the standard deviation of normal RR-intervals (SDNN), the root mean square of successive RR differences (rMSSD), and the root mean square (RMS) (Ganio et al., 2009) index to assess cardiac autonomic reactivation. The authors concluded that both rMSSD and RMS of 30 s segments adequately reflected the autonomic changes occurring during the 5 min period of recovery (compared to 15 and 60 s segment lengths) and eventually expanded on these findings by evaluating these methods following a submaximal exercise bout (Ng et al., 2009). These studies show that indices of HRR and HRV immediately following the cessation of exercise reflect changes in cardiac autonomic modulation-as was also discussed in Michael et al. (2017).
However, measures of complexity-approximate entropy (ApEn) and sample entropy (SampEn), were not assessed in previous studies, but reflect important changes in heart rate dynamics that are not detected by traditional linear or frequency based measures of HRV (Goldberger Ary et al., 2000). Briefly, complexity measures stem from evaluation of nonlinear dynamics, which assess the structure of a time series and reflects the system's ability to adapt to internal or external perturbations. For instance, entropy indices quantify the degree of regularity or the appearance of repetitive patters within a time segment vs. traditional linear measures that examine the magnitude of variability. Thus, a low amount of entropy corresponds to a high amount of regularity within a signal (i.e., low adaptability of the system) vs. a higher degree of complexity that reflects a less-ordered signal (i.e., increased adaptability of the system) (West and Goldberger, 1987;Manor et al., 2010). A recent study by Berry et al. (2020), expanded on previous research by examining the influence of non-stationarity on indices of variability and complexity as well as various methods of detrending these data. These findings clearly show the biasing effects of non-stationarity on measures of variability and complexity as well as the effects of various methods of detrending on these indices.
To further elucidate the changes in heart rate dynamics immediately post-exercise, this study assesses the complexity during the entire 5 min recovery period immediately postexercise using both a previously reported method of assessing HRV during the recovery from exercise (Goldberger et al., 2006) and the methods outlined in Berry et al. (2020). In addition, we examine the relations between various physiological measures with indices of HRV from each of the respective methods.

MATERIALS AND METHODS
This analysis was performed on data that had been collected as part of the larger longitudinal RIGHT Track (RT) Health study. The overarching aim of RT Health is to examine the influence of early childhood self-regulation on the development of cardiovascular risk in adolescence and young adults. This study was approved by the Institutional Review Board at the University of North Carolina at Greensboro (IRB #11-0360). Prior to participating, all subjects provided written and informed consent.
The methods outlined below are specific to the HRV assessments and analyses; methods on the entire study have been published elsewhere (Wideman et al., 2016). Participants arrived at the laboratory between the hours of 0900 and 1400 at least 2 h post-prandial and having abstained from moderate to vigorous physical activity for 24 h, caffeine for 12 h, alcohol for 48 h. Body composition was assessed via air displacement plethysmography using COSMED's BODPOD R system. An incremental test to exhaustion was performed on a GE 2100 treadmill to assess peak oxygen uptake (VO 2peak ) (Parvo Medics TrueOne 2400). Preexercise, during exercise, and post-exercise RR intervals were recorded using a Polar V800 watch and downloaded following the completion of the visit. During the exercise test, participants self-selected their running speed while grade was increased by 3% every 2 min. Participants were asked to choose a speed that they felt comfortable maintaining for approximately 15 min so that the incremental increases in grade would result in maximal volitional fatigue within 8-12 min.
All mathematical and statistical procedures were performed using R 3.5.0. Artifact correction using linear interpolation of preceding-and proceeding-heart periods was completed prior to any additional analyses being performed using the open source "RHRV" package (Rodríguez-Liñares et al., 2008). Analysis of variance was performed using the "car" package (Fox and Weisberg, 2019) while multivariate adaptive regression splines were performed using the "earth" package (Milborrow, 2017). A example of an entire HR (bpm) and RR-interval (msec) profile is provided in Figures 1A,B, respectfully. Goldberger et al. (2006) The first 5 min of RR-interval data following the cessation of exercise were identified (Figures 1B,C). These segments were broken apart into twenty 15s, ten 30s, and five 60s segments that spanned the entire 5-min of recovery (Figure 2). Following linear regression analysis on the RR-intervals for each of the three different segment lengths, the residuals were computed and corresponding SDNN, rMSSD, and RMS indices were calculated.

Differencing and Polynomial Detrending of RR rec
Two separate detrending approaches were performed and comparisons between these two methods were evaluated. These methods of detrending non-stationary time-series are provided in Berry et al. (2020). Difference( diff ) Calculations. Differencing ( Figure 1E) was used to detrend the time-series and the resultant stationary time-series was subsequently analyzed to compare differences in variability and complexity. Residual( res ) Calculations. The optimal model was determined by purposefully overfitting the initial 5 min of RR rec data for all subjects. After modeling the initial 5 min of RR rec from 1 to 10 orders for each individual, a third-order polynomial was determined to be the optimal regression line ( Figure 1C), since it was flexible enough to account for individualities in this response and reduced the autocorrelation among the residuals without significantly overfitting the data. The residuals from these (individualized) lines were plotted as a secondary time-series ( Figure 1D) for analysis.

Variability and Complexity of the Detrended RR rec Data
The aforementioned transformations of the data were completed with the purpose of creating a time-series with a stationary mean to analyze the variability and complexity surrounding RR rec ; allowing us the ability to characterize the patterns of variability around the non-stationary decline in HR (rise in RR-intervals). Variability Measures. The equivalent of common time-domain measures of HRV, including SDNN of normal RR intervals and rMSSD, were used to assess the variability of the detrended-RR rec . Complexity Measures. ApEn and SampEn were used to assess the complexity of the residuals surrounding the non-linear trend in RR rec ; a higher number indicates a greater degree of complexity, whereas smaller values characterize more regular signals. While SampEn is more robust against changes in data-length compared to ApEn (Yentes et al., 2013), both indices were included in this analysis for continuity with previously published data.

Statistical Methods
With respect to our replication of Goldberger et al. (2006), repeated measures analysis of variance (ANCOVA) were used to assess changes in SDNN, rMSSD, and RMS across time between sex after controlling for race, body fat (BF), VO 2peak , and HR rest . Separate univariate ANCOVAs were used to compare differences in SDNN diff , rMSSD diff , SampEn diff , SDNN res , rMSSD res , and SampEn res between sex and race after controlling for changes in Frontiers in Physiology | www.frontiersin.org BF, VO 2peak , and HR rest . After these statistical comparisons were performed, exploratory statistical learning-methods were used to better understand some of the relations among these indices with various physiologic/demographic data.
Multivariate adaptive regression splines (ARS) produce continuous derivatives and are not only more powerful but also more flexible in determining the interaction among variables (Friedman, 1991). The original dataset was split, via random sampling, into a training-dataset (n = 46) and a testing-dataset (n = 46). Subsequently, the selected model was tested on the remaining data (the testing-dataset) to determine the validity of the models produced with the training-dataset. Physiological measures considered for each model included: height (Ht), weight (Wt), body fat (BF), body mass index (BMI), maximal oxygen uptake (VO 2peak ), resting HR (HR rest ), maximal HR (HR max ), HR at 1 min (HR one ), and 5 min (HR five ) following the cessation of exercise, recovery ( bpm) at 1 min (HRR one ) and at 5 min (HRR five ). These measures were scaled and centered prior to being used within the ARS analysis.

RESULTS
Data from N = 92 young adults (male = 46, female = 46) were included in this study. Subject demographics, separated by sex, are provided in Table 1.

Replication of Previously Established Methods of Assessment of RR rec
Results from replication of the methods established in Goldberger et al. (2006) are provided in Table 2 and visual representation of means from the segmented-analysis throughout the first 5 min of recovery are provided in Figure 2. Values are presented Mean (± SD). BF, Body fat; BMI, body mass index; VO 2max, Maximal oxygen update; HR rest, resting heart rate; HR max, maximal HR; HR one, HR at 1 min; HR five, HR at 5 min; HRR one, HR recovery at 1 min; HRR five, HR recovery at 5 min. Replication of these previously established methods in our young adults provided similar results to those of the healthy middle-aged adults reported in Goldberger et al. (2006). Extending on the comparisons previously made between healthy individuals and those with coronary artery disease (Goldberger et al., 2006), we examined the relationships between sex and race after controlling for BF, VO 2peak , and HR rest . While neither SDNN, rMSSD, or RMS were different between sex or race with any of the three segment-lengths, we did observe significant relationships between VO 2peak and HR rest with each of these analyses. More interestingly were the significant interactions between VO 2peak and HR rest with the changes in SDNN, rMSSD, and RMS across time and the lack of consistency in these findings across different segment-lengths ( Table 2).

Differencing and Polynomial Detrending of RR rec
Pearson correlations among demographic and the detrending (differencing-and residual-approaches) measures for the entire 5 min time series immediately post-exercise are presented in Table 3. The correlations between demographic characteristics and SDNN, rMSSD, and RMS of differenced and residual methods produced correlations ranging from −0.67 to 0.65. Similarly, the correlations between traditional indices of HRR with SDNN, rMSSD, and RMS ranged from −0.83 to 0.75.
Means from the newly proposed detrending methods are provided in Table 4 and results from the ANCOVAs are provided in Table 5. Similar to our findings with the segmenting-method, we did not observe significant group differences based on sex or race but we did observe significant associations between VO 2peak and HR rest with SDNN, rMSSD, ApEn, and SampEn for all but SampEn diff .
Results from ARS models are provided in Table 6. The residual-method provided a significantly (p = 0.04) lower Standard deviation of the normal RR-interval (SDNN); root mean square of successive differences (rMSSD); approximate entropy (ApEn); sample entropy (SampEn).

DISCUSSION
The major findings from these analyses were (1) reproducibility of previously established methods (Goldberger et al., 2006) for immediate post-exercise HRV in a young-adult population, (2) a unique dichotomy in the simple correlations between traditional measures of HRR and the complexity measures from the differencing-and residual-methods of RR rec , and (3) reliable models that demonstrate the complex nature of these indices. Additionally, these results provide preliminary insight into what physiologic variables may contribute to differences in these complexity indices when calculated on the detrended time series from the 5 min immediately post-exercise compared to traditional indices of HRR. The replication of the segmenting-method, as described in Goldberger et al. (2006), (used to assess acute changes in RR rec ) in our young adult population mirrored the patterns of change in SDNN, rMSSD, and RMS previously reported by Goldberger et al. (2006), however, the true raw values assessed for SDNN, rMSSD and RMS across the 5 min segment appeared to be different.  This discrepancy may be associated with the difference in age between the two samples as our sample included adolescents and young-adults compared to Goldberger et al. (2006) which included adults and older-adults. The individuals included in our sample have completed puberty and linear growth by this age, and these data fall within the normal range of values for young healthy adults. However, these individuals are still developing psychologically and physiologically, which may contribute to these differences. The ANCOVAs testing for effects with the segmentingmethod, as well as the new residual-method (from the entire 5 min period immediately following the cessation of exercise) proposed in this paper, produced similar results: no effects for sex or race were observed for any of the tests, while both VO 2peak and HR rest were significantly associated with SDNN, rMSSD, and RMS from both methods.
The correlations among various demographic variables and the variability and complexity indices from the differencingand residual-methods in Table 3 show a unique relation between VO 2peak and HR rest with HRR one , SampEn diff , and SampEn res . While the variability measures (SDNN diff , SDNN res , rMSSD diff , and rMSSD res ) were each significantly correlated with traditional metrics of HRR (VO 2peak , HR rest , HR max , HR one , HR five , HRR one , and HRR five ), SampEn diff and SampEn res had either no correlation with, or much-lower correlations with these indices. These findings suggest that indices of variability calculated on the entire 5 min segment using either differencing or residual methods, share large amounts of variance with traditional indices of HRR-in other words, the information provided by assessing either SDNN diff/res or rMSSD diff/res does not provide innovative or novel information beyond the traditional indices of HRR investigated in this, and many other, studies. Considering the physiological mechanisms associated with sympathetic withdrawal and the role of vagal input following the cessation of exercise to reduce HR, these observations are in-line with what is expected.
However, the complexity metrics obtained surrounding the increase in RR rec (i.e., SampEn diff , and SampEn res ) immediately following the cessation of exercise are differentially correlated with demographic measures and various traditional indices of HRR. This suggests that the complexity indices may provide distinct and unique information about the changes in HR immediately following exercise. While the ANCOVAs further supported these observations, the best support for this idea came from the ARS models (Table 6). Contrary to the ANCOVA analyses, the ARS models are exploratory and do not test specific hypotheses. These models progressively add predictors (and any combination of predictors), to the model and test for various non-linearities among them. The optimal model is determined by comparing changes in model fit. The optimal model was calculated on a training-dataset (a randomly selected 50% of the entire dataset) and then tested on the test-dataset.
Both SampEn diff and SampEn res were significantly associated with VO 2peak and neither was significantly correlated with HR rest , HRR one , or HRR five , but these variables (or some combination of these variables) were repeatedly selected in the ARS models as significant predictors of the complexity metrics (ApEn diff , ApEn res , SampEn diff , and SampEn res ). Similarly, HR rest , HRR one , and HRR five were commonly identified by ARS as significant terms in other models, suggesting that the variability surrounding the trend in RR rec immediately following the cessation of exercise does not provide information about changes in cardiac autonomic regulation that is not already observed in traditional indices of HRR. However, the repeated selection of HR rest , HRR one , and HRR five by ARS models to predict SampEn diff , or SampEn res despite the lack of significant correlations between any of these variables is an interesting finding that highlights the novelty of these detrending methods. These findings suggest that the complexity surrounding the trend in RR rec provides information about cardiac autonomic regulation not detected through traditional HRR indices or the variability surrounding the trend in RR rec .
Although the segmenting-method can be used to separate out the fast-phase and the slow-phase of the cardiac autonomic response to the cessation of exercise, we suggest that assessing the entire 5 min together may provide a more holistic and comprehensive assessment of the cardiac dynamics associated with RR rec . Furthermore, calculating the residuals from a third-order polynomial fit to the first 5 min of RR rec is less computationally intense and likely to be more robust against acute changes in RR rec compared to the segment approach previously described (Goldberger et al., 2006). When considering data length, it's important to remember that although time is a constant, there could potentially be drastic differences in the number of RR-intervals that occur within each segment being analyzed. Furthermore, the influence of irregular and/or erroneous beats would have a larger impact on the variability (McNames et al., 2003;McNames and Aboy, 2006;Thuraisingham, 2006) and complexity (Rhea et al., 2011) score calculated across a shorter segment when compared to a longer, 5 min, sampling period. Thus, we recommend the residual-method (as outlined in our methods) as a method of assessing the cardiac dynamics associated with recovery from exercise. In line with the findings of Berry et al. (2020) this method also produced more accurate and reliable models compared to the differencingmethod-observed through a significantly higher MSE in the ARS models calculated on the differenced data compared to the residual data.
Although the time course for parasympathetic reactivation post-exercise has been investigated previously, the exact timing of these changes is likely to be highly impacted by an individual's health and training status and thus, to vary from person to person. As such, the interindividual variability in postexercise parasympathetic reactivation may provide a nuanced assessment of cardiac autonomic regulation and subtle changes in this exercise-induced responsiveness may reflect the early stages of disease-an important target for future studies. While the development of an algorithm capable of reliably parsingout these points in the recovery phase is feasible, automating this analysis would be computationally intensive-more-so than the method (SampEn res ) outlined here, which we have shown can be predicted by non-linear relations among demographic information and traditional measures of HRR.
Several other methodological approaches have been suggested in the scientific literature to address HRV assessment immediately post-exercise, but a variety of concerns and confounds (Task-Force, 1996;Rhea et al., 2011;Peçanha et al., 2017), such as segment length, significantly limit the utility of these metrics in the exercise literature. Importantly, while many of these metrics are mathematically rigorous and computationally more feasible than ever before, the crux of their utility in exercise physiology hinges on their ability to translate into meaningful physiologic contexts. Herein lies one of the major advantages to complexity assessments surrounding the RR rec , since it translates easily into a meaningful physiological construct. Highly adaptive systems show increasing entropy (greater complexity) and are known to be associated with healthy responses (Shaffer and Ginsberg, 2017) and increased adaptability of the system (West and Goldberger, 1987;Manor et al., 2010).
In summary, we assessed the variability and complexity surrounding the non-stationary trend in RR rec through two methods including differencing and residual approaches. Results of SDNN, rMSSD, and SampEn calculations from the two detrending methods were similar, however, the multivariate ARS models produced from the residual approach had a lower MSE when compared to the differencing approach. This suggests that information about the individual recovery response immediately following the cessation of maximal exercise may be better represented through (third-order) polynomial regression detrending methods compared to differencing approaches. The complexity surrounding the trends in RR rec immediately following the cessation of exercise provides unique information and novel context related to cardiac autonomic regulation and the dynamics of RR rec not observed through traditional measures of HRR. Whereas the variability surrounding RR rec does not provide additional information beyond traditional HRR metrics already utilized in the literature. Further research is needed to establish the utility of this approach in other settings (i.e., field) and in other populations (i.e., older adults, clinical populations and athletes) as well as investigating other psychophysiological factors that may contribute to the dynamic regulation of cardiac control following the cessation of exercise. In addition, the relations of other commonly utilized metrics (e.g., spectral analyses) within the HRV literature should be applied to these methods.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Institutional Review Board at University of North Carolina at Greensboro. The patients/ participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
NB: data collection, data analysis, statistical analysis, and writing. EB: writing and editing. LS, SK, and SC: project design and editing. LW: project design, data collection, and editing. All authors contributed to the article and approved the submitted version.

FUNDING
This project was funded by the NICHD R01HD078346.