1 Introduction

Aerobic metabolism, measured universally via oxygen uptake () [1], is the principal mechanism by which humans generate energy from ingested foodstuffs for life. Physical activity demands additional O₂ to working muscles, which is matched by O₂ delivery from the cardiopulmonary system to limit reliance on the less efficient anaerobic pathways. The kinetics, the maximal attainable (i.e. ) and the required for sub-maximal activities, are highly related to health, fitness and exercise performance [2,3]. Direct measurement of requires expensive, invasive and fragile instrumentation, such as metabolimeters [4]. As a consequence, the study of exercising is mostly confined to laboratories and clinics. During outdoor activities, wearing and carrying a metabolimeter can put the athletes and the instrumentation in danger. Therefore, estimating without reliance on a metabolimeter would be highly useful for a number of performance assessment applications.

Typically, when a metabolimeter is not available, the steady value of (i.e. ) is estimated from heart rate (HR). However, this methodology has limitations [5]. For example, for very low and very high intensity exercises, the relationship is not linear. Furthermore, heart rate is affected by a high day-to-day variability [6]. However, another method we might use to directly estimate is through its relationship with mechanical power output (P). Indeed, cycling exercise is a repetitive and easily testable activity in which the mechanical power output can be measured directly and reliably using a power meter [7] and even estimated using simple energetic relationships [8].

However, like heart rate, does not respond promptly to mechanical power output variations and dynamics must be taken into account [9]. The three distinct phases involved with dynamics include: 1) a cardio-dynamic phase-I, 2) a fundamental phase-II and 3) a slow phase-III component. Whilst phase-I and II are always present in responses to step-changes in the workload, the phase-III only becomes discernible at heavy and severe exercise intensities [10]. If the dynamic is considered to be linear (this assumption has been questioned multiple times [11–13]), a first-order model can be used to roughly approximate the at the next instant k+1 (i.e. (k+1)) from and mechanical power output (in Watt) at the current instant k (i.e. and P(k)) (Fig 1A): where G is the “gain” [14], is the resting [15] and Δt(k) is the time that separates the two instants k and k+1. This formulation has some shortcomings, e.g.: 1) changes of G and τ across exercise intensity domains [16] (or with transitions from greater baseline intensities [17]) and cadences (ω) [18], 2) prolonged exercise affects the relationship between G and P [19] and 3) response to exercise is affected by recent exercise history [20]. Such a description can be improved including those features known to be relevant or related to , e.g.: current and past values of mechanical power output, pedalling cadence, heart rate and respiratory frequency (RF).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1.

A schematization of the analytical equations approach (A) and artificial intelligence (AI) approach (B) is given. A: the current value of the power (P) and oxygen uptake () (i.e. P(k) and (k) respectively) are used in a formula (i.e. Δe^-t/τ, with Δ calculated as the difference between the steady state and the current ) to forecast future values of (i.e. ). B: In an AI approach to the time-series problem, current (k) and past values (k-1, k-2, … k-n) of heart rate (HR), P and cadence (ω) are used to forecast future values .

https://doi.org/10.1371/journal.pone.0229466.g001

The problem of forecasting data starting from observations of other variables taken sequentially in time can be considered as a time series prediction problem. Analytical equations of the dynamics have limited capacity to accurately model such complex data without requiring very large and complex formulations. An alternative approach may be found among the artificial intelligence technologies [21]. Particularly, machine learning algorithms can be used to learn from data and individuate patterns of variation between variables (Fig 1B). A machine learning algorithm that considers time sequences can be implemented by means of the so-called artificial neural networks, a biologically-inspired computational system that mathematically formalizes the connections among and within layers of artificial neurons [22]. Artificial neurons receive one or more inputs and sums them to provide an output. Inside the neuron, each input is weighted, and the sum is passed through a non-linear activation function. Given a sufficient number of layers and neurons, a neural network can always be trained (i.e. the weights of the neural network are calibrated) to approximate a real relation between inputs and outputs [23].

Examples of the application of artificial intelligence to time series problems include financial time series forecasting [24], as well as arrhythmia detection from ECG signals [25]. Recurrent neural network and, in particular long-short term memory [26], are suited for time series forecasting problems and sequences [27]. Unlike feed-forward neural network, recurrent neural networks make use of an internal memory to process sequence of inputs. This is a very important property when the prediction of the neural network depends on the historical context of inputs.

With respect to the field of cycling performance, for example, an artificial intelligence approach has been proposed for training data processing [28]. In the field of exercise physiology, a neural network has been developed [29] to model the heart rate versus mechanical power output relationship. With this neural network, it was possible to find the heart rate associated with the anaerobic threshold non-invasively in soccer players. In data estimation, Laitinen & Rasanen [30] used a neural network to estimate in children with congenital heart disease from inputs like heart rate and blood pressure. However, the accuracy attained suggested that the predictive power of their neural network was “insufficient” at that time. In 2017, Gonzalez et al., [31] presented an accurate mathematical description for dynamics during high-intensity variable cycling exercise, and the same authors suggested that a neural network could “perform even better” than their analytical model [32]. More recently, machine-learning has been used [33] to predict accurately during walking and with different daily activities [34], including cycling [35].

In light of the promises of the artificial intelligence technologies [21], the purpose of this pilot study was to predict the individual response of during high-intensity cycling exercise starting from easy-to-obtain inputs. We hypothesised that a recurrent neural network could provide accurate individualised predictions across a variety of exercise conditions. To highlight the potential of this methodology, we compared the predictive accuracy of the neural network with a first-order kinetics equation and a previously published higher-order model.

2 Methods

2.1 Experimental data

Seven recreational cyclists (6 males, Table 1) participated in the study and visited the laboratory on three separate occasions. The ethics committee of the Department of Neurological and Movement Sciences of the University of Verona approved the study.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Participants’ characteristics: Mean (SD) of the weight, the maximal oxygen uptake (), the peak power output (PPO) and the three intensity-levels adopted in the second and third tests (P₁, P₂ and P₃).

https://doi.org/10.1371/journal.pone.0229466.t001

The participants gave informed consent and the research was conducted in accordance with the declaration of Helsinki. All tests were performed on an electromagnetically-braked bicycle ergometer (Excalibur Sport, Lode). Measurements of mechanical power output and pedalling cadence were collected continuously. Respiratory measurements, such as and respiratory frequency were collected using breath-by-breath methods from an automated open-circuit gas analyser (Quark CPET, Cosmed). Immediately before every test session, the gas analyser and the flow meter were calibrated. Invalid breaths (i.e. those lying outside the ranges of respiratory frequency b/min 2–90 (min-max); ventilation (L) 0.100–10000, FeO₂% (%) 5–20; FeCO₂ (%) 1–10) were automatically removed in real-time by the CPET software. Heart rate was recorded continuously (beat-by-beat) during the test with a heart rate monitor incorporated into the gas analyser. Heart rate was interpolated and provided at the breath-by-breath time sequence by the metabolimeter.

During the first visit, participants underwent a graded exercise test (GXT) for aerobic assessment. and respiratory frequency data were averaged over 4 min at rest, to give a resting metabolic rate () and resting respiratory frequency (RF_R) respectively. Participants warmed up for 10 minutes at 85 W and freely chosen pedalling frequency. The graded exercise test started at a workload of 100 W for 4 min and, subsequently, the workload increased by 40 W every 4 min until exhaustion. The pedalling cadence during the test was kept constant at 90 rpm, using a monitor that provided participants with visual feedback. The peak power output (PPO) of the test was determined using the power of the last completed stage and the time of the last uncompleted stage [36]. The and the maximal respiratory frequency RF_MAX were defined as the highest value of and respiratory frequency registered during the test over a 20-s rolling average [37]. The first ventilatory threshold (VT1) was determined from visual inspection of: 1) the first disproportionate increase in minute ventilation (VE); 2) an increase in with no increase in (where is the exhaled volume of CO₂); 3) an increase in end-tidal O₂ with no consequent fall in end-tidal CO₂ tensions. The second ventilatory threshold (VT2) was determined from: 1) the second disproportionate increase in minute ventilation; 2) the first systematic increase in ; 3) the first systematic decrease in end-tidal CO₂ tension[38,39]. To account for the differences that exists in the versus power output relationships from graded versus constant exercise [9], the power values expected to elicit the associated with the first and second ventilatory thresholds were estimated using the equations established by Kuipers et al. [40]. We are aware that graded exercise testing protocols can influence the versus power output relationship [41], hence reducing the validity of the specific power values associated with first and the second ventilatory threshold, i.e.: P_VT1 and P_VT2, respectively. P_VT1 and P_VT2 were obtained considering the power output of the previously completed stage and the time of the completed stage when the ventilatory threshold occurred (e.g. 18 min = 240 W; 220 W (16 min) + 2 min/4min x 40 W).

Three different mechanical power output levels were defined for every participant as follows: moderate intensity P₁ = 0.5*P_VT1, heavy intensity P₂ = 0.5*(P_VT2-P_VT1)+ P_VT1 and severe intensity P₃ = 0.5*(PPO-P_VT2) + P_VT2 (Table 1). After a recovery period of 1 hour, participants performed a 30” Wingate test (WG) on a mechanically braked cycle ergometer (Ergomedic 894-Ea, Monark). During this test, the highest mechanical power P_MAX and the maximal cadence ω_MAX were determined.

During the second and third visit to the laboratory, athletes performed a warm-up for 10 min at a constant power of 85 W and rested for 4 minutes, before performing two different protocols on separate days. The first protocol (Test 1) consisted of a constant mechanical power output of 100 W for 4 minutes, followed by three repetitions of three constant bouts of 4 minutes at P₂, P₃ and P₁ (please see the additional material for a graphical representation of the protocol). The second protocol (Test 2) started with a linear increase in the mechanical power output (i.e. a ramp) from P₁ to P₃ in 4 min. The initial ramp was followed by: a 1-min bout at P₃, a 4-min bout at P₁, a 1-min ramp from P₁ to P₂, a constant bout of 3 min at P₂, a 4-min bout at P₁, a 2-min ramp from P₁ to P₂, a constant bout of 2 min at P₂ and a 4-min bout at P₁ (please see the additional material for a graphical representation of the protocol). These two arbitrary protocols were designed to elicit different dynamics behaviours and facilitate the convergence of the parameter estimation.

2.3 Neural network design

An artificial intelligence regressor was developed and used to predict values of . The neural network was coded and implemented using Python (ver. 3.6, Python Software Foundation), a high-level programming language for general-purpose programming. In particular, the open-source library called Keras was adopted to specifically design and test the neural network. The neural network was created using a Tensorflow backend with cuda support (2xNVidia GT750M i74xxx). A summary of the model is given in Table 2.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. A total of 21717 parameters have been included in the LSTM NN designed in this study.

Three LSTM layers of 32 neurons were used with 1 hidden layer of 10 feed-forward neurons and 1 output layer of 1 neuron. Input shape for LSTM layers were determined form the batch size (10), the number of past inputs considered in the time series (70) and the number of neurons of the layer.

https://doi.org/10.1371/journal.pone.0229466.t002

The neural network was composed with long-short term memory neurons [26], best suited for time-series analyses and sequence detection [27]. A total of 3 long-short term memory layers of 32 neurons each were formed, plus 1 hidden layer of 10 neurons and 1 output layer of 1 neuron. The neural network was trained with a stochastic gradient method (adagrad) that optimises a categorical cross entropy loss. The training dataset entries were shuffled and the whole dataset was crossed in 20 epochs. The weights were initialized as random values, while biases were initialised as random positive values. The batch size (that defines the number of samples propagated in the neural network) was set to 10.

There are no specific and scientifically proven steps to be followed in the design of the neural network. However, we know that the choice of the number of layers, the number of neurons, the number of epochs and the batch size affect the accuracy of the output and the computational time. Therefore, to select these parameters, we proceeded by trial and error, until the best combination of accuracy and computational time was found. The final architecture with 3 long-short term memory layers, has been shown to work well in other time-series classification problems using physiological data [43].

3 Results

Training the neural network required approximately 30 min (PC equipped with 2xTitan i75900), while Model 1 and 2 calibration (particle swarm optimization) required 10 min and 20 min respectively (MacBook Pro, 2.8 GHz Intel Core i7). Testing the models required few seconds for every simulation.

In trial 1, after the particle swarm optimization, mean values of the parameters of Model 1 were: G = 10.07 (0.85) mlO₂min^-1W^-1 and τ = 45 (3) s. Values for Model 2 are reported in the Appendix. For the neural network and Models 1 and 2, results of the residual and Bland-Altman analyses are presented collectively in Table 3 for both the experimental Trials. The performances of the neural network in Trial 1 and 2 are shown in Fig 2A and 2B for a representative participant.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. A) The performance of the regressor is shown for a single representative athlete.

Experimental data (circles) of oxygen uptake () are reported. Predicted values of (solid line) are superimposed on experimental data. In this example, MAE was 0.028 (i.e. 164 mlO₂min^-1), with a RMSE of 0.04 (i.e. 229 mlO₂min^-1). B) The performance of the regressor is shown for a single representative athlete during a WG test. Experimental data (circles) of oxygen uptake () are reported. Predicted values of (solid line) are superimposed on experimental data. In this example MAE was 0.03 (i.e. 176 mlO₂min^-1), with a RMSE of 0.05 (i.e. 294 mlO₂min^-1). Please see the S1 Material for the other individuals’ responses.

https://doi.org/10.1371/journal.pone.0229466.g002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Results are reported for the residual and the Bland-Altman analyses for the three different models (AI reg. (i.e. our AI regressor), model 1 (i.e. the first-order model) and 2 (i.e. the Gonzalez’s model)).

In Trial 1 we compared predicted and experimental data during a variable high-intensity exercise. In Trial 2 we compared predicted and experimental data during a brief 30” “all-out” Wingate test.

https://doi.org/10.1371/journal.pone.0229466.t003

The residual analysis in Trial 1 shows that the predictive power of the neural network was significantly superior to that of the other models, as seen by the smaller mean absolute error and root mean square error and higher correlation coefficient and variance explained (Table 3). For both our neural network and Models 1 and 2, the Bland-Altman analysis for the measured versus predicted showed no proportional error rate, with differences unrelated to the magnitude of the measurement error. In the case of the neural network the bias was not significant. For model 1, the equality line fell outside the confidence intervals of the mean bias of the sample and outside the limits of 200 mlO₂min^-1. Model 2 performed slightly better than Model 1: the equality line fell outside the confidence intervals of the mean bias of the sample but inside the limits of 200 mlO₂min^-1. For the neural network, the autocorrelation analysis suggested that there was no significant autocorrelation between observations and lagged observations. In fact, autocorrelation consistently stayed within the confidence bands. In the case of Model 1 and 2, the autocorrelation fell outside the confidence bands for the initial portion of the signal.

During Trial 2, residual analysis highlighted that the neural network could accurately predict the actual data both during the ascending and descending phases of the evolution. On the contrary, both Models 1 and 2 did not predict the additional required after high-intensity exercises. This is confirmed by the high values of correlation coefficient and variance explained for the predictions of the neural network (Table 3). Bland-Altman analysis suggested that, in the case of the regressor, the bias was not significant. On the other hand, for Model 1, the equality line fell outside the confidence intervals of the mean bias of the sample and outside the limits of 200 mlO₂min^-1. Model 2 performed slightly worse: the equality line fell outside the confidence intervals of the mean bias of the sample and outside the limits of 200 mlO₂min^-1. Bland-Altman analysis suggested that, in the case of Models 1 and 2, the biases were significant. The autocorrelation analysis for the predicted values of the neural network showed that there was no significant autocorrelation between observations and lagged observations. When predictions were made with Models 1 and 2, autocorrelation analysis highlighted that the other models failed to incorporate important relationships between current and past input values. In fact, a consistent portion of the autocorrelation lied outside the confidence bands.

4 Discussion

We hypothesized that a recurrent neural network approach could be successfully used to accurately predict individual cycling data from easy-to-collect inputs [21]. In fact, the mechanical power output and the pedalling cadence are both easily collectable by portable power-meters [7]). Indeed, heart rate and respiratory frequency are both measurable with chest belts [47,48] and have already been successfully used by Beltrame et. al. [49] for the prediction of from wearable sensors.

The ability of neural networks to model complex data was already known and other more basic learners could have been used (e.g. k-nearest-neighbour or support vector machine). While simpler learners like Hammerstein-Wiener models have been already tested [32,50], we are not aware of any existing example of the application of k-nearest-neighbour or support vector machine methods in the prediction of during high-intensity cycling.

However, the major innovation of our neural network lies in the long histories of values of the inputs (reflected in the number of input neurons). Latinen & Rasanen [30] used a neural network with 14 input neurons, one hidden layer of 4 neurons and one output neuron for . Beltrame et al. [33] used a neural network with 7 input neurons, one hidden layer of 11 neurons and one output neuron. Both studies only used current inputs and not past values. Beltrame et al. [49] used only 1 sec of lag to include “dynamic changes” of heart rate. Very recently, Borror et. al. [35] presented a neural network that can predict cycling in cycling, with a workflow similar to ours. They included body mass, mechanical power output, pedalling cadence and heart rate as inputs. However, in their work, no past input values are passed to the neural network, and the heart rate dynamics is only considered by means of its “time derivative”. In our neural network, there were 3 hidden layers of 32 long-short term memory neurons each, one hidden layer of 10 neurons and one output neuron. The neurons adopted were long-short term memory neurons, particularly suited for time-series analysis [26]. To the best of our knowledge, we are the first to apply recurrent neural networks to the prediction of cycling .

The predictive power of the neural network was very high during Trial 1, as measured and predicted showed a nearly perfect agreement (MAE = 229 (35) mlO₂min^-1, r = 0.94). The performances of models 1 and 2, although inferior, were still good in Trial 1 (MAE = 355 (86) mlO₂min^-1, r = 0.83 and 273 (49) mlO₂min^-1, r = 0.9 for model 1 and 2 respectively). This means that the performances of our neural network and Models 1 and 2 were very close. However, a robust model should predict data across a wide range of scenarios. To this end, we tested our models using a short “all-out” sprint effort (Wingate test). Importantly, Models 1 and 2 were not designed specifically for the Wingate test and have a limited number of parameters that can be tuned. However, the neural network, due to the considerable number of parameters used, can work well across a wider range of exercising scenarios. In fact, the number of parameters used may limit the number of physiological mechanisms that can be mathematically described.

In Model 1, a single phase is included and characterised by the parameter τ. The time-constant τ and the oxygen gain G, in Model 1, are constant across all the exercising intensities. Therefore, it becomes difficult to predict experimental values of during brief “all-out” exercises [51]. In Model 2, the parameter T_I (see Appendix) has been included to account for the delayed component that sum to the principal element. However, if the time duration of the exercise is shorter than T_I (e.g. the Wingate test lasts 30”), then this additional component is not activated. In Model 2, the high number of parameters affected the confidence of parameter estimation and this is confirmed by the large variability of the parameter estimates. Mean and standard deviation of the variables found with our experimental data, were compatible to those reported in the original article [31] (Appendix).

On one hand, the high predictive power of the neural network, although reduced, was remarkably conserved during Wingate test (Trial 2: MAE = 304(150) mlO₂min^-1, r = 0.89). This means that the dataset that we used to train the neural network (in terms of duration of the exercises) for every participant, was large enough to provide good predictive power. Further studies are needed to establish the minimal amount of data that should be used to train a neural network and retain a high predictive ability. On the other hand, the performances of Models 1 and 2 deteriorated during Wingate test (Trial 2: MAE = 391(71) mlO₂min^-1, r = 0.75 and MAE = 463(112) mlO₂min^-1, r = 0.48 for the Model 1 and 2 respectively). It can be noticed (Fig 2B) that a small lag is present between the measurements and predictions. This might be because two of the inputs used by the neural network (i.e. respiratory frequency and heart rate) did not promptly react to abrupt changes in power output. However, an autocorrelation analysis showed that our neural network could incorporate relevant relationships between current and past input values, whereas Models 1 and 2 could not. Due to the reduced number of parameters of Model 1 and 2, the predictive power does not heavily depend on the amount of data used to calibrate the parameters. We suggest that the performance of the analytical models, although inferior, is guaranteed even if smaller datasets are used for their calibration. We did not investigate the influence of the dimension of the training set on the performance of the neural network, but we believe that the performance would deteriorate with smaller and smaller training datasets. This is a first limitation of a neural network approach: we rely on large datasets.

The second main limitation of the neural network method lies in its “black box” approach. In fact, it is unlikely we can understand how the non-linearities of the dynamics are represented inside the neural network. Additionally, our exercises were carried out in the laboratory environment and they were limited in time (max duration ~1400 s). In practical settings (e.g. training and races), a cardiovascular drift could mislead our estimations. The use of long-short term memory neurons makes it difficult to understand the variables that contribute the most to the total estimation. In our study, the pedalling cadence was kept constant, so it is likely that the contribution of this variable may be limited in our study. As well, respiratory frequency indicating the ventilatory response to exercise has an important link with , while heart rate has additional known associations with exercising [5,52].

Even though we investigated a few different exercising conditions (i.e. moderate, heavy and severe intensity and “all-out” efforts), we should the results of this pilot study with caution. In fact, more work is needed before this algorithm could be embedded in a portable system able to estimate cycling in real-time: the verification of the predictive ability of the neural network on a larger sample (7 cyclists is a small sample) and on different environmental conditions (e.g. outdoor). Also, including input parameters like body mass, gender and fitness level, may provide in the future even better predictive outcomes for the estimation of the aerobic performance. Importantly, the ability of the neural network in predicting the values for an individual who was not included in the training dataset, has yet to be assessed.

5 Conclusions

In the context of forecasting values, the results of our pilot study suggested that a recurrent neural network can use great quantities of information from other mechanical (such as mechanical power output and pedalling cadence) and physiological markers (such as heart rate and respiratory frequency), as well as past input values, to attain accurate predictions of cycling . Results suggest that this algorithm has the potential to be, in the next future, embedded in a portable system and to provide real-time assessment of individual cycling during training or racing.

Appendix

Dynamics equations of Model 2 (see [31,53]) are reported with a formulation that can be readily implemented in a spreadsheet. The main difference between this model and Model 1 is that phase-II (“fast” phase) and III (“slow” phase) of dynamics are considered in this model. Gonzalez et al. included these two additional phases with two delayed components that become active only after a given period. The principal governing equation is:

Where is the principal component that is active after T_II and is characterized by a time-constant τ_II.

Where A_II(k) can be computed as:

Where s is the gain for the fast phase. is the slow component of that activates after T_III and is characterized by a time constant τ_III.

Where A_III(k) can be computed as:

Where Pc is a “critical power” threshold.

In trial 1, after PSO, the values of the parameters were (mean(SD)): = 397(398) mlO₂min^-1, Pc = 359(39) W, Δ = 79(72) W, s = 8.67(0.49) mlO₂min^-1W^-1, τ_I = 43(1.38) s, τ_II = 199(52) s, T_I = 10(6.72) s, T_II = 113(27) s. In trial 2, after PSO, the values of the parameters were: = 779(445) mlO₂min^-1, Pc = 383(15) W, Δ = 64(34) W, s = 9.03(0.9) mlO₂min^-1W^-1, τ_I = 42(1.7) s, τ_II = 183(30) s, T_I = 11(4.5) s, T_II = 103(34) s.

Supporting information