Assessment of Heart Rate Variability during an Endurance Mountain Trail Race by Multi-Scale Entropy Analysis

Abstract

The aim of the study was to analyze heart rate variability (HRV) response to high-intensity exercise during a 35-km mountain trail race and to ascertain whether fitness level could influence autonomic nervous system (ANS) modulation. Time-domain, frequency-domain, and multi-scale entropy (MSE) indexes were calculated for eleven mountain-trail runners who completed the race. Many changes were observed, mostly related to exercise load and fatigue. These changes were characterized by increased mean values and standard deviations of the normal-to-normal intervals associated with sympathetic activity, and by decreased differences between successive intervals related to parasympathetic activity. Normalized low frequency (LF) power suggested that ANS modulation varied greatly during the race and between individuals. Normalized high frequency (HF) power, associated with parasympathetic activity, varied considerably over the race, and tended to decrease at the final stages, whereas changes in the LF/HF ratio corresponded to intervals with varying exercise load. MSE indexes, related to system complexity, indicated the existence of many interactions between the heart and its neurological control mechanism. The time-domain, frequency-domain, and MSE indexes were also able to discriminate faster from slower runners, mainly in the more difficult and in the final stages of the race. These findings suggest the use of HRV analysis to study cardiac function mechanisms in endurance sports.

1. Introduction

The analysis of heart rate variability (HRV) has emerged as a powerful tool for assessing the status of the cardiovascular autonomic function in different diseases, both cardiological and non-cardiological [1,2,3,4,5]. Time-domain and frequency-domain HRV analyses have often been proposed for the quantitative assessment of cardiovascular control mechanisms [6]. Primarily because of the non-linear interplay of different physiological control loops in generating the heart rate, linear methods are not adequate to fully describe such a complex system. Several non-linear HRV methods such as fractal scaling analysis, higher order spectra analysis, multi-scale entropy analysis, power law analysis, complexity analysis, symbolic dynamics analysis and heart rate turbulence analysis have been studied for various diseases [7,8,9,10,11,12,13,14,15]. It must be taken into account that while the collection of heart rate (HR) was initially only possible with expensive laboratory-based electrocardiograph recorders, the recent availability of specifically designed portable recorders has substantially boosted the use of HRV monitoring [16,17].

HR varies according to individuals based on several factors, heredity, fitness level, exercise mode (e.g., endurance or resistance training), skill (economy of exercise), environmental (temperature [18], humidity, altitude [19]), state of mood [20] and diet, among others. Nevertheless, while HRV analysis is mainly associated with the prediction of sudden cardiac death [21,22] and with assessing cardiovascular [23] and metabolic illness [24] progression, recent observations have suggested it is also applicable to exercise and training [25]. During physical exercise, interactions occur between the cardiovascular system, the respiratory system, and the muscles. The autonomic nervous system (ANS) integrates their functions and plays a role in adapting to exercise effectively. HR is dependent on, and adjusted by ANS.

Although there are notable physical and physiological differences between athletes training for different sport activities, HRV is becoming one of the most used training and recovery monitoring tools in Sport Sciences [25,26]. Despite such variety, it is possible to apply HRV because cardiovascular autonomic regulation is an important determinant of training adaptations [27]. According to some authors, sustained training alters HRV caused by a modification of vagal activity with increased parasympathetic modulations [25,28]. Another study demonstrates that autonomic balance would be one of the first to be affected in case of fatigue-induced homeostatic disturbance [29]. Plews et al. [26] argued that HRV opens the door to effective monitoring of training adaptation. In accordance with the physiological adaptation theory, some experts affirm that carrying out higher intensity training sessions on days when parasympathetic activity has rebounded back to normal level or higher has a potential to increase training adaptation [30]. These studies suggest a training plan of decreasing the training stimulus when morning averaged HRV is lower and of training harder on days when HRV is normal or above normal.

The study of HRV is a potentially powerful basic scientific tool for better understanding cardiovascular system regulation. Exercise produces a physiological perturbation that significantly influences ANS activity. According to previous reports that examined HRV during exercise, HR increases immediately after starting an exercise because of withdrawal of the parasympathetic nervous activity and an increase in sympathetic nervous activity [31]. As the exercise load increases, HRV decreases [32]. Such behaviour has also been observed in ultra-endurance mountain trail races [33].

A comparison between recreationally active subjects and athletes of differing sports demonstrates that athletes exhibit an overall increase in HRV and parasympathetic cardiac modulation, while evidence suggests that high-intensity training can chronically lead to a shift from vagal to sympathetic cardiac modulation [27]. Recent laboratory studies confirm that HRV changes during exercise and recovery are affected by both the intensity and the physiological impact of exercise [34].

All previous studies except one [33], have analysed the effect of high-intensity training on HRV measured at rest, or on HRV measured during exercise in controlled laboratory conditions. However, it is unclear whether the adjustment of cardiac ANS in the laboratory reflects actual exercise in the field, particularly in adverse environments such as at moderate altitude. Also, during an endurance mountain trail race, runners are subjected to very diverse and continuously changing physiological demands. Consequently, the main purpose of the present study was to analyse the HRV responses to high-intensity exercise during a 35-km mountain trail race, as well as their relation to exercise intensity, fatigue level, terrain slope and running speed during the race. In addition, HRV response was compared based on race performance to ascertain whether fitness level could influence ANS cardiac modulation. As cardiovascular control is carried out by several regulatory mechanisms interacting across multiple temporal scales, the present work addressed the analysis of inter-beat interval (RR) time series by means of the evaluation of entropy rate over different time scales using multi-scale entropy (MSE). Traditional time-domain and frequency-domain analysis were also applied.

2. Materials and Methods

2.1. Data Analysed

The data were recorded as part of the SUMMIT Project, whose main goal is analysing whether the population engaged in high intensity physical activities for an extended period have a higher health risk that those who practice moderate or low physical activity. The participants were 15 male mountain trail runners that participated in the “Gran Volta Cerdanya Ultrafons”, a 35-km mountain trail race held in the Cerdanya region of the Spanish Pyrenees on 9 June 2013.The elevation profile of the race was: starting elevation of 1143 m above sea level (a.s.l.), minimum elevation of 1110 m a.s.l., peak elevation of 2066 m a.s.l., finish elevation of 1218 m a.s.l., peak elevation gain of 956 m, cumulative elevation gain of 1985 m, cumulative descent of 1910 m, and slopes between 20% and 30%. The local ethics committee approved the study, and all participants provided written informed consent.

After excluding data from four runners due to technical problems with the data collection, 11 runners were finally included in the study (age 38.9 ± 2.6 years, weight 76.0 ± 6.8 kg, height 178.2 ± 5.9 cm). All runners trained more than 10 h/week. Their race times ranged from 3:20 to 4:42 hh:mm and the average speed for the whole group of runners was 10.8 km/h. To compare different levels of performance, the participants were classified after the race into faster and slower runners (those that completed the race in <3.5 h and in >4 h, respectively). During the race all runners wore beat-to-beat HR monitoring devices (Polar RS800CX, Polar Electro Oy, Kempele, Finland) equipped with GPS to record distance and elevation data and RR-interval time series (intervals between consecutive heart beats). HR monitors of the same series have been validated [17,31].

2.2. Pre-Processing of the RR-Interval Time Series

The RR-interval time series (RR tachograms) were obtained from the HR monitors with a resolution of 1 ms. Filtering was required to remove anomalous RR interval values due to ectopic beats or artefacts related to detection errors, usually due to unusual movements or sweat during exercise. In this case, an adaptive filter was applied to RR series to detect each RR interval value differing by more than 6% of the five previous RR averaged intervals and to replace the RR interval by a normal-to-normal (NN) interval.

HRV analysis can be performed from 24-h recordings (long-term HRV) or from 5-min segments recordings (short-term HRV) [6]. In this study, a short-term HRV analysis was carried out to characterize HRV during different exercise loads along the race. For this, 36 consecutive segments of 5-min duration were selected from the RR time series of the faster runner, and the GPS positions corresponding to the initial values of each 5-min RR-segment were stored. These GPS positions were the reference points for the segmentation of the RR time series of the remaining runners. In this way, all the 5-min RR-segments of the different runners were synchronized with the GPS position and, therefore, the beginning of each segment corresponds to the same race conditions (distance from the start, elevation, and terrain slope). Figure 1 shows the elevation profile of the route and the 36 segments characterizing different exercise levels during the race according to the distance travelled from the start, elevation, and terrain slope.

2.3. Characterization of Each 5-Min Segment of the Course

2.3.1. Exercise Load Characterization during the Race

In order to characterize each of the 5-min segments during the race (henceforth s1, s2, …, s36), several performance variables were calculated—terrain slope, exercise load or exercise intensity, fatigue interval, and running speed. The terrain slope $TS$ (%) was calculated using the final and the initial altitudes, as well as the horizontal distance travelled, available from the GPS data. The exercise load $EL$ (a.u.) was calculated based on the “training impulse” method TRIMP for monitoring exercise load [35] by means of different updated procedures [36]. The fatigue interval FI (s) was calculated as the time elapsed since the start of the race up to a given segment. The running speed $RS$ (km/h) was determined by considering the distance travelled, calculated from the GPS data, during any given 5-min interval. Figure 2 shows the changes in these variables along the RR segments, calculated from the data of the fastest runner.

2.3.2. Heart Rate Variability Analysis

Standardized Time-Domain and Frequency-Domain Analysis

Based on the standards of the Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology [6], the following HRV indices from the linear time-domain were calculated for each of the 5-min RR-segments:

• meanNN (ms): Mean value of normal-to-normal (NN) interval time series.
• sdNN (ms): Standard deviation of each NN interval.
• rMSSD (ms): Square root of the mean squared differences of successive NN intervals.
• NNxx (n.u.): Number of interval differences of successive NN intervals greater than xx = {10, 20, 30, 40, 50} ms.
• pNNxx (%): Proportion of the NNxx in the total number of NN intervals.

Standardized frequency-domain indices [6] were determined from the Welch power spectral density estimate of equidistant RR time series, obtained by linearly interpolating the NN interval series at a resampling rate of 2 Hz. The following standard frequency-domain indices for short-term recordings were extracted:

• TP (ms²): Total power of the density spectra in the range ≤0.4 Hz.
• VLF (ms²): Power in very-low frequency range ≤ 0.04 Hz.
• LF (ms²): Power in the low frequency band 0.04–0.15 Hz.
• LFn (%): Normalized LF power LFn = LF/(LF + HF) × 100.
• HF (ms²): Power in the high frequency band 0.15–0.4 Hz.
• HFn (%): Normalized HF power HFn = HF/(LF + HF) × 100.
• LF/HF: ratio between LF and HF.

Multi-Scale Entropy Analysis

The multi-scale entropy (MSE) approach [10] was applied to quantify the degree of irregularity over a range of time scales. The method involves the construction of coarse-grained NN series and the quantification of the degree of irregularity of each of these.

MSE analysis uses the following procedure, given the NN intervals {NN₁, …, NN_i, …, NN_N}, consecutive coarse-grained time series ${y_j}^{\left(\tau \right)}$, corresponding to a scale factor τ, are constructed. First, the original times series are divided into non-overlapping windows of length τ. Second, the data points inside each window are averaged according to Equation (1):

$${y_j}^{\left(\tau \right)}=\frac{1}{\tau } {\sum }_{i=\left(j-1\right)\tau +1}^{j\tau }N{N}_i,1\le j\le N/\tau$$

(1)

Finally, the sample entropy (SampEn) [37,38] was applied for each coarse-grained constructed time series as a function of the scale τ varying from 1 to 20. The two parameters involved in the calculation of SampEn, the tolerance for accepting matches of two patterns (r) and the dimensional phase space (m), were set to 0.2 and 2, respectively.

2.4. Statistical Analysis

The Wilcoxon signed-rank test was applied to show differences between the HRV indexes calculated from the 36 segments of the race associated with different exercise load. The Student unpaired t-test was used to detect group differences between the HRV indexes obtained for the faster and slower runners (those that completed the race in less than 3.5 h and in more than 4 h, respectively). Correlations were assessed by Spearman rank-order correlation coefficients between the HRV indexes obtained from time-domain analysis, frequency-domain analysis, MSE analysis and the race performance variables of the runner. For each test, p-values are indicated.

3. Results

The following subsections include boxplots of the HRV indexes of the 36 RR-segments considered to analyse the changes in HRV indexes over the course. For each boxplot, the central mark is the median and the edges of the box are 25th and 75th percentiles. Values beyond the end of the whiskers are considered outliers and marked with “o”. The race profile is also shown to allow quick access to the test profile.

3.1. HRV Analysis During the Race

3.1.1. Standardized Time-Domain and Frequency-Domain Analysis of HRV

Among the indices, the meanNN and sdNN showed the most change across the different RR-segments during the race. Figure 3 and Figure 4 present the results for these indexes. Figure 3a gives the boxplot of the meanNN values for all runners. Figure 4 presents a 36 × 36-segment cross-tabulated Wilcoxon test p-values chart (as-log₁₀(p-values)) for the meanNN index. Cells containing significant p-values are coloured dark red (p-value < 0.005), medium red (0.005 ≤ p-value < 0.01), or light red (0.01 ≤ p-value ≤ 0.05). The diagonal of the matrix, which is symmetric, is coloured with white cells. Figure 3b is the boxplot of the sdNN values for all runners.

The meanNN index, which is associated with the HR itself, tended to increase over the course of the race implying that the average HR varied from higher to lower values (Figure 3a). In fact, the meanNN ranged from 363 to 381 ms in the first fourteen segments (high-ascent race interval), reflecting greater cardiovascular load, whereas in the last segments the values ranged from 375 to 402 ms, mirroring the less demanding profile and accumulated fatigue. The different response is observed in the p-values obtained when comparing segments s3–s8 and s15–s36 (Figure 4).

Segments s1–s14 correspond to a high-ascending part of the course. Even during this high-ascending interval, Figure 4 shows significant differences among segments, specifically segments s3–s7 compared with segments s9–s14. Then, during the descending stretch including segments s15–s18, heart rate tends to decrease (meanNN values are higher than previous ones) as is observed in Figure 3a. In this way, observing the significant levels in Figure 4, statistically high differences were found between segments s1–s12 (increasing slope hill) and segments s15–s18 (decreasing slope hill).

In segments s19–s30 (Figure 3a), which correspond to an ascending and followed by a flat interval, the mean heart rate increased to a value that lies between the values of segments s15–s18 (negative slope) and the values of segments s1–s13 (positive slope). Finally, the meanNN values of the last course intervals (segments s30–s36) increased again at segment s30 with values even higher than the meanNN of segments s15–s18 which correspond to the highest negative slope, and decreased to values similar to segments s19–s26 which correspond to an ascending slope.

The plot of the sdNN index, which is directly related to total HRV, is shown in Figure 3b. The sdNN index decreased at the beginning of the race and abruptly increased at s9 (midpoint of the first ascent). Then it decreased and increased again up to s17 (near to the first descending stretch). Then it stabilized up to s26 (onset of the flat stretch). After that point, sdNN showed an oscillatory trend up to the end of the race. The oscillations observed from s1 to s9, from s10 to s17 and from s18 to s26 showed a strong variability in the runners’ HR over the course of the race. This may show that their heart and cardiovascular system were capable of responding to the varying demands imposed by the changing terrain with associated changes in exercise load.

Changes in the sdNN index were oscillatory, but covered a wide range of values denoting a large variance among individuals. However, sdNN tended to decrease from s1 to s4 and then increase up to s9 when it reached the highest values of the entire course at the mid-point of the first ascending interval. In those first high ascending stretches (s1, s2 and s10), sdNN was lower than at s4–s8 (p-value < 0.005). Then, sdNN values decreased and increased again until segment s17 (near to first descent). This response was repeated up to s26 (onset of flat stretch), but with lower values as compared to s17. Segments s11–s12 (ascent) showed lower values compared to s18 segment (descent). The mean values of sdNN were about 10 ms from s26 but clearly increased in the last two segments (s35 and s36). Moreover, all segments from s3 to s8 were differentiated from segment s36 (p-value < 0.005).

sdNN values during some segments of the first high ascending stretch (s1–s8) showed many inter-segmental differences, especially between s4 to s8. Its values from s1 and s2 differed from those following the first ascending stretch but not from the remaining segments (except s18, s23 and s36). This index also differed between s3–s8 and s9–s10 (first steep ascent), s16–s21 (descending-ascending interval), and s26–s36 (flat-descending interval). Segments s11–s12 (also steep ascent) showed different sdNN indexes than s1–s2, s9–s10 (ascent), s18–s20 (descending-ascending stretch), and s29–s36 (final descent). Finally, there were also inter-segmental differences during the last descending interval. We note that the standard deviation of individual sdNNs varied along the segments of the race, had a wide range of values within each segment, and showed large inter-individual variability.

The rMSSD index, associated with the parasympathetic activity of the autonomic nervous system, showed significant inter-segmental variability between s1–s6 (lower rMSSD values around 4 ms) and s16–s18 (higher rMSSD values around 6ms). There were also differences between rMSSD values from s1–s6 and s29–s36 (lower rMSSD values from 5 to 7 ms). The terrain slope of s1–s6 was higher and so was HR compared to s16–s18 and s29–s36 (Figure 3a), which were descending stretches. The changes of NN10 and pNN10 over the course of the race were similar to those of the rMSSD index, as were the differences between segments. Considering that parasympathetic activity leads to a decrease in HR, these responses suggest a role for parasympathetic activity in descending intervals that is likely to be related to lower exercise load.

The LF and LFn indexes are associated with the overall ANS activity, including both sympathetic and parasympathetic activities, and with vasomotor oscillations responsible for the relative distribution of blood to different regions of the body. The LF index in segments s1–s3 was significantly different compared to s4–s6, s12, s14 and s21–s23. Consequently, the ANS responses differed between the first and second ascents. Segment s4 showed differences with s27–s28, s5 with s27 and 29, and s6 with the last segments of the first ascending stretch (s9–s10 and s14) and the final interval from the flat (after ascent) stretch (s26–s29). The LFn index showed few differences between segments.

Parasympathetic activity is responsible for HF and HFn indexes. Sympathetic cardiac nerves increase cardiac activity, whereas vagus nerves (parasympathetic branch) slow cardiac activity. The lower the HF values are, the higher the parasympathetic activity. Parasympathetic activity tended to decrease at the final stages of the race (in parallel with higher meanNN values) as compared to the initial stages (lower meanNN values). In fact, the main differences in HF index were found between the first high-ascending stretch (when the runners’ HR was highest) and the segments of the final stretch (when HR was lowest). The LF/HF ratio assesses the global sympathovagal balance. Figure 3c shows clear changes in LF/HF ratio between different stretches. There were differences between s1 (initial ascent) and s8–s9, s13 and s19 (ascent) and s28 (flat). There were also differences between s31 (descent after flat interval), and s7 (ascent), s26 and s29 (flat after ascent). It is important to point out that these segments correspond to stretches with varying external load.

3.1.2. Non-Linear Time-Domain Analysis of HRV: Multi-Scale Entropy

MSE indexes are related to complexity or regularity of the systems. The higher the MSE value, the higher the complexity, which means that there are more interactions between the ANS and the cardiac function, and that ANS activity is less regular. Figure 5a presents the evolution of the MSE indexes (mean ± standard error) as function of the scale τ, calculated from all thirty-six RR-segments for all runners. In addition, the number of these MSE indexes that were statistically differentiated between the thirty-six RR-segments is shown in Figure 5b. MSE indexes decrease from scale τ = 1 (SampEn, mean ± standard deviation = 1.20 ± 0.505) to scale τ = 6 (SampEn = 0.587 ± 0.246) and then monotonically increase to higher scales as at the scale τ = 20 (SampEn = 0.936 ± 0.419). When all segments of the course were compared, scale τ = 14 presented the highest number of indexes (#98) with p-value < 0.05 (Figure 5b), followed by scale τ = 6 (#93) and scale τ = 1 (#83).

The changes in MSE indexes from all runners at the scales τ = {1, 6, 14} are shown in Figure 6 by boxplot representations. A repeated pattern was present along the trail segments for each of the scales, as it can be observed in Figure 6a–c. The statistical analysis detected significant inter-segmental differences during the race, also shown in Figure 6 for selected τ scales. The SampEn responses at scale τ = 1 (Figure 6a) showed differences when comparing segments s1–s2 and s4–s8 (decreased regularity of RR signals), and between s4–s8 and s9 (increased regularity). At the first steep descents (s15–s20) there was a decrease of complexity, particularly between s4 and s18–s19 (p-value = 0.01). The final segment of the plateau s29 was different from s4–s8, and s30 from s10. Finally, the last s36 shown has the lowest complexity value of the SampEn and differed from s4–s8. In Figure 6b, the SampEn responses at scale τ = 6 were different when comparing s1–s3 and s4–s6 (decreased regularity), and between s4–s8 and s9–s10 (increased regularity). When analysing the final stretch of the first steep descent and subsequent initial ascent, s18 was different from s4–s6, and s19–s20 from s4 (increased complexity). However, no differences were found in the final segment s36 at this scale τ = 6. At scale τ = 14 (Figure 6c), there were differences between s2–s3 and s4–s6 (decreased regularity), and between s1–s8 and s10 (increased regularity). Related to segments at the plateau, s24–s29 showed a decreased complexity compared with s4.

Segments s32–s34, placed at the last part of the course, showed decreased regularity when compared with s10 (first ascent). Finally, the SampEn value of the final s36 differed from s1, s4–s8, s12, s15, s23 and s32, which presented greater complexity.

3.2. HRV Comparison Based on Different Race Performance

Table 1. Comparison of HRV index values between the faster and slower groups of runners. Significant differences are shown (p-value <0.05).

Index	Race Segment	Faster Runners	Slower Runners	p-Value (Unpaired t-Test)
TD analysis
sdNN (ms)	s9	24.6 ± 2.71	8.64 ± 1.89	0.001
	s25	6.19 ± 0.743	18.1 ± 6.48	0.034
	s36	7.48 ± 3.59	22.9 ± 5.58	0.016
rMSSD (ms)	s35	5.82 ± 0.321	3.80 ± 0.588	0.006
NN10 (n.u.)	s35	58.0 ± 22.9	5.00 ± 1.00	0.016
pNN10 (%)	s35	7.52 ± 3.09	0.651 ± 0.167	0.018
FD analysis
LFn (%)	s21	25.8 ± 11.1	61.3 ± 3.74	0.006
	s34	20.0 ± 11.7	61.5 ± 9.93	0.009
HFn (%)	s21	74.1 ± 11.1	38.6 ± 3.74	0.006
	s34	79.9 ± 11.7	38.4 ± 9.93	0.009
LF/HF	s21	0.370 ± 0.212	1.60 ± 0.263	0.003
	s34	0.270 ± 0.201	1.70 ± 0.621	0.019
MSE analysis
SampEn (τ = 1)	S9	0.277 ± 0.056	1.02 ± 0.172	0.002
SampEn (τ = 2)	S25	1.13 ± 0.220	0.212 ± 0.114	0.003
	S31	0.999 ± 0.091	0.479 ± 0.162	0.008
SampEn (τ = 3)	S9	0.142 ± 0.031	0.466 ± 0.097	0.005
	S31	0.795 ± 0.075	0.335 ± 0.052	0.001
SampEn (τ = 4)	S31	0.799 ± 0.120	0.297 ± 0.054	0.003
SampEn (τ = 9)	S31	0.963 ± 0.137	0.413 ± 0.078	0.004
SampEn (τ = 14)	S17	0.419 ± 0.123	0.881 ± 0.113	0.009

Data are mean ± standard deviation. TD, indexes from time-domain analysis; FD, indexes from frequency-domain analysis; MSE, indexes from multi-scale entropy analysis.

shows the HRV indexes that present statistical significant differences, in any of the 36 RR-segments analysed over the course, when comparing the faster and slower groups of runners. Table 1 presents, from left to right, HRV index, race segment (from s1 to s36), index values (mean ± standard deviation) of faster runners, index values (mean ± standard deviation) of slower runners and p-value (unpaired t-test). Regarding the time-domain HRV indexes, the sdNN values in segments s9, s25 and s36 differed between the two groups. The sdNN values were higher for the fast runners in s9 (steep ascent). However, at s25 (beginning of the plateau) and s36 (final race segment), sdNN was lower for the faster runners. Values of rMSSD, NN10 and pNN10 indexes were higher for the faster group at s35 (flat-descending stretch).

LFn, HFn and LF/HF indexes differed between groups at s21 and s34. In both segments, HFn was higher while LFn was lower in the faster runners. The lower values of LF/HF in the faster group indicate an increase of the spectral power in the HF band and an opposite response in the slowest group of runners. Furthermore, the rank of LFn, HFn and LF/HF at s21 and s34 were respectively similar.

The MSE analysis gave indexes with higher number of significant differences between the two groups of runners during the race. Statistical differences between both groups were observed in segments s9, s17, s25 and s31. Figure 7 shows the behaviour of the SampEn function of scale τ, where grey circles correspond to MSE values from the faster runners and black circles correspond to MSE values from the slower runners. The lines in black are averages of the SampEn at each τ and for each of the two groups. In segments s9 and s17, the SampEn values were lower for the faster runners. The opposite was observed in segments s25 and s31, where the SampEn values were lower for the slower runners. The largest group differences were observed in s9 (mid steepest ascent) at scales τ = 1 and τ = 3 (p-value ≤ 0.005), in s17 (end of steepest descent) at scale τ = 14 (p-value = 0.009), in s25 (start of the plateau) at τ = 2 (p-value ≤ 0.005), and in s31 (last descent) at scales τ = {2, 3, 4, 9} (p-value ≤ 0.008).

3.3. Course Characterization and HRV Analysis

Table 2. Statistically significant correlation coefficients (ρ) between HRV indexes and course variables (p-value < 0.0005).

	Exercise Load $\mathit{E}\mathit{L}$	Fatigue Interval FI	Terrain Slope $\mathit{T}\mathit{S}$	Running Speed $\mathit{R}\mathit{S}$
meanNN (ms)	−0.433	0.400	−0.300	0.228
sdNN (ms)	−0.315	0.206	-	-
rMSSD (ms)	−0.253	-	-	-
NN10 (n.u.)	−0.226	-	-	-
pNN10 (%)	−0.229	-	-	-

presents the statistical significant results from the correlation analysis between the HRV indexes and performance variables: EL (exercise load); FI (fatigue interval); TS (terrain slope); RS (running speed). There was a moderate positive relationship between the FI and the meanNN, as well as moderate negative relationships between EL and the HRV indexes meanNN and sdNN. A weak positive relationship was obtained between FI and sdNN, as well as between RS and meanNN. Weak negative relationships were obtained between EL and some HRV indexes: rMSSD, NN10 and pNN10, as well as between TS and meanNN.

4. Discussion

The present study analysed the HRV response to high-intensity exercise during a 35-km mountain trail race to ascertain whether fitness level could influence ANS cardiac modulation as characterized by time-domain, frequency-domain, and multi-scale entropy indexes. In the only study we are aware of in which HRV was investigated during a mountain trail race [33], neither frequency domain indexes nor non-linear indexes were considered. During an endurance mountain trail race the runners are subjected to very different conditions of exercise intensity, fatigue level, terrain slope and running speed. The incidence of these conditions to the HRV responses, which has not been previously reported, was also analysed.

There were nine main findings. First, the meanNN (inversely linked to HR) and sdNN (related to total HRV) increased over the course of the race in conjunction with intervals of greater exercise intensity and fatigue. Second, the rMSSD, NN10 and pNN10, associated with the parasympathetic activity, suggested a role for parasympathetic activity in descending intervals likely related to lower exercise load. Third, the LF and LFn, frequency indexes associated with both sympathetic and parasympathetic activity and blood flow distribution, suggested that ANS cardiac modulation varied greatly over the race and among individuals. Fourth, the HF and HFn, associated with parasympathetic activity, largely varied over the race and tended to decrease at the final stages. Fifth, most changes of the LF/HF ratio, index of global sympathovagal balance, corresponded to intervals with varying exercise load. Sixth, the MSE indexes, related to complexity or regularity of systems, indicated the existence a very large amount of interactions between the heart and its neurological control mechanism. Seventh, the meanNN, sdNN, rMSSD, NN10 and pNN10 (time-domain), LFn, HFn, and LF/HF (frequency-domain), together with MSE indexes at difference scales, were indicators that differentiated faster from slower runners, mainly at the harder and final stages of the race. Eight, the MSE indexes provided the largest amount of inter-segmental differences between the two groups of runners, mostly in the first steep-ascending stretch and in the last flat and steep-descending intervals. Finally, the only noteworthy correlations between race performance variables and HRV indexes were provided by time-domain indexes (meanNN, sdNN, rMSSD, NN10 and pNN10).

The meanNN index, mirroring HR, increased from the beginning to the end of the race following changes in elevation (i.e., lower during ascents and higher during descents). This suggests that meanNN, indicator of the sympathetic activity of the ANS, varies mainly as a function of the exercise intensity, fatigue level, terrain slope and speed, with which it was significantly correlated (Table 2). Considering the sign of the associations, while exercise intensity and terrain slope were responsible for the increase in HR (decreased meanNN), fatigue and running speed had the opposite effect (increased meanNN). In fact, the comparison between the first ascent (s1–s12) and the subsequent descent (s15–s18) showed that the highest meanNN values occurred in the negative slope segments. This was also observed when analysing the last descent of the race (s30–s36). The changes in sdNN, which is related to the total HRV, contained oscillations (decreasing values followed by increasing values and vice-versa) showing a varying HRV during the course. Its association with exercise intensity and fatigue (Table 2) suggests that this index is also an indicator of the sympathetic activity.

The time-domain indexes rMSSD, NN10 and pNN10, associated with the parasympathetic activity, showed large inter-segmental variability between ascent and descents stretches. They were inversely correlated with exercise intensity, suggesting a role for parasympathetic activity in descending intervals related to decreased exercise load. In the only study we are aware of in which HRV was investigated during a mountain trail race, a decrease of meanNN, sdNN, rMSSD, and pNN50 was reported, indicating an increase in sympathetic modulation and a decrease in parasympathetic activity [33]—this is in line with our results. Considering only resting HRV analysis, similar changes were obtained when comparing before versus after a 14-h mountain running event [39] and an Ironman triathlon event [40]. It has been suggested that these changes may be a factor involved in cardiac adaptations in ultra-endurance athletes [33].

LF and LFn frequency indexes are associated with the overall activity of the ANS, including sympathetic and parasympathetic activities, and with vasomotor oscillations that are responsible for the relative distribution of blood to different regions of the body. The LF index indicated that the ANS responses differed between the first and second ascents, while the LFn index showed few inter-segmental differences. These indexes suggested that ANS cardiac modulation varied greatly over the race and among individuals because of varying terrain conditions and training level. Also, HF and HFn, associated with parasympathetic activity, varied greatly over the race, and tended to decrease at the final stages, confirming their role as indicators of vagal modulation. Most changes observed in the LF/HF ratio, an index of global sympathovagal balance, corresponded to intervals with varying exercise load (i.e., changes of terrain slope), suggesting that the autonomic control finely tunes cardiac function according to the exercise load variations met during the race.

The MSE indexes are markers of irregularity and unpredictability (complexity) of the system under investigation [37]. The higher the MSE value, the higher the system complexity, which means that there are more interactions between the ANS and the cardiac function and less regular activity of the ANS. The present results suggest that there are many interactions between the heart and its neural control mechanism. The analysis of MSE indexes at the scales with the highest number of significant differences when comparing race segments (τ = 1, τ = 6 and τ = 14) showed the same oscillations (i.e., decreasing values followed by increasing values) as described for the sdNN index. The magnitude of estimated exercise intensity during the course may justify this behaviour as described in the sdNN analysis. However, the increase of MSE indexes from s1–s2 to s4–s8 denotes an increase of the complexity synchronized with a decrease in exercise intensity. On the other hand, the decrease of MSE indexes from s4–s8 to s9 indicates a decrease of the complexity synchronized with an increase in exercise load. The MSE analysis at different scales τ also provided indexes differentiating both groups of runners (Table 1). At the first ascent, SampEn values were lower for the faster runners indicating a lower ANS complexity. Conversely, the same group exhibited higher SampEn at s25 and s31 values (flat and final descent) denoting a higher ANS complexity. This response supports the concept that MSE indexes reflect the interplay between sympathetic and parasympathetic cardiac modulation during intense exercise.

This analysis of the HRV indexes comparing the faster and slower groups of runners also showed that, in the time-domain analysis, sdNN index in s9 (mid first ascent) was higher for the faster group of runners, showing higher total HRV. Conversely, sdNN in s25 (beginning of the plateau) and s36 (last segment) was lower for the faster group, showing a reduced total HRV. rMSSD, NN10 and pNN10, indicators of vagal modulation, were higher for the faster group at the end of the race. The frequency-domain analysis showed differences between both groups in s21 and s34. HFn, indicator of vagal modulation, was also higher in the faster group of runners. The lower value of LF/HF in the faster group indicated a lower global sympathovagal balance at the middle and at the end of the race. These results showed that better runners exhibited a lower HRV and greater parasympathetic activity when confronted with higher exercise loads during the race.

Concerning the effects of the low-moderate altitude in which the race took place on HRV (1143 to 2066 m a.s.l.), previous work has shown that acute effects of altitude exposure and exercise on HRV were found during exercise >2500 m a.s.l) and mainly for the harmonic components of HRV [41].

Globally, these complex cardiac responses during varying intense exercise suggested that the runners’ heart and cardiovascular systems can adapt to the external changes of the terrain, and the physiological loading imposed by the heavy exercise, through the neural control of the heart and cardiovascular system.

Considering the reduced number of runners that wore beat-to-beat HR monitoring devices during the race, these results need further validation with a higher number of subjects. The obtained results encourage the use of other nonlinear analysis methods as future work.

5. Conclusions

The autonomous neural cardiac modulation analysed by the multi-index characterization of HRV during a 35-km endurance mountain trail race shows multiple differential changes, mostly related to exercise intensity (linked to terrain and slope) and fatigue. The changes of the frequency domain and non-linear domain HRV parameters during an endurance trail race have not been previously described. The changes were characterized by increased meanNN (inversely linked to HR) and sdNN, associated to sympathetic activity, and by decreased rMSSD, NN10 and pNN10, related to parasympathetic activity. LF and LFn, associated with sympathetic and parasympathetic activity and blood flow distribution, suggested that ANS cardiac modulation varied greatly over the race and among individuals. HF and HFn, also associated with parasympathetic activity, varied greatly over the race, and tended to decrease at the final stages, whereas changes of the LF/HF ratio, index of global sympathovagal balance, corresponded to intervals with varying exercise load. MSE indexes, related to complexity of systems, indicated the existence of a very large amount of interactions between the heart and its neurological control mechanism. MeanNN, sdNN, rMSSD, NN10 and pNN10 (time-domain), LFn, HFn, and LF/HF (frequency-domain), together with MSE indexes at difference scales could discriminate between faster from slower runners, mainly at the harder and final stages of the race. The only noteworthy correlations between race performance variables and HRV indexes were provided by time-domain indexes (meanNN, sdNN, rMSSD, NN10 and pNN10). These findings open a door to the use of HRV analysis to study the mechanisms of cardiac function in endurance sports. Their main clinical implications reside in the ability of different HRV indexes to reflect the cardiovascular load during high-intensity exercise performed at low-moderate altitude-beyond the simple heart rate response, and to relate some of these indexes of autonomic control with the performance capacity of endurance athletes.