Normative Values for Heart Rate Variability Parameters in School-Aged Children: Simple Approach Considering Differences in Average Heart Rate

Background: Heart rate variability (HRV) analysis is a clinical tool frequently used to characterize cardiac autonomic status. The aim of this study was to establish normative values for short-term HRV parameters by considering their main determinants in school-aged children. Methods: Five-minute electrocardiograms were taken from 312 non-athlete children (153 boys) at age of 6 to 13 years for computation of conventional time- and frequency-domain HRV parameters. Heart rate (HR), respiratory rate, age, body mass index, and sex were considered as their potential determinants. Multiple regression analysis revealed that HR was the principal predictor of all standard HRV indices. To develop their universal normative limits, standard HRV parameters were corrected for prevailing HR. Results: The HRV correction for HR yielded the parameters which became independent on both sex and HR, and only poorly dependent on age (with small effect size). Normal ranges were calculated for both time- and frequency-domain indices (the latter computed with either fast Fourier transform and autoregressive method). To facilitate recalculation of standard HRV parameters into corrected ones, a calculator was created and attached as a Supplementary Material that can be downloaded and used for both research and clinical purposes. Conclusion: This study provides HRV normative values for school-aged children which have been developed independently of their major determinants. The calculator accessible in the Supplementary Material can considerably simplify determination if HRV parameters accommodate within normal limits.

Another important factor, usually not accounted for in HRV normative studies but significantly influencing both HR and HRV, is breathing (Quintana et al., 2016a;Shader et al., 2017). The respiratory rate declines from birth to early adolescence (Fleming et al., 2011), and consequently such changes may influence HRV (Billman, 2011;Quintana and Heathers, 2014;Quintana et al., 2016a,b;Laborde et al., 2017).
The aims of the present study were to determine the main determinants of HRV in healthy school-aged children, and by incorporating the significant determinants to define normative values for both time-and frequency-domain HRV indices.

Study Population
The whole study group consisted of 346 children of both sexes. The inclusion criteria to this study were as follows: (i) age between 6 and 13 years, (ii) absence of diseases and/or regular use of medications affecting the cardiopulmonary system and/or interfering with the autonomic nervous system, and (iii) not being an active athlete in any sports. The parents/legal guardians were interviewed about children's diseases and/or medications (the school health records concerning health status were additionally verified). Fifteen subjects were excluded from the analysis due to suspicion of cardiac/non-cardiac diseases. The group participated in an earlier study, its details and description of the procedures imposed before or during the proper electrocardiographic examinations have been published elsewhere (Gasior et al., 2015). Since there are data indicating an association of cardiac autonomic function with physical activity in children and adolescents (Oliveira et al., 2017), we excluded 19 active athletes (Araújo and Scharhag, 2016) from the present study providing normative values. Consequently, 312 healthy Caucasian children of both sexes between the ages of 6 and 13 years were enrolled in the final analysis. The group was divided into four subgroups according to age: (I) 6 ≤ age < 8 years; (II) 8 ≤ age < 10 years; (III) 10 ≤ age < 12 years; and (IV) 12 ≤ age < 14 years. This classification is in accordance with other studies where participants younger than 10 years were considered as children and those between 10 and 14 years as preadolescents or early adolescents (Frigerio et al., 2006;Cysarz et al., 2011;Varlinskaya et al., 2013). The study was approved by the University Bioethical Committee and followed the rules and principles of the Helsinki Declaration, as well as all parents or legal guardians gave their informed written consent.

Body Mass Status Measurement
The body mass status was measured using Body Mass Index (BMI) defined as body mass in kilograms divided by height in meters squared. BMI is widely used for assessing the weight status of pediatric participants since it can identify those who are overweight, at risk of being overweight, or underweight based on age and sex (Barlow and Dietz, 1998).

Electrocardiographic (ECG) Recordings and HRV Analysis
Twelve-lead, 6-min ECG recordings (sampling frequency = 1000 Hz) were performed using a portable PC with integrated software system (Custo cardio 100 12-channel PC ECG system; Custo med GmbH, Ottobrunn, Germany) during regular school days between 8 a.m. and 2 p.m. in a supine position. For heart rate stabilization, children were placed in a supine position for 5 min before the examination. All ECG recordings were preprocessed before HRV analysis, i.e., they were visually inspected for potential non-sinus or aberrant beats and such erroneous beats were properly corrected (interpolation of degree zero) (Peltola, 2012). Time-and frequency-domain HRV analyses were performed from 5-min ECG recordings by using Kubios HRV Standard 3.0.2 software (University of Eastern Finland, Kuopio, Finland) (Tarvainen et al., , 2017. Standard deviation of R-R intervals (SDNN), root mean square of successive R-R interval differences (RMSSD), and pNN50, which denotes percent of R-R intervals differing >50 ms from the preceding one, were determined. Before calculating spectral HRV parameters, smoothness priors based detrending approach was applied (smoothing parameter, Lambda value = 500) (Tarvainen et al., 2002), and then R-R interval series were transformed to evenly sampled time series with 4-Hz resampling rate. The detrended and interpolated R-R interval series were used to compute HRV spectra by employing fast-Fouriertransform (FFT) with Welch's periodogram method (300 s window width without overlap) and by using the autoregressive model (AR) with model order set at 16 (Boardman et al., 2002). The following spectral components were distinguished: very low frequency (VLF, 0-0.04 Hz), low frequency (LF, 0.04-0.15 Hz), high frequency (HF, 0.15-0.50 Hz), and total power [(in two versions, i.e., TP 1 , 0-0.5 Hz (VLF+LF+HF) and TP 2 , 0.04-0.50 Hz (LF+HF)] in absolute units (ms 2 ), as well as nLF, nHF in normalized units (nu) (Task Force, 1996).

ECG Derived Respiration
The respiratory rate (RespRate) was calculated from ECG recordings according to Sinnecker et al. (2014). The time series of the QRS amplitude for all 12 ECG leads were calculated and then local maxima (i.e., data points with values greater than both the preceding and the following data point) were established. Then, the mean interval between consecutive local maxima for each time series was computed and the reciprocal value of the mean maximum-to-maximum interval was obtained. The median of these reciprocal values over all-time series represented RespRate used in this analysis (Sinnecker et al., 2014).

HRV Correction
Since HRV is dependent on HR due to both physiological (i.e., the autonomic nervous system influence) and mathematical (i.e., the inverse non-linear relationships between HR and R-R interval) reasons, it is hard to ascertain if clinical significance of HRV comes from the variability or from HR (Sacha and Pluta, 2008) it is even more important in populations where average HR is changing, like in children during their grow and development. Therefore, to investigate the "objective" variability, one should get rid of the HR impact on HRV. To do so, one should correct the HRV for the prevailing HR. The correction procedure employed in this study relied on the division or multiplication of standard HRV indices by different powers of their corresponding mean R-R interval (mRR) -the higher HRV dependence on HR, the higher power of mRR must be used to get rid of this dependence. Specifically, if HRV parameters revealed negative relation with HR (i.e., SDNN, RMSSD, pNN50, VLF, LF, HF, TP 1 , TP 2 , and nHF), they should be divided by the suitable power of mRR in order to become HR independent, however, if the HRV parameters are positively related with HR (i.e., nLF and LF/HF), they must be multiplied by the appropriate power of mRR -this way we were able to remove both physiological and mathematical HRV dependence on HR (Sacha et al., 2013b,c).

Statistical Analysis
The Kolmogorov-Smirnov test was used to assess the normality of the data distribution. Normally distributed values were presented as mean ± standard deviation, but others as median and interquartile range (IQR). Accordingly, Pearson's correlation coefficient (r) or Spearman's rank correlation coefficient (R) were used to assess the relationship between variables, and Student's t-test or non-parametric Mann-Whitney test were employed to determine differences. The Pearson's chi-squared test was applied to determine differences between variables prevalence in multiple groups. Values of the standard and corrected HRV parameters were compared between four age subgroups by using Kruskal-Wallis (K-W) test and analysis of variance (ANOVA), respectively. Multiple regression analysis was carried out to identify independent determinants of standard and corrected HRV parameters -variables with a skewed distribution were logarithmically transformed (ln). To avoid multicollinearity, redundant variables were removed from the multivariate regression models in the case of pairwise correlations between continuous variables (Slinker and Glantz, 2008). The effect size was measured by calculating Cohen's f 2 within a multiple regression model (Selya et al., 2012). The f 2 were calculated for a combined prediction of the model, and Cohen's f 2 variation was computed to measure the local effect size (Cohen, 1992;Selya et al., 2012). According to Cohen's guidelines, f 2 ≥ 0.02, f 2 ≥ 0.15, and f 2 ≥ 0.35 represent small, medium, and large effect sizes, respectively (Cohen et al., 2003). The threshold probability of p < 0.05 was taken as the level of statistical significance. Since all HRV parameters were not normally distributed, normative ranges have been given as medians with the 5th and 95th percentiles. Statistical calculations were performed by using the software STATISTICA 10 (StatSoft, Inc., Software, Tulsa, OK, United States).
Heart rate, RespRate, age, BMI and sex were initially chosen as potential determinants of standard HRV parameters ( Table 1 presents their correlations with HRV). All these variables were associated with each other with one exception, i.e., there was no significant correlation between RespRate and BMI (Figure 1). Since the association between RespRate and HR was stronger than between RespRate and HRV, as well as, the association between BMI and age was stronger than between BMI and HRV, (Figure 1 and Table 1), both RespRate and BMI were considered to be redundant for the multivariate analysis (collinearity tests confirmed this redundancy). Accordingly to the statistical assumptions of the multiple regression analysis (Kraha et al., 2012;Ray-Mukherjee et al., 2014) and taking into account the practicality of normative data in a clinical context, the following variables were included in the regression analysis: HR, age, and sex. Table 2 exhibits results of this analysis for time-and frequency-domain HRV parameters calculated with FFT -the results for HRV spectra estimated with AR were similar and are presented in Supplementary Table S1. The models accounted for 9-60% of the entire HRV variance, but importantly, HR was the strongest predictor for all analyzed standard HRV parameters. HR was the most powerful determinant in the whole group ( Table 2 and Supplementary Table S1) and also within the age subgroups: 6-7, 8-9, 10-11, and 12-13 years (Supplementary Tables S2-S9).
Frontiers in Physiology | www.frontiersin.org FIGURE 1 | Association between variables selected as potential independent predictors of standard HRV parameters. BMI, body mass index; HR, heart rate; HRV, heart rate variability; RespRate, respiratory rate; R, Spearman's rank correlation coefficient; r, Pearson's correlation coefficient.
However, there was a significant overlap in the normative ranges between consecutive HR quartile subgroups which makes such data impractical. Therefore, considering the fact that HR was the most influential factor determining all standard HRV parameters, we excluded its impact by correcting the standard HRV parameters for their prevailing HR. This way we were able to get normative values for so called corrected HRV, i.e., the variability independent on average HR. Indeed, standard timedomain HRV parameters lost their dependence on HR after dividing SDNN, RMSSD, and pNN50 by mRR to the power: 2.2, 3.0, and 5.0, respectively. In spectral analysis, VLF obtained by FFT and AR lost its dependence on HR after dividing by mRR to the power 3.0 and 4.0, respectively. Other frequency-domain HRV indices lost their association with HR after dividing LF, HF, TP 1 (VLF+LF+HF), TP 2 (LF+HF), and nHF by mRR to the power: 4.0, 5.0, 5.0, 5.0, and 0.5, respectively (for both FFT and AR). The nLF and LF/HF (calculated with either FFT and AR) stopped being dependent on HR after multiplying by mRR to the  (Sacha et al., 2013b). Correlation coefficients between HR and all corrected HRV parameters were not statistically significant and ranged between −0.09 and 0.08 (p > 0.10 for all).
In the multiple regression analysis, age turned out to be the strongest determinant of the corrected time-and frequencydomain HRV parameters calculated with FFT (Table 4) -the results were the same when HRV spectra were estimated with AR (the AR data are presented in Supplementary Table S10). However, the calculated models accounted for only 2-9% of the whole variance of the corrected HRV indices. The age contribution to this model was quite low (i.e., β-value ranged between −0.13 and −0.29), and the effect size for the combined model, as well as for the variables taken individually (age and sex) was small (f 2 ≤ 0.091 for all). Therefore, practically there was no reason to establish normative values according to age, and we calculated normal limits for the overall study group ( Table 5). Of note, due to division by high powers of mRR, the numbers in Table 5 are long (i.e., they present many zeros after the decimal place), consequently calculations and interpretations of such results may be uninterpretable in clinical situations, hence in the Supplementary Material one may find the Excel sheet calculator which recalculates standard HRV parameters into corrected ones and determines whether a given value of corrected HRV parameter is within normal limits (Supplementary Table  S11, Sheet "Corrected HRV Calculator").
values for HRV parameters need to be established. This study provides such normative quantities for short-term (i.e., 5 min) time-and frequency-domain HRV parameters based on a group of 312 school-aged children, i.e., at age of 6-13 years. Since HRV is under influence of a number of factors, we conducted the multiple regression analysis which revealed that HR is the strongest predictor for all standard HRV parameters in the overall study group and in each analyzed age subgroups ( Table 2 and Supplementary Tables S1-S9). Age was the second independent predictor but accounted only for a small percentage of the entire HRV variance. Interestingly, in the multivariate analysis sex did not independently determine HRV in our study population ( Table 2 and Supplementary Table S1). The reason for this possibly comes from the observation that boys usually present significantly lower HR than girls and therefore sex differences in HRV may result from differences in HR. When considering HR, we could ignore differences between the sexes and establish normative HRV values only for subgroups categorized according to the quartiles of HR (Table 3). However, among such specified subgroups, confidence intervals significantly overlapped for most of HRV parameters, which may contest practicality of this approach. Hence, we calculated HRV independent on HR by correcting standard HRV indices for prevailing HR. In this method, we obtained so-called corrected HRV parameters, which reflect "objective" variability regardless whatever the average HR is. In the multiple regression analysis, age was significant predictor for such corrected HRV indices (Table 4 and Supplementary Table S10), nevertheless, its contribution to the whole variance of HRV was small, i.e., the models only accounted for up to 9% of the HRV variance and the age role in these models was modest (β values ranged between −0.29 and 0.13) as well as its effect size turned out to be small (f 2 ≤ 0.091) ( Table 4 and Supplementary  Table S10). Therefore, we omitted age impact in these 6 to 13-year-old children and established normal limits for the entire study group ( Table 5). The calculation of corrected HRV parameters and their assessment become simpler by using the file provided in the Supplementary Material, which recalculates standard HRV indices into corrected ones and automatically determines whether they accommodate within normal limits (Supplementary Table S11, Sheet "Corrected HRV Calculator").
The influence of HR on HRV in children has also been recognized in other studies (Goto et al., 1997;Massin and von Bernuth, 1997;Jarrin et al., 2015;Bjelakovic et al., 2017;Herzig et al., 2017). Indeed, Jarrin et al. (2015), accounted the observation that HR is the strongest factor determining HRV and presented their normative values adjusted for HR, however, only in participants with a very narrow age range, i.e., 9-11 years (mean ± SD: 10.2 ± 0.3 years) (Jarrin et al., 2015). Yet, the problem of the strong association between HR and HRV can simply be resolved by correcting HRV parameters for prevailing HR (Sacha and Grzeszczak, 2001;Pluta, 2005a,b, 2008;Sacha, 2013Sacha, , 2014aSacha et al., 2013aSacha et al., ,b,c, 2014Gasior et al., 2016) and this method turns out to be feasible and convenient also in children populations (Gasior et al., 2015).
More and more clinical and methodological papers highlight that respiration should be taken into account, or at least monitored, in HRV studies (Brown et al., 1993;Williams and Lopes, 2002;Eckberg, 2003;Song and Lehrer, 2003;Yildiz and Ider, 2006;Beda et al., 2014;Heathers, 2014;Lehrer and Gevirtz, 2014;Quintana and Heathers, 2014;Quintana et al., 2016a;Laborde et al., 2017). In population with known fast breathing rate, like children, respiratory depth, and frequency are particularly associated with heart rate and its fluctuations and may influence HRV data independently of cardiac autonomic activity (Eckberg, 2009;Quintana et al., 2016b). Also in our study, we observed significant correlation between children's respiratory rate and HR, nevertheless, the correlation between HRV and HR was stronger than the association between HRV and breathing rate, and the statistical analysis revealed that respiratory rate was redundant in the multiple regression analysis. In our recent study, we have shown that the influence of breathing rate on HRV appears to be, at least in part, HR dependent, i.e., after HRV correction for prevailing HR, the correlation between breathing rate and HRV decreases (Gasior et al., 2016).
In the present study, we calculated VLF power from shortterm ECG recordings (i.e., lasting 5 min) which may be problematic from methodological point of view (Task Force, 1996;Heathers, 2014). However, this parameter was also considered in the previous study concerning reference values for short-term (5 min recordings) standard HRV in children (Michels et al., 2013). Shaffer and Ginsberg in an overview paper concerning HRV metrics and norms state that "the VLF band requires a recording period of at least 5 min" (Shaffer and Ginsberg, 2017). Moreover, in our previous studies addressing large populations of patients after myocardial infarction, VLF parameter calculated from 512 R-R intervals turned out to be the strongest predictor of different causes of death (Sacha et al., 2013a(Sacha et al., , 2014. Therefore, we included this parameter in the present analysis despite the short-term nature of ECG recordings. This may enable future studies aiming to explore its clinical value in pediatric populations with various disorders and consequently may help to understand the nature of this parameter. To make our data pragmatic and comparable for other studies: (i) we calculated short-term time-and frequency-domain HRV parameters by using actual version of a very popular and free tool: Kubios HRV Standard ver. 3.0.2 (Tarvainen et al., , 2017; (ii) in spectral analysis, we used frequency bands up to 0,5 Hz which is appropriate for children populations; (iii) we presented frequency-domain HRV parameters calculated with both fast Fourier transform and autoregressive method; and (iv) we provide the Excel sheet calculator to facilitate computation and assessment of corrected HRV parameters.

Limitations
A limitation of this study is the relatively narrow age range for which these normal limits can be applied, i.e., 6-13 years. However, this age period corresponds to the considerable developmental changes in both circulatory and neurological systems, and therefore the HRV normative values for this life period may be very useful in different clinical scenarios.
Finally, the norms can only be applied in Caucasian children. Caution is needed in applying the normative HRV values for other populations as the significant ethnic differences of HRV in children and youth were shown (Wang et al., 2005;Reed et al., 2006;Eyre et al., 2013). Future studies using approach considering differences in average HR may help to identify ethnic differences in corrected HRV indices.

CONCLUSION
This study provides HRV normative values for school-aged children which have been developed independently of their major determinants, especially related to average HR. An Excel sheet calculator is provided to increase accessibility (Supplementary Table S11, Sheet "Corrected HRV Calculator") and simplifies determination of HRV parameters from an individual child and if values are within normal limits.

AUTHOR CONTRIBUTIONS
JG, JS, and MP conceived and designed the experiment. JG and MP acquired the data. AT, TK, and BW analyzed the data for ECG assessment. JZ analyzed the data for respiratory rate calculations. JG and JS analyzed the data for HRV and statistical analysis. JG, JS, MP, JZ, PJ, AT, TK, BW, and MD interpreted the data. JG, JS, MP, JZ, PJ, AT, TK, BW, and MD drafted the work and revised it critically for important intellectual content. JG, JS, MP, JZ, PJ, AT, TK, BW, and MD approved the final version of the manuscript to be published.

FUNDING
The authors and clinical organizations with which the authors are affiliated or associated did not receive any grants or outside funding in support of the research for or preparation of the manuscript, did not receive payments or other benefits from a commercial entity, do not have any professional relationships with companies or manufacturers who will benefit from the results of the present study. The aforementioned disclosure also applies to the authors' immediate families.

ACKNOWLEDGMENTS
We gratefully acknowledge Professor Craig Anthony Williams for his constructive comments and language correction. We thank the children for participation in the study and their parents/legal guardians for permission to conduct the measurements. We also thank the teachers, nurse, and school staff for help and cooperation during the study.