Development and validation of dynamic bioenergetic model for intermittent ergometer cycling

Purpose The aim of this study was to develop and validate a bioenergetic model describing the dynamic behavior of the alactic, lactic, and aerobic metabolic energy supply systems as well as different sources of the total metabolic energy demand. Methods The bioenergetic supply model consisted of terms for the alactic, lactic, and aerobic system metabolic rates while the demand model consisted of terms for the corresponding metabolic rates of principal cycling work, pulmonary ventilation, and accumulated metabolites. The bioenergetic model was formulated as a system of differential equations and model parameters were estimated by a non-linear grey-box approach, utilizing power output and aerobic metabolic rate (MRae) data from fourteen cyclists performing an experimental trial (P2) on a cycle ergometer. Validity was assessed by comparing model simulation and measurements on a similar follow-up experimental trial (P3). Results The root mean square error between modelled and measured MRae was 61.9 ± 7.9 W and 79.2 ± 30.5 W for P2 and P3, respectively. The corresponding mean absolute percentage error was 8.6 ± 1.5% and 10.6 ± 3.3% for P2 and P3, respectively. Conclusion The validation of the model showed excellent overall agreement between measured and modeled MRae during intermittent cycling by well-trained male cyclist. However, the standard deviation was 38.5% of the average root mean square error for P3, indicating not as good reliability.


Introduction
The skeletal muscles responsible for human locomotion are dependent on energy released from adenosine triphosphate, which can be replenished by the anaerobic alactic and lactic systems and the aerobic system, respectively.Due to the different rates and capacities of these systems, the alactic system is the main source of adenosine triphosphate during short sprints, whereas the lactic system is most prominent during longer sprints, and the aerobic system is the main contributor during prolonged exertion.Despite this general pattern, research shows that the anaerobic system is also important in endurance sports characterized by varying intensity, such as cycling and cross-country skiing (Skiba et al. 2012;Gløersen et al. 2020).Over varying terrain, the anaerobic system is utilized during periods of greater exertion, such as uphill sections, where the aerobic system cannot meet the metabolic demand, resulting in oxygen deficits that are repaid during periods of lower exertion (Gløersen et al. 2020).The anaerobic system is also used as a buffer during the ramp-up of the aerobic system, which has a much slower response.The alactic system is limited to short contributions due to the limited amount of available phosphocreatine, but it can also be rapidly recovered during periods of lower intensity to be utilized again (McCann et al. 1995).In theory, the glycogen stored in skeletal muscle cells could fuel the lactic system for a considerable amount of time.However, during high intensity exercise glycolysis causes metabolites to accumulate, which can disturb the homeostasis of the muscles over time, ultimately limiting further glycolysis and therefore greatly reducing the work capacity of the muscles (Stanley and Connett 1991).
Numerous attempts to mathematically describe the behavior of the human energy supply systems have been made.Morton and Billat (2004) presented a model applicable to intermittent exercise, which was further developed by Skiba et al. (2012Skiba et al. ( , 2015)).These models allow a finite energy store (comparable to anaerobic capacity) to be utilized and recovered, with the utilization and recovery rates dependent on metabolic demand.However, the models make no distinction between energy generated by the lactic and alactic systems, nor do they factor in the dynamic behavior of the aerobic system.Moreover, according to these models, energy from the aerobic system is provided at the critical power level at all times.
Regarding the dynamic behavior of the aerobic system, it is widely accepted it can be described as the sum of a constant baseline and three exponential functions with their own separate error signals, time constants, and time delays (Poole and Jones 2012).The first exponential function (the cardio dynamic component) is of small magnitude and duration and is therefore often incorporated into the second exponential function (the primary component).For the third exponential function (the slow component), several underlying mechanisms have been suggested, e.g., a drift in metabolic demand, reduced efficiency in recruited muscle fibers, and successive recruitment of less efficient muscle fibers (Poole and Jones 2012).
Several mathematical models that incorporate the described dynamic behavior of the aerobic system have been suggested.A model proposed by Artiga Gonzalez et al. (2019) gives a mean absolute percentage error of 2.1% when fitted to a race-like road cycling exercise, but this increases to 10.8% when the model was fitted to another semi-racelike protocol and then applied to the race-like exercise.Furthermore, the model relies wholly on power output measured by an ergometer, and as such any deviation from a purely linear relationship between power output and oxygen uptake is captured only by adjustments to oxygen dynamics modeling.
In reality, purely linear demand is unlikely.For instance, there is likely to be additional metabolic demand related to the acidification of the muscles (Gaesser and Brooks 1984).As a result, the model proposed by Artiga Gonzalez et al. (2019) gives no further insight into the underlying control mechanisms of oxygen uptake, nor can it provide estimates of anaerobic metabolic rates.Another model proposed by Stirling et al. (2005) offers the possibility of calculating both oxygen demand and oxygen deficit (e.g., anaerobic contribution) but no way of separating lactic and alactic contributions.In terms of oxygen uptake response and underlying oxygen demand, this curve-fitting model provides no additional insights into the physiological mechanisms of oxygen uptake dynamics.Moxnes et al. (2014) tested this model for several constant power outputs with subsequent recovery and showed poor results at higher exercise intensities.The study also tested three additional models describing the dynamics of the aerobic system.The most complex model agreed well with measured aerobic metabolic rate although it did not capture the exact dynamics.However, the model used the lactic metabolic rate as input, approximated from measured blood lactate.This is a crude estimation and since lactic metabolic rate is used as input, the model offers no additional information on the dynamics of the lactic system nor the alactic system.
Several attempts to collectively describe the dynamics of all three energy supply systems and their interactions mathematically have been made based on the hydraulic tank approach proposed by Margaria (1976).In his model, the available energy from each of the energy supply systems is depicted as fluid in a tank and the utilization of energy as fluid flow from each of the tanks.The relative size and positioning of the tanks and the flow capacity of their interconnecting pipes define the dynamic behavior of the system.Morton (1986Morton ( , 1990) ) further generalized this approach and derived restrictions based on empirical data.Sundström (2016) further developed the hydraulic tank approach with separate supplies of fat and carbohydrates and found it to show good agreement regarding relative aerobic contribution but overestimations of times to exhaustion.Lidar et al. ( 2021) calibrated Morton's model, and a similar model with a combined anaerobic system, using individual data from simulated sprint cross-country skiing races and achieved good overall agreement.However, neither of the models were able to fully capture the dynamic behavior of the aerobic system.Moreover, the protocols used did not include any periods of significant recovery.Had this been the case, anaerobic recovery would likely have been exaggerated by the model.Firstly, because all excess aerobic energy, in relation to metabolic demand, is added to the anaerobic stores, while in reality, recovery efficiency of 100% is highly unlikely.Secondly, there is likely to be residual metabolic demand during passive (and possibly active) recovery following heavy exertion (Gaesser and Brooks 1984) which was not included in the speed-dependent linear relationship for metabolic demand used by Lidar et al. (2021).
Thus, to better understand metabolic processes during variable-intensity endurance exercise and predict human performance in endurance exercise, a robust mathematical model is needed-one that describes the interaction between metabolic energy demands and the three metabolic energy supply systems.Therefore, the aims of this study were (1) to develop a bioenergetic model describing both metabolic energy demands and the aerobic and anaerobic lactic and alactic metabolic rates from readily measurable quantities during intermittent cycling and (2) to assess the validity and reliability of the developed model based on individual aerobic metabolic rate.

Method
A novel bioenergetic supply and demand model was developed from the bottom up.To assess the validity and reliability of the model, athletes performed cycle ergometer exercises while following purposive variable-power protocols.One protocol was used to individually adapt the bioenergetic model using non-linear grey-box parameter estimation.The individualized models were then applied to data from another similar protocol.

Experimental protocols
Each subject completed power-controlled tests on three separate days on a cycle ergometer with 2-7 days in between test days (there were 13 days between test day 2 and 3 for one subject due to illness).The first day comprised a submaximal incremental protocol (P1a) and a maximal incremental protocol (P1b).The second and third days each consisted of an intermittent protocol (P2 and P3, respectively) with individualized constant ergometer power outputs.All subjects could view their power output, heartrate, and cadence on the ergometer monitor during the tests and were instructed to maintain a cadence of 90 rpm throughout all protocols.

Incremental protocols and individual power output calculations
Power output in P1a started at 80 W and increased by 20 W every 3 min until monitored ventilatory data indicated that a respiratory quotient (RQ) above 1.0 had been reached.A venous lactate sample was taken 30 s before each power output increment.Participants rested for about five minutes between P1a and P1b.Power output in P1b started at 100 W and increased 40 W every minute until volitional exhaustion (subject stopped cycling or cadence dropped below 70 rpm for more than 3 s).Participants were verbally encouraged during the last few minutes of the protocol.
From P1a, the venous blood lactate concentration ([bLa − ]) vs power output was fitted to a third degree polynomial.A blood lactate baseline was established by calculating the mean of the first two values of [bLa − ], then adding one additional value at a time to the mean calculation until the difference between the mean value and the next [bLa − ] value exceeded 0.2 mmol L −1 .This mean value was used as the baseline for [bLa − ] and the lactate threshold (LT1) was the value 0.2 mmol•L −1 above the baseline.This [bLa − ] difference of 0.2 mmol•L −1 was the rounded-up value of the typical error (0.12 mmol•L −1 ), determined with a selection of duplicate samples (n = 40) from the beginning and end of the study.The power output at LT1 (P LT1 ) and onset of blood lactate (P OBLA ) were taken from the third degree polynomial at the LT1 [bLa − ] and [bLa − ] = 4 mmol•L −1 , respectively.
Peak oxygen uptake ( VO 2,peak ) for the individual power output calculations was calculated as the maximal 30 s mean during P1b.Spirometry data from the submaximal test was used to calculate MR ae according to where RQ is the nonprotein respiratory quotient (McArdle et al. 2009), VO 2 is the oxygen uptake in L•min −1 , and 4184 60 is the transformation from kCal•min −1 to W. A linear regression of VO 2 vs MR ae for the submaximal power outputs in P1a (RQ ≤ 1) was extrapolated to VO 2,peak in P1b and used to extract the corresponding power output (P VO2,peak ).Finally, the maximum power output (P max ) was calculated as the power output at exhaustion during P1b minus 40 W times the fraction of the minute left until the next power increment (e.g., exhaustion with 15 s cycling left at 320 W would give P max = 320 − 40 15 60 = 310 W).The individual values of P LT1 , P OBLA , P VO2,peak , and P max were used to decide the individual power output levels for P2 and P3.

Intermittent protocols
Seven power output levels for P2 and P3 were purposively chosen to engage the three metabolic energy supply systems to various extents.The power output levels are shown in Table 1.The designs of P2 and P3 are shown in Fig. 1.All exercise and rest periods in the protocols were of specific, predefined durations, except the last periods of L4.2 and L4.3, respectively, which were continued until volitional exhaustion (subject stopped cycling or cadence dropped below 70 rpm for more than 3 s).Participants were verbally encouraged during the last few minutes prior to exhaustion.Participants were allowed to remove the spirometry mouthpiece for a short time (< 30 s) during active rest at two specified instances during each of the intermittent tests but were otherwise instructed to keep the mouthpiece in throughout the test (also at the point of exhaustion).Additionally, two subjects removed the mouthpiece for < 20 s at exhaustion.

Bioenergetic model description
The developed bioenergetic model mathematically describes the behavior of the alactic, lactic, and aerobic energy supply systems.To enable this, it also models the underlying energy demand rates and a representation of active muscle lactate concentration ([mLa − ]).The lactic system is the part of the glycolysis that causes an accumulation of lactate from the pyruvate that is not oxidized in the aerobic system.The model is given by Eqs. ( 2)-( 13) as a state-space model where u denotes input (independent variables), x denotes state (mediating variables), and y denotes output (dependent variable).The bioenergetic model must adhere to the principle of energy conservation described as where MR dem is the metabolic demand rate and MR sup is the metabolic supply rate.

Energy supply
The energy supply of the bioenergetic model can be expressed, in terms of metabolic rate, as where MR al = ̇x1 is the alactic system metabolic rate, MR la = x 2 is the lactic system metabolic rate, and is the aerobic system metabolic rate.MR rest is both the metabolic demand and the aerobic metabolic rate at sitting rest, but without the metabolic demand due to ventilation, which is treated separately.x 3 is the primary component of aerobic metabolic rate.The bioenergetic system metabolic rates ( ̇x1 , x 2 , x 3 ) and [mLa − ] are governed by the following system of time-dependent differential equations: The rate of change of the aerobic metabolic rate ̇x3 is con- trolled by an error signal MR dem − MR rest − x 3 , and a time constant τ ae determining the response time.The error signal is the difference between the total metabolic demand rate MR dem and the current aerobic metabolic rate MR rest + x 3.In a dynamic model formulation, this is equivalent to a monoexponential response without time delay.We hypothesize that this mono-exponential in combination with drifting metabolic demand and high metabolic demands will capture the complete dynamic behavior of the aerobic system. (3) To our knowledge, the only proposed models of the lactic system are those within the hydraulic tank models that include a separate lactic system.Given the close relation between the lactic and aerobic systems, we propose that the rate of change of the lactic system ẋ2 is similarly regulated by an error signal (MR dem − MR rest − x 3 − x 2 ) and a time constant τ la .Here, the error signal is the difference between the total metabolic demand rate and the sum of the current aerobic metabolic rate and current lactic metabolic rate.Additionally, the response of the lactic system is dampened by the factor K times the rate of change of aerobic metabolic rate.
Furthermore, in Eq. ( 5), ̇x1 is the alactic metabolic rate normed with the alactic energy capacity E al,max .Since the regulation of the alactic system acts directly on the alactic metabolic rate, the alactic system will instantly meet the metabolic demand that is not supplied by the lactic and aerobic systems.By norming with E al,max , the state variable x 1 will have a feasible range of 0-1, where 0 equates to a fully recovered alactic system and 1 equates to completely depleted alactic system.
Finally, ̇x4 is the rate of change of [mLa − ] normed with the muscle lactate storage capacity V m , which gives the state variable x 4 a feasible range of 0-1, where 0 equates to no accumulated muscle lactate (i.e., fully recovered lactic system) and 1 equates to maximum accumulation of muscle lactate (i.e., further utilization of the lactic system is limited).Muscle lactate is accumulated at a rate of x 2 /V m and reduced by a maximum amplitude A red /V m times the [mLa − ] x 4 times a demand factor Z dem , given by Eq. ( 6).The demand factor is equal to 1 when MR dem is equal to MR rest and decreases linearly to 0 when MR dem equals the metabolic rate at the lactate threshold MR lt .This means that an increased [mLa − ] as well as a decreased MR dem will give a faster rate of reduction in muscle lactate.Furthermore, to capture the behavior of an elevated steady state [bLa − ] (and hence [mLa − ]), the reduction of [mLa − ] will cease from the lactate threshold (LT1).Since MR ae may reach an elevated steady state above LT1 and since the lactic system in the proposed model describes the part of glycolysis resulting in lactate accumulation, this would cause the lactic glycolysis to cease.Hence, muscle lactate removal above LT1 would not result in a steady state [mLa − ].
The given equations are supplemented by several conditions.Should the alactic metabolic rate ẋ1 fall below zero, only part of the energy diverted to the alactic recovery is converted into usable alactic energy due to the inevitable hysteresis cost of recovery.This is expressed in Eq. ( 7), where η al is the efficiency of alactic recovery and thus the hysteresis cost is (1 − al )• ̇x1 .
Recovery of the lactic system is linked to the reduction of accumulated muscle lactate given in Eq. ( 5).The anaerobic glycolysis itself is not allowed to fall below 0, as stated by Eq. ( 8).
Finally, the MR ae cannot exceed MR ae,max (the energetic equivalent of VO 2max ), hence x 3 is limited by MR ae,max − MR rest , as given by Eq. ( 9).

Energy demand
The remaining model equations define the components of metabolic demand.Adding the metabolic demands together, including MR rest , the total metabolic demand is given by Eq. ( 10).
The fundamental metabolic demand rate MR f is modelled by Eq. ( 11) as a linear function of the measured power output from the ergometer u 1 .( 6) in theory, is the metabolic demand rate of unloaded cycling and B f the increase in metabolic demand per unit increase in power output.Additionally, MR f is set to 0 when the power output is 0, i.e., during passive rest.The metabolic demand due to ventilation is not included in MR f .Hence, MR f is different from the often-used linear relationship between metabolic rate and power output that can be established at submaximal workloads.
The metabolic demand rate due to ventilation MR ve has been shown to vary non-linearly with minute ventilation but with large inter-individual differences (Vella et al. 2006).Here, it is assumed that a single relationship between MR ve and the measured minute ventilation u 2 is valid for all measured u 2 and that in theory u 2 = 0 would result in MR ve = 0. Hence, MR ve is given by Eq. ( 12), where A ve is the maximum metabolic demand rate due to ventilation, VE,max is the maxi- mum measured ventilation during P2, and B ve is a distributing factor between a linear and a quadratic response.Both u 2 and VE,max are calculated as body temperature, pressure, water vapor saturated (BTPS).
The assumption that Eq. ( 12) is valid for all measured u 2 is the reason metabolic demand due to ventilation is not included in MR f and MR rest .
Metabolic demand associated with the accumulation of metabolites in the muscles MR acc is given by Eq. ( 13).Though the accumulation of lactate itself may be the cause of a rise in metabolic demand due to lactate shuttling, the [mLa − ] in this model is used as a proxy for the concentration of various accumulated metabolites that may cause additional metabolic demands (Gaesser and Brooks 1984).The metabolic demand associated with accumulated metabolites is given the same form as the metabolic rate from ventilation.
A acc is the maximum amplitude of metabolic demand from accumulated metabolites and B acc is the distributing factor between a linear and a quadratic response.

Model individualization through parameter estimation
In total, the described bioenergetic model includes 17 parameters that can be adjusted to individually adapt the model.Of these, 14 were individually estimated using the non-linear grey-box parameter estimation algorithm ( 11) (nlgreyest) in the Systems Identification Toolbox in MAT-LAB (R2020b, Mathworks Inc., Natick, MA, United States), based on Lidar et al. (2021).This algorithm minimizes the mean squared error between measured and model data by iteratively adjusting the model parameters and solving the system of ordinary differential equations in the bioenergetic model (Eq.4).To prevent feasible-range violation of the x 1 and x 2 state variables, a multiplicative penalty function was applied to the dependent variable (MR ae ).If either of these state variables came close to its limit, MR ae would rise rapidly, causing a higher error between measured and modeled data.To enforce a [mLa − ] close to the maximum of 1 at exhaustion, [mLa − ] was forced to exceed 0.85, 100 s prior to exhaustion during P2.A multi-start algorithm was implemented to increase the likelihood of global optima convergence, i.e., multiple initial values were set for selected parameters (see Table 2) and the parameter estimation algorithm was executed with all combinations of the initial values.The set of parameters that rendered the smallest mean squared error from the multi-start parameter estimation was considered the global optimal solution.
Three of the 17 parameters were treated separately.The maximum minute ventilation VE,max was considered an individual constant with values set as the maximum measured VE during P2 for each subject.The alactic system is the buffer of the bioenergetic system and, as such, it does not affect the behavior of the lactic or aerobic systems, nor metabolic demand.Therefore, the alactic capacity E al,max was adjusted after the parameter estimation to the initial value times the maximum achieved value of x 1 during P2.Finally, the smoothing factor K was set to a constant value of 0.8 since preliminary parameter estimations showed that lower and higher values resulted in unreasonably low and high contributions from the lactic system respectively.
To avoid infeasible parameter estimations, upper and lower limits were set for the 14 parameters in the parameter estimation.For MR rest , MR lt , and MR ae,max , the initial values (Table 2) were approximated from P2 measurements (except MR lt from P1a).Since the approximation of MR ae,max was the measured maximal 60 s mean during P2, it was considered a reasonable but underestimated approximation, and hence was limited to a maximum increase of 5% (no allowed decrease).Preliminary parameter estimations yielded unreasonably low and high values of MR rest , but with absolute limits preventing these extreme values there were only minor deviations from the initial values.MR lt was limited to a variation between 70% and 130% of its initial value.
The aerobic time constant τ ae was limited to a variation between 10 and 100 s (di Prampero 1981) with an initial estimate of 25 s (Artiga Gonzalez et al. 2019).The initial value of τ la was established from Fig. 5 in Gastin (2001), in which energy demand is ~ 73 ml O 2 equivalent•kg −1 •min −1 .If glycolysis were to continue toward this value without any contribution from the aerobic system, it would attain 63% after one time constant, which is about 46 ml O 2 equivalent•kg −1 •min −1 .Sketching this hypothetical curve, the value of 46 ml O 2 equivalent•kg −1 •min −1 would be reached after approximately 12.5 s.Large variations from this approximation yielded unreasonable results in preliminary parameter estimations, and hence τ la was limited to between 10 and 15 s.
The initial values of A ve were set to 8.8% of the approximated MR ae,max according to the mean oxygen cost of ventilation in the study by Vella et al. (2006).The allowed range was set to 5-15% of MR ae,max , where the lower limit was the lowest value according to Vella et al. (2006), but the upper limit was the approximate maximum cost of ventilation found by Aaron et al. (1992), since the highest value from Vella et al. (2006) stood out compared to similar studies.It was assumed that the imposed metabolic demand related to accumulated metabolites would be less than that related to ventilation.As such, initial values of A acc were set to 3.8% or 4.2% of MR ae,max (about half of A ve ), but with the same limits as A ve , since no reasonable indications of its magnitude could be found in the literature.Initial values of A f , B f , B ve , B acc , A red , η al , and V m , were all chosen from preliminary parameter estimation results.The limits for A c and B c were applied to avoid unreasonable results, but even so, no values outside these ranges were obtained during the preliminary parameter estimations.The lower and upper limits of B ve and B acc were set to 0 and 1, respectively, since these parameters give the distribution between a linear and quadratic response in their respective equations.The limits of A red and V m were both set as percentages of V m .Since the values of these parameters were uncertain, A red was allowed to vary more than ± 30% and V m was allowed to vary ± 50% of their respective initial values.
All parameters' initial values, estimated values, and allowed ranges are summarized in Table 2.

Data collection
All tests were performed on a belt-braked cycle ergometer (Monark LC7 TT Novo) with the power output and cadence sampled at 1 Hz.The power output was smoothed by setting the power output of each period to be equal to the mean power output for the same period.To reduce the number of data points and hence the calculation time for parameter estimation, the power output was resampled to 0.5 Hz using linear interpolation.An example of the measured and smoothed power output is shown in Fig. 2. All subjects used their own personal cycling shoes and were allowed to use their own pedals.Prior to the first protocol, all subjects adjusted saddle height and handlebar position (height and forward/backward distance) to their preference.These measurements were recorded and used for all subsequent test protocols.
Venous blood samples for analyses of [bLa − ] were taken using 1-ml syringes from a 10 cm extension set connected to a catheter (Discofix® C; Vasofix® Safety; 1.1 × 25 mm, B. Braun Melsungen AG, Melsungen, Germany) in the cephalic vein or, in some cases, the mediana cubiti vein.Between the samples, the system was flushed with isotonic saline to avoid coagulation.Thus, each sampling started with discharging a volume greater than 2 ml before the actual sample was taken.The [bLa − ] samples were analyzed using laboratory equipment (Biosen S-line; EKF-Diagnostic, Cardiff, UK) within 10 min of completion of P1b, P2, and P3 respectively.
The subjects' breath-by-breath expiratory minute ventilation and fractions of expired oxygen and carbon dioxide, were measured using an automated metabolic measurement system (Moxus Modular Metabolic System; AEI Technologies Inc., Pittsburg, USA) described in a previous paper (Ainegren et al. 2018).For the individual power output calculations, 60-s averages of the VO 2 and RQ reported by the system were used.For the parameter estimation and model performance assessment, data from the metabolic measurement system was preprocessed in MATLAB.The inspiratory minute ventilation was calculated from the expiratory minute ventilation and fractions of inspired and expired nitrogen and expired oxygen and carbon dioxide using the Haldane transformation.The minute ventilations and gas fractions were smoothed using locally weighted linear regression (smoothdata with lowess method) with a 40-s time window to reduce the noise and obtain the underlying mean curve (Cleveland 1979).Gaps in the data during periods in which the participants removed the spirometry mouthpiece were filled using autoregressive modeling (fillgaps).From this smoothed data, VO 2 , carbon dioxide output, and RQ were calculated.Finally, MR ae was calculated from VO 2 and RQ according to Eq. ( 1).Expiratory minute venti- lation as BTPS (u 1 ) was used as an input for the calculation of MR ve , since it conveys actual expired air with current air temperature and density, but for all other calculations standard temperature, pressure, dry air (STPD) was used.
To remove errors introduced by smoothing at the beginning and the end of the protocols, the first 100 s and last 60 s of the data was omitted and the data was resampled to 0.5 Hz using linear interpolation.This was also done to reduce the number of data points and hence the calculation time for parameter estimation.The first 100 s contained about half of the initial sitting rest period and the last 60 s contained a third of the concluding period (L2 for P2 and AR for P3).The breath-by-breath VO 2 and expiratory minute ventilation as STPD ( VE,STPD ), as reported by the metabolic measurement system, is compared to the smoothed data for one subject in Fig. 3.The individualized bioenergetic models were used for numerical simulation of P2 and P3.In those simulations y (MR ae ), MR la , MR al , x 1 , MR dem , MR f , MR ve , MR acc , and x 4 ([mLa − ]) were calculated with measured power output and minute ventilation as BTPS as inputs.

Statistical analysis
Data are presented as mean ± SD unless otherwise specified.As a measure of the overall model agreement with measurements, the root mean squared error (RMSE) in W between measured and modeled MR ae was calculated with the timeresolved data (0.5 Hz) for each subject.Using the same data, the mean absolute percentage error (MAPE) was calculated with measurement data as reference and used as an overall measure of the model's agreement with the measurements in relative terms.
The calculated mean ± SD of RMSE and MAPE are considered metrics of the validity, while SD of RMSE and MAPE alone are considered metrics of the reliability.Data from P2 is considered to show the validity and reliability reflecting the combined error of the model formulation, parameter estimation, data acquisition, and data preprocessing, while data from P3 also includes errors from day-to-day variation and protocol differences.The method of parameter estimation has been applied in several studies (Artiga Gonzalez et al. 2019;Lidar et al. 2021) and is believed to find close-to-optimal solutions.Also, the validity of the metabolic measurement system has been tested (Beltrami et al. 2014).Hence, data from P2 is believed to primarily to show validity and reliability of the model formulation and data preprocessing, and the difference between P2 and P3 is believed to primarily highlight the day-to-day variation and protocol specificity.
To further evaluate the model-to-measurement agreement, Bland-Altman plots (Bland and Altman 1999) were drawn with all data points from the time-resolved data (0.5 Hz) for P2 and P3 respectively.This identifies the mean difference ± 95% limits of agreement of MR ae .
Statistical parametric mapping (SPM) has been used for statistical inference in neuroimaging (Friston et al. 2007), but has also been shown applicable to one-dimensional data (Pataky 2010).In the case of an SPM t-test; residuals, variances, and finally t-statistics are calculated from a general linear model for each node in the data (Pataky 2010).Since one-dimensional data (e.g., time series) usually exhibit spatiotemporal smoothness, a corrected critical t-value is established from the geometry of the data and the estimated level of smoothness using random field theory (Pataky 2010, Appendix A).In this study, SPM was used to perform twotailed paired t-tests on the measured and modeled MR ae normalized with the individual parameter MR ae,max for the entire time series for P2 (1379 nodes) and P3 (1079 nodes), respectively.Calculations were made with the open-source SPM 1D software (Pataky 2012) in MATLAB.An alpha level of < 0.05 was chosen with a null hypothesis of there being no difference between measured and modeled data.For the statistical inference, p-values are reported for each cluster of adjacent nodes where the t-statistics exceed the established critical t-value for P2 and P3, respectively.The results from the SPM analysis were also used to investigate if any specific parts of the intermittent protocols show signs of greater disagreement.

Results
The RMSE of the individual MR ae between measured and modeled data was 61.9 ± 7.9 W and the MAPE was 8.6 ± 1.5% for P2.For P3, the corresponding values were 79.2 ± 30.5 W and 10.6 ± 3.3%, respectively.
Results for the individualized parameters are shown in Table 2. Expressed as percentage of MR ae,max , A ve was 11.0 ± 1.1% and A acc was 6.9 ± 1.3%.In Table 3 the modeled overall metabolic energy contributions up to exhaustion from the three energy supply systems are shown for P2 and P3, respectively.Additionally, the maximum value of x 1 during P3 was 1.12 ± 0.03 (by definition 1.0 ± 0.0 during P2) meaning 12% more than the estimated maximum alactic capacity was used during P3.Also, the maximum value of x 4 during P3 was 1.21 ± 0.10 (0.99 ± 0.00 during P2) meaning on average 21% higher [mLa − ] than the estimated maximum from P2 was achieved during P3.
The differences between measured and modeled MR ae in relation to mean MR ae are shown with a Bland-Altman plot for P2 in Fig. 4A and for P3 in Fig. 4B.The mean difference shows a small underprediction by the model in P2 and an increased underprediction in P3.The confidence intervals show there are greater overall differences in P3 compared to P2.There is no obvious negative or positive trend in the data across the range of mean MR ae .
Visually, Fig. 5A shows overall excellent agreement in the time resolved MR ae /MR ae,max averaged across the subjects (n = 11) from measurement and model data, respectively.Still, significant differences between the model and measurements occur at several different stages during P2 (Fig. 5B). Figure 6A shows the corresponding results for P3, also with good overall agreement.The number of clusters with a t-value above the critical t-value are of the same magnitude as in P2 but less prominent (Fig. 6B), which is likely explained by the larger SD of measurements in P3 (Fig. 6A).
Figures 7A and 8A show the modeled energy systems contributions for one representative subject during P2 and P3, respectively.The contributions from different terms to MR dem are shown in Figs.7B and 8B during P2 and P3 respectively.There is an obvious difference between MR dem during the initial passive rest periods and the passive rest following exhaustion caused by the differences in MR ve and MR acc .

Model validity and reliability
The proposed bioenergetic model shows good overall agreement with the measurements of MR ae .The estimated values   6).Data for the period leading to volitional exhaustion has been stretched/squeezed to allow comparison across the subjects.The dashed lines show the critical t-value, indicating the limit for significant differences (p < 0.05) between measured and modeled data as do the estimated values of MR ve (Vella et al. 2006).The anaerobic energy contributions were on average 59.6 and 75.6 kJ during P2 and P3 respectively.The anaerobic energy contribution for one well-trained cyclist during a mass start cycling race was estimated to about 35 kJ (Skiba et al. 2012, Fig. 4).On the other hand, the mean anaerobic capacity is 93.2 kJ when estimated from the body masses for the subject in this study according to the findings of Andersson et al. (2022) regarding 3 min cycling time trials.Since the glycolysis is affected by the rate of accumulation and reduction of muscle metabolites (Stanley and Connett 1991), the exercise intensity and protocol will likely affect the lactic energy system contribution.Hence, it is problematic to define a general anaerobic capacity.The results in this study at least do not stand out as unreasonable.
The combined validity of the model formulation, parameter estimation, data acquisition, and data preprocessing (MAPE 8.6 ± 1.5% when fitted to P2 data) is not as good as the model proposed by Artiga Gonzalez et al. (2019), which, when fitted to an intermittent protocol, showed an MAPE of 3.4 ± 1.0%.However, applicability of the proposed model to a similar protocol, assessed as the difference in MAPE between P2 and P3, is seemingly better.MAPE was 10.6 ± 3.3% when the model fitted to P2 data was applied to P3; the model proposed by Artiga Gonzalez et al. (2019) showed MAPE of 10.8 ± 4.7% when the fitted model was applied to a race-like protocol.It should be noted that the protocols in this study are more similar, which is likely favorable for the proposed model.On the other hand, the smoothing applied to the data by Artiga Gonzalez et al. (2019) is likely favorable to their model in comparison to the smoothing applied in this study (this is further discussed below).
From the measured data in this study, it was visually obvious (not visible in this article) that three subjects exhibited larger systematic differences between P2 and P3.Quantifying the day-to-day variation, the RMSE between P2 and P3 of the last 60 s mean of the first 6 (submaximal) periods of each protocol was 46.5 ± 32.1 W for the whole group (n = 14), but 31.3 ± 7.9 W with these three subjects omitted (n = 11).The mean difference of the first 6 periods for these three subjects was 73.4,-70.9, and 101.7 W, respectively.With these subjects omitted, the overall RMSE of P3 was 65.2 ± 10.2 W (n = 11), compared to 79.2 ± 30.5 W for the whole group (n = 14).This may suggest greater potential for the proposed model as the RMSE of P3 and P2 are close (n = 11).However, mostly it emphasizes the limitations of predicting endurance performance with numerical modeling, with day-to-day variations having a greater influence than the relative accuracy of any proposed model.5).Data for the period leading to volitional exhaustion has been stretched/squeezed to allow comparison across the subjects.The dashed lines show the critical t-value, indicating the limit for significant differences (p < 0.05) between measured and modeled data

Performance quantification and comparison
Regarding RMSE and MAPE as metrics of model performance, it should be noted that the errors in both cases are calculated as the difference between the model predictions and measurements for each specific time sample.A minor misalignment in time between the results may give a comparatively large error, while the amplitude and shape of the curve could be perfectly matched.The implications of this will be greatest during transitions, when MR ae is rapidly changing.This issue will also be present in the Bland-Altman plots, as these too are applied to matched time samples.The equivalent mean squared error measure is used when fitting the model parameters, so these problems should be minimized.Nevertheless, this obscures the comparability of results and argues for an alternative quantification of curve agreement.One possible suggestion could be the mean of the closest distance from each respective data sample from the evaluated time series to the entire reference time series.
An additional issue with quantifying the model performance is the smoothing of measured data, which will have greatest effect in sharp transitions.In this study specifically, a rapid rise or fall of MR ae with a sudden change in demand will be smoothed to a somewhat more gradual change that starts earlier.Since the model is formulated to react to changes in demand, this will undoubtedly result in an error right before every major transition.A more gradual approach toward steady state is applied to the end of the transitions and, as such, these will not be as heavily affected.This phenomenon will directly affect the comparison of results between this study and the work of Artiga Gonzalez et al. (2019), since a greater degree of smoothing was applied to the dependent variable in the latter case.In addition, Artiga Gonzalez et al. (2019) smoothed the ergometer power output used as the model input (independent variable), which applied the above-mentioned smoothing effect in the transitions.Both these processes potentially reduce the error between the model and measurements, thus introducing some false validity.On the contrary, the power output in this study was smoothed in a fashion that retained the sharp transitions.

Physiological aspects
The proposed model provides possible descriptions of both the underlying mechanisms of the metabolic demand and the underlying control mechanisms of each of the energy supply systems.Unlike no previous model that we know of, the proposed model introduce feedback mechanisms from the metabolic energy supply system to the corresponding demand as  (Noordhof et al. 2015).Even so, the model must be considered a crude simplification of the real underlying mechanisms.Several other possible model variations and combinations were tested during preliminary parameter estimations.Reduced muscle efficiency was implemented as a factor proportional to [mLa − ] and multiplied by MR f .Even though such an effect has been reported (Hopker et al. 2017), this resulted in no improvement in RMSE but gave a somewhat random relative distribution between MR ve , MR acc , and the metabolic demand due to reduced muscle efficiency.This could be a sign of the system of equations being underdetermined in relation to input data.
The model proposed in this study uses Watts as the standard unit for metabolic rate, and, as such, measured oxygen uptake is recalculated into aerobic metabolic rate before model fitting or comparison.As a first note, the recalculation takes RQ into account, but this is not treated separately by the model, which could lead to small errors, especially at lower exertions and/or at the beginning of an exercise.No ramifications of this problem have been tested.Another related issue concerns including part of the glycolysis in the aerobic system i.e., the part of the glycolysis that produces the right amount of pyruvate to fuel mitochondrial respiration without increasing the concentration of accumulated lactate.In the model, this is included as a constant fraction of the aerobic system, but the level of glycolysis likely varies.Specifically, glycolysis is likely reduced in favor of mitochondrial respiration using the accumulated lactate as fuel during recovery following intense exercise (Gaesser and Brooks 1984).Aerobic energy output per oxygen input without glycolysis should be less than 10% lower than with full glycolysis, calculated using adenosine triphosphate output from the involved processes.This effect would be highest during passive or active rest e.g., when MR ae is less than 300 W resulting in an error less than 30 W. Preliminary testing was performed with a model compensating for this effect using an additional linearly dependent cost on MR dem below the maximum lactate steady state and with a maximum amplitude of 10% of MR ae .However, as this effect is limited, it was removed from the proposed model for the purpose of simplicity.
A two-compartment lactate accumulation model was tested initially, where lactate accumulated in the muscles could transfer into another compartment defined as blood and other tissue.Reduction of lactate was also possible in this second compartment to reflect the fact that the cardiopulmonary system favors lactate as fuel when available  (Gaesser and Brooks 1984).However, this introduced an additional ordinary differential equation that needed to be solved numerically as well as several additional parameters, again leading to an underdetermined system in relation to input data in which several significantly divergent parameter value combinations could yield similarly small values of RMSE.In addition, the division of glycolysis between the lactic and aerobic systems, in combination with a twocompartment lactate system, either fails to establish a steady state lactate concentration or calls for a pure curve fitting controlling mechanism rather than a physiologically plausible one.The two-compartment lactate accumulation model made it possible to include an additional metabolic rate related to the transfer of lactate.This aspect is not included in the proposed model, which has only one compartment for lactate accumulation and does not distinguish between lactate transfer and lactate reduction through oxidation in the muscles.
As P2 and P3 are intermittent protocols of varying power output, it is uncertain if and how the slow aerobic component will be visible.The most characteristic deviation from the typical exponential behavior of MR ae with a time constant around 25 s, can be seen in the example data for P2 in Fig. 7A.At the end of period 9, when the power output level drops, and during the following increase up to exhaustion, MR ae more closely resembles a linear decrease and increase.It is also visible in the mean curves in Fig. 5A, as the same trend was seen for several subjects.We interpret this as being related to the slow component.It is also clear that the proposed model does not capture this behavior well.This behavior is also seen in P3, starting from the end of the first 1-min period and at least until the start of the last period leading up to exhaustion.The underprediction of MR ae during this part of P3 leads to an overprediction of MR la , which partly explains why both x 1 and x 4 ([mLa − ]) increases above 1 in P3 for all subjects.Especially during the 1-min periods of lower intensity more pronounced recovery than modeled may occur.
The behavior of the slow component has previously been captured by applying a second exponential function with another, significantly larger, time constant and a time delay (Özyener et al. 2001).However, the concept of a time delay only really works with protocols characterized by a single step response.With an intermittent protocol, several questions arise.For instance, at what exercise intensity should the time delay countdown start, and should the time delay restart after a low intensity period?Nevertheless, there are still underlying control mechanisms not captured with sufficient accuracy by the proposed model.Instead of introducing a second exponential, testing was conducted in which the time constant of the primary exponential τ ae was increased exponentially (with varying degrees of power) as [mLa − ] increased.This unfortunately rendered the model too unstable to provide reasonable solutions.However, exhaustion might be the result of metabolic demand starting to rise in an unstable manner.Thus, the idea of an exertion-dependent time constant might be worth investigating further.

Conclusions
In summary, the present study proposes a new bioenergetic model which attempts to describe the dynamic behavior of the metabolic alactic, lactic and aerobic systems as well as contributions to the metabolic demand rates relating to both the fundamental ergometer cycling work, ventilation, and accumulation of muscle metabolites.A method of adapting the proposed model to reflect the bioenergetics of a specific athlete through grey-box parameter estimation is shown to be satisfactory with 14 estimated parameters.The validation of the model shows excellent overall agreement between measured and modeled aerobic metabolic rate during intermittent cycling by well-trained male cyclist.However, SD is 38.5% of the mean RMSE for P3, indicating not as good reliability.In part this is explained by day-to-day variation, but the model fails in the capturing the full details of the slow component and recovery periods, which likely makes it sensitive to protocol differences between parameter estimation and application.This argues for further development and validation to establish the full potential of the model.
need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http:// creat iveco mmons.org/ licen ses/ by/4.0/.
Fig.1Intermittent test protocols used during test day 2 (P2) and test day 3 (P3) where the numbers reflect the duration in minutes of each period of constant power.The constant power levels are shown to the left where PR is passive rest (no pedaling), AR is active rest, and L1,

Fig. 2
Fig. 2 Sampled and smoothed ergometer power output for one representative subject for P2.The data was smoothed by setting the power output equal to the mean power output for each section.Additionally, the 100 first seconds and final 60 s were deleted to avoid smoothing error and reduce calculation times

Fig. 3
Fig. 3 Sampled and smoothed breath-by-breath oxygen uptake (Panel A) and expiratory minute ventilation as STPD (Panel B) for one representative subject for P2.Data was smoothed using locally weighted linear regression with a 40-s time window and resampled to 0.5 Hz using linear interpolation

Fig. 4 Fig. 5
Fig. 4 Bland-Altman plot showing the mean of and difference between measured and modeled MR ae for all data samples from all subjects during P2 (Panel A) and P3 (Panel B) including the model-

Fig. 6
Fig. 6 Time resolved mean and SD (Panel A) and statistical parametric mapping from the two-tailed paired t-test (Panel B) of measured and modeled individual MR ae normalized with MR ae,max for P3 (n = 11).Results from three subjects are omitted due to deviations in duration during one or more periods (same subjects as in

Fig. 7
Fig. 7 Metabolic supply rates (Panel A) and metabolic demand rates (Panel B) for P2 in a representative subject (S8).MR al is alactic metabolic rate, MR la lactic metabolic rate, MR ae modeled aerobic metabolic rate, MR dem total metabolic demand rate and Measured MR ae is the smoothed measured aerobic metabolic rate.MR acc is meta-

Fig. 8
Fig. 8 Metabolic supply rates (Panel A) and metabolic demand rates (Panel B) for P3 in a representative subject (S8).MR al is alactic metabolic rate, MR la lactic metabolic rate, MR ae modeled aerobic metabolic rate, MR dem total metabolic demand rate and Measured MR ae is the smoothed measured aerobic metabolic rate.MR acc is meta-

Table 1
Ergometer power output levels used in the intermittent protocols (P2 and P3)

Table 2
Bioenergetic model parameters and constants with optimized individual values from the parameter estimation as mean ± SD (n = 14), initial values used in the parameter estimation, allowed ranges, and a brief description All combinations of the initial values were used in the parameter estimation, totaling 64 initial value combinations.Initial values reported as mean ± SD indicate that individual initial values based on P1a and/or P2 were used.K was set to 0.8 for all subjects based on preliminary parameter estimations a Individual values taken from P2 measurements b Individual values taken from P1a measurements c Percentage of initial parameter value d Percentage of MR ae,max initial value e Percentage of V m initial value

Table 3
Modeled metabolic energy contributions up to exhaustion from the two intermittent test protocols (P2 and P3) respectively, reported as amount of energy (kJ) and as percentage of total metabolic energy output The first 100 s of passive rest in each protocol was omitted (see Data Collection) Gonzalez et al. 2019), and MR rest are similar to their estimated values from measurements.The estimated values of τ ae agree with earlier findings (ArtigaGonzalez et al. 2019),