# The Role of a Heavy Neutrino in the Gamma-Ray Burst GRB221009A

Kingman Cheunga,b a Department of Physics and Center for Theory and
Computation, National Tsing Hua University, Hsinchu 30013, Taiwan
b Division of Quantum Phases and Devices, School of Physics, Konkuk
University, Seoul 143-701, Republic of Korea

###### Abstract

Recently, several telescopes, including Swift-BAT, GBM, and LHAASO, have
observed the ever highest-energy and long-duration gamma-rays from a gamma-ray
burst named as GRB221009A (located at a red-shift of $z=0.151$) on October 9,
2022. Conventional understanding tells us that very high-energy photons
produced at such a far distance suffer severe attenuation before reaching the
Earth. We propose the existence of a sub-MeV to $O(10)$ MeV heavy neutrino
with a transitional magnetic dipole moment, via which the heavy neutrino is
produced at the GRB. It then travels a long distance to our galaxy and decays
into a neutrino and a photon, which is observed. In such a way, the original
high-energy photon produced at the GRB can survive long-distance attenuation.

Recently, the ever highest-energy and long-duration gamma rays were detected
on Oct. 9, 2022, first by the Swift Burst Alert Telescope (BAT) BAT  and the
Fermi Gamma-ray Burst Monitor (GBM) FERMI . Subsequently, observations of the
very highest gamma rays were recorded from the same source by the FERMI-LAT
LAT-1 ; LAT-2 , LHAASO LHAASO , and the Carpet-2 Carpet-2 . Such a gamma-ray
burst (GRB) was measured at a red-shift of $z=0.151$ redshift ; redshift-2 ,
which corresponds to about 720 Mpc ($\approx 2\times 10^{25}\,{\rm m}$). The
energy of the very high energy (VHE) photons is determined to be at least 10
TeV and perhaps as high as 250 TeV.

The detection of very high-energy photons at such a far distance cannot be
explained by conventional physics. The current understanding of such VHE
photons is due to the explosion of a supermassive star and the creation of a
black hole. In such a violent environment, it is not difficult to imagine very
high-energy collisions taking place and thus creating a lot of hadrons such as
pions, kaons, etc. The energetic neutral pions then decay into photons and
charged pions into muons and neutrinos. However, since the GRB is at a very
far distance from us, the VHE photons will lose most of the energies along
their path to us by pair creation, Compton scattering, and other mechanisms.
Thus, the detection of such VHE photons by satellite experiments or
experiments on the Earth is beyond our understanding.

An explanation of photon-axion-photon conversion was put forward to explain
the anomaly Galanti:2022pbg ; Baktash:2022gnf ; Lin:2022ocj ; Nakagawa:2022wwm
; Zhang:2022zbm . The VHE photons so-produced are converted into axions by the
interaction $f_{a}aF_{\mu\nu}\tilde{F}^{\mu\nu}$. The axions then travel a
long distance to near our galaxy without interacting with the intergalactic
space and are converted back to photons in presence of a magnetic field of our
galaxy. In this way, the original energy of the VHE photons is preserved.
Other possible interpretations include Lorentz inverse violation Finke:2022swf
; Li:2022wxc ; Zhu:2022usw , modification of ultrahigh-energy cosmic ray
spectrum Das:2022gon ; AlvesBatista:2022kpg .

The idea of photon-axion-photon conversion is based on the fact that axions
rarely interact along the path of propagation. Another well-known particle
that shares this property is the neutrino. How does a VHE photon related to a
neutrino? We propose the existence of a heavy neutrino (denoted by $N$) of a
mass around $O(10^{-2})-O(10)$ MeV, which has a transitional magnetic dipole
moment with an active neutrino. In the violent environment around the GRB, it
is not difficult to imagine that there are numerous high-energy hadronic
collisions, which produce a large number of pions, kaons, etc. Thus, we
suggest that the neutral pion can decay into a photon, an active neutrino, and
a heavy neutrino, and also the charged pion into a muon, a photon, and a heavy
neutrino

$\pi^{0}\to\gamma\gamma^{*}\to\gamma\nu
N\;,\qquad\pi^{\pm}\to\mu^{\pm}\nu^{*}\to\mu^{\pm}\gamma N\;,$ (1)

in which the flavor of the active neutrino is not important. The branching
ratio of this decay is small, but it is not critical in our interpretation,
because we do not know exactly the number of collisions that can happen there.
The more important is that the lifetime of the heavy neutrino is long enough
such that it can survive the path coming towards our galaxy. We give more
detail in the following. The key ingredient of this interpretation is that
once the heavy neutrino $N$ is produced near or at the GRB, it acquires an
energy as high as tens or hundreds of TeV. It then travels a cosmological
distance without decay or attenuation, followed by its decay $N\to\nu\gamma$
when it comes close to our galaxy. In this way, we can observe VHE photons of
energies $O(10)-O(100)$ TeV.

The transitional magnetic dipole moment of the heavy neutrino (HN) is
parameterized as

${\cal
L}=\frac{1}{2}\mu_{\nu}\overline{\nu_{L}}\sigma^{\mu\nu}N\,F_{\mu\nu}\;+\;{\rm
h.c.}\;,$ (2)

where $F_{\mu\nu}$ is the field strength of the photon field and $\mu_{\nu}$
is the magnetic dipole moment of the active-to-heavy-neutrino transition. This
interaction is responsible for both the production and decay of the HN.
Production of $N$ can proceed in meson decays, such as
$\pi^{\pm}\to\mu^{\pm}\nu^{*}\to\mu^{\pm}(\gamma N)$ and
$\pi^{0}\to\gamma\gamma^{*}\to\gamma\nu N$, or Primakoff upscattering
$\nu\gamma^{*}\to N$ Magill:2018jla .

Since there are large uncertainties in the hadronic environment and energy
profile of particles at the GRB, we do not have precise knowledge about the
number of high-energy pions that can be produced there. Nevertheless, we know
the total amount of energy given off in a typical supernova explosion is about
$10^{44}\,{\rm J}\approx 6.2\times 10^{62}\,{\rm eV}$ 111 The energy of the
GRB221009A was estimated to be about $2\times 10^{54}\,{\rm erg}\equiv
1.2\times 10^{66}\,{\rm eV}$ redshift . and 99% of this energy is carried away
by neutrinos. We perform a naive estimate that only 0.1% of the total energy
is going into high-energy hadronic collisions and only 10% of this collision
energy goes to pions. The number of high-energy pions ($~{}\sim 100$ TeV) that
can be produced is as many as $(10^{58}\,{\rm GeV})/(10^{14}\,{\rm GeV})\sim
O(10^{44})$. Even though the required magnetic dipole moment $\mu_{\nu}$ is
smaller than $10^{-10}\,{\rm GeV}^{-1}$ (shown later in Fig. 1) due to the
SN1987A constraint Magill:2018jla , the number of $N$ produced should be large
enough to account for the observation.

On the other hand, the decay of $N\to\nu\gamma$ is more important in our
discussion. We calculate the decay width of $N$ based on the Lagrangian in Eq.
(2). In the rest frame of $N$,

$\Gamma(N\to\gamma\nu)=\frac{|\mu_{\nu}|^{2}}{8\pi}M_{N}^{3}\;,$ (3)

where we have ignored the mass of the active neutrino, which is tiny compared
to $M_{N}$. The decay length is then given by

$L_{\rm decay}=\beta\gamma c\tau\;,$ (4)

where $\beta\simeq 1$, $\gamma\simeq E_{N}/M_{N}$ and $\tau=1/\Gamma$.
Substituting into the above equation we obtain

$L_{\rm decay}=8\pi cE_{N}\frac{1}{M_{N}^{4}}\frac{1}{|\mu_{\nu}|^{2}}\;.$ (5)

To obtain a representative value of $L_{\rm decay}$ we take $E_{N}=100$ TeV,
$M_{N}=10^{-1}$ MeV, and $\mu_{\nu}=10^{-9}\,{\rm GeV}^{-1}$, we obtain

$L_{\rm decay}\simeq(5\times 10^{24}\,{\rm m})\,\left(\frac{E_{N}}{100\,{\rm
TeV}}\right)\,\left(\frac{10^{-1}\,{\rm
MeV}}{M_{N}}\right)^{4}\,\left(\frac{10^{-9}\,{\rm
GeV}^{-1}}{\mu_{\nu}}\right)^{2}\;.$ (6)

With this choice of $E_{N}$, $M_{N}$ and $\mu_{\nu}$, the decay length of $N$
falls in the right ballpark of the distance required to reach from the GRB to
our galaxy. Other choices of $M_{N}$ and $\mu_{\nu}$ can be obtained by the
scaling in Eq. (6).

Figure 1:  Contours of constant decay length $L_{\rm
decay}=10^{24},10^{25},10^{26}$ m as a function of $(M_{N},\mu_{\nu})$. The
constraint from the supernova SN1987A is extracted from Ref. Magill:2018jla .

We show a contour plot of decay length $L_{\rm decay}=10^{24},10^{25},10^{26}$
m for the heavy neutrino with $E_{N}=100$ TeV as a function of $M_{N}$ and
$\mu_{\nu}$ in Fig. 1. There exist a number of constraints on the transitional
magnetic moment of a heavy neutrino for mass between $10^{-2}$ and 100 MeV
Brdar:2020quo . Among the constraints, the relevant one to our study is from
supernova SN1987A Magill:2018jla . From the figure, when the magnetic dipole
moment is below $\sim 10^{-10}\,{\rm GeV}^{-1}$ the decay length is in the
correct ballpark of the required distance for $M_{N}\simeq 0.3-10$ MeV.

To conclude VHE photons produced at the GRB221009A, if traveled directly to
us, would suffer severe attenuation. We have proposed the existence of a sub-
MeV to $O(10)$ MeV heavy neutrino with a transitional magnetic dipole moment,
via which the heavy neutrino is produced in the neutral and/or charged pion
decays. It then travels a long distance to our galaxy and decays into a
neutrino and a photon, which is observed. In such a way, the original high-
energy photons produced at the GRB can survive long-distance attenuation.

###### Acknowledgements.

We thank D. Alexander Kann for communication on the distance of the GRB. This
work was supported by MoST with grant nos. MoST-110-2112-M-007-017-MY3.

## References

  * (1) S. Dichiara, J. D. Gropp, J. A. Kennea, N. P. M. Kuin, A. Y. Lien, F. E. Marshall, A. Tohuvavohu, and M. A. Williams. GCN Circ. 32632 (2022).
  * (2) P. Veres, E. Burns, E. Bissaldi, S. Lesage, and . Robert. GCN Circ. 32636 (2022).
  * (3) E. Bissaldi, N. Omodei, and M. Kerr. GCN Circ. 32637 (2022).
  * (4) R. Pillera, E. Bissaldi, N. Omodei, and F. L. G. La Mura. GCN Circ. 32658 (2022).
  * (5) Y. Huang, S. Hu, S. Chen, M. Zha, C. Liu, Z. Yao, and Z. Cao. GCN Circ. 32677 (2022).
  * (6) D. Dzhappuev and et al. The Astronomer’s Telegram 15669, 1 (2022).
  * (7) A. de Ugarte Postigo. GCN Circ. 32648 (2022).
  * (8) D. Malesani. GCN Circ. 32765 (2022).
  * (9) G. Galanti, M. Roncadelli, and F. Tavecchio, Explanation of the very-high-energy emission from GRB221009A, arXiv:2210.05659.
  * (10) A. Baktash, D. Horns, and M. Meyer, Interpretation of multi-TeV photons from GRB221009A, arXiv:2210.07172.
  * (11) W. Lin and T. T. Yanagida, Electroweak axion in light of GRB221009A, arXiv:2210.08841.
  * (12) S. Nakagawa, F. Takahashi, M. Yamada, and W. Yin, Axion dark matter from first-order phase transition, and very high energy photons from GRB 221009A, arXiv:2210.10022.
  * (13) G. Zhang and B.-Q. Ma, Axion-photon conversion of LHAASO multi-TeV and PeV photons, arXiv:2210.13120.
  * (14) J. D. Finke and S. Razzaque, Possible Evidence for Lorentz Invariance Violation in Gamma-ray Burst 221009A, arXiv:2210.11261.
  * (15) H. Li and B.-Q. Ma, Lorentz invariance violation induced threshold anomaly versus very-high energy cosmic photon emission from GRB 221009A, arXiv:2210.06338.
  * (16) J. Zhu and B.-Q. Ma, Light speed variation from GRB 221009A, arXiv:2210.11376.
  * (17) S. Das and S. Razzaque, Ultrahigh-energy cosmic-ray signature in GRB 221009A, arXiv:2210.13349.
  * (18) R. Alves Batista, GRB 221009A: a potential source of ultra-high-energy cosmic rays, arXiv:2210.12855.
  * (19) G. Magill, R. Plestid, M. Pospelov, and Y.-D. Tsai, Dipole Portal to Heavy Neutral Leptons, Phys. Rev. D 98 (2018), no. 11 115015, [arXiv:1803.03262].
  * (20) V. Brdar, A. Greljo, J. Kopp, and T. Opferkuch, The Neutrino Magnetic Moment Portal: Cosmology, Astrophysics, and Direct Detection, JCAP 01 (2021) 039, [arXiv:2007.15563].


∎

11institutetext:  Fabian C. Weigend 22institutetext: School of Computer, Data
and Mathematical Sciences and School of Health Sciences
Western Sydney University
22email<EMAIL_ADDRESS>33institutetext: David C. Clarke
44institutetext: Department of Biomedical Physiology and Kinesiology and the
Sports Analytics Group
Simon Fraser University 55institutetext: Oliver Obst 66institutetext: School
of Computer, Data and Mathematical Sciences
Western Sydney University 77institutetext: Jason Siegler 88institutetext:
College of Health Solutions
Arizona State University

# A hydraulic model outperforms work-balance models for predicting recovery
kinetics from intermittent exercise

Fabian C. Weigend  David C. Clarke  Oliver Obst  Jason Siegler

(Received: date / Accepted: date)

###### Abstract

Data Science advances in sports commonly involve “big data”, i.e., large
sport-related data sets. However, such big data sets are not always available,
necessitating specialized models that apply to relatively few observations.
One important area of sport-science research that features small data sets is
the study of recovery from exercise. In this area, models are typically fitted
to data collected from exhaustive exercise test protocols, which athletes can
perform only a few times. Recent findings highlight that established recovery
such as the so-called work-balance models are too simple to adequately fit
observed trends in the data. Therefore, we investigated a hydraulic model that
requires the same few data points as work-balance models to be applied, but
promises to predict recovery dynamics more accurately.

To compare the hydraulic model to established work-balance models, we
retrospectively applied them to data compiled from published studies. In
total, one hydraulic model and three work-balance models were compared on data
extracted from five studies. The hydraulic model outperformed established
work-balance models on all defined metrics, even those that penalize models
featuring higher numbers of parameters. These results incentivize further
investigation of the hydraulic model as a new alternative to established
performance models of energy recovery.

###### Keywords:

Human Athletic Performance Mathematical Models Cycling Critical Power W’
Recovery Bioenergetics High-Intensity Interval Training

††journal: Currently Under Review

##  1 Introduction

Emerging technologies that enable the real-time monitoring of athletes in
training and competition have fostered interest in methods to predict and
optimize athlete performance. Predictive models for how much an athlete “has
left in the tank” enable the investigation of pacing strategies (Behncke,,
1997; Sundström et al.,, 2014; de Jong et al.,, 2017) and to dynamically
adjust strategies to optimize the outcome of a competition (Hoogkamer et al.,,
2018). They can be described as a digital athlete, i.e., a computer-based
model for enhancing training programming or strategy optimization. A
foundation for such advances is research in performance modeling, which can be
understood as the mathematical abstraction of exercise physiology.

###  1.1 Critical-power-based approaches to performance modeling

One of the seminal models in the area of performance modeling is the critical
power model, which relies on the notions of a critical power ($CP$) and a
finite energy reserve for work above critical power ($W^{\prime}$) (Hill,,
1993; Whipp et al.,, 1982). Monod and Scherrer, (1965) defined $CP$ as: “the
maximum work rate a muscle can keep up for a very long time without fatigue”.
Thus, $CP$ can be considered as a threshold for sustainable exercise.
$W^{\prime}$ represents a capacity for work to be performed at a rate above
$CP$ and is conceptualized as an energy storage. Using the definitions of time
to exhaustion (TTE) and constant power output ($P$) the critical power model
can be summarized in the relationship

$\lx@glossaries@gls@link{acronym}{tte}{\leavevmode\text{TTE}}=\frac{\lx@glossaries@gls@link{acronym}{w^{\prime}}{\leavevmode
W^{\prime}}}{\lx@glossaries@gls@link{acronym}{p}{\leavevmode
P}-\lx@glossaries@gls@link{acronym}{cp}{\leavevmode CP}}.$ (1)

To determine $CP$ and $W^{\prime}$, an athlete has to conduct between three
and five exercise tests until exhaustion at various constant exercise
intensities. $CP$ and $W^{\prime}$ are then fitted to these distinct TTE and
$P$ observations and their relationship can be used to predict the time to
exhaustion at other intensities. Hill, (1993) emphasized that the
attractiveness of the model lies in its coarse simplicity and it should not be
employed if highly accurate predictions are required. Nevertheless, its
straightforward application and its elegant abstraction have led to improved
understanding and prediction of performance dynamics (Poole et al.,, 2016;
Sreedhara et al.,, 2019; Vanhatalo et al.,, 2011).

While the critical power model predicts energy expenditure at high
intensities, it does not consider the recovery of $W^{\prime}$ after exercise
has ended or during exercise at low intensities. Formally, exercise protocols
that alternate between intensities below $CP$ and above $CP$ constitutes as
intermittent exercise. In order to predict performance capabilities of
athletes during intermittent exercise, models need to predict recovery of
$W^{\prime}$ during phases of exercise below the $CP$ intensity.

One of the most widely covered approaches to predict recovery of $W^{\prime}$
during intermittent exercise is the $W^{\prime}_{\mathrm{bal}}$ model
(Sreedhara et al.,, 2019; Jones and Vanhatalo,, 2017). Since the first
publication by Skiba et al., (2012), an updated form of
$W^{\prime}_{\mathrm{bal}}$ was introduced by Skiba et al., (2015) and another
alternative form was proposed by Bartram et al., (2018).
$W^{\prime}_{\mathrm{bal}}$ models have been used to search for optimal
drafting strategies in running (Hoogkamer et al.,, 2018) or to predict phases
of perceived exhaustion during cycling exercise (Skiba et al.,, 2014).

Despite these advances, research into energy recovery modeling is an evolving
field, and $W^{\prime}_{\mathrm{bal}}$ models have been scrutinized for their
limitations. Similar to energy expenditure dynamics, recovery dynamics are
derived from exhaustive exercise tests and therefore data for model fitting
and validation are sparse (Vanhatalo et al.,, 2011; Sreedhara et al.,, 2019).
Recent findings suggest that current $W^{\prime}_{\mathrm{bal}}$ models overly
simplify energy recovery dynamics, and that model modifications that account
for characteristics of prior exhaustive exercise (Caen et al.,, 2019) as well
as bi-exponential recovery dynamics (Caen et al.,, 2021) might improve
recovery predictions. These modified models feature additional parameters,
which introduces challenges in fitting them to small data sets. Indeed, the
search for models that optimally balance complexity with applicability to few
data points is a primary challenge of energy recovery modeling.

###  1.2 Hydraulic models of human performance

Hydraulic models offer an alternative to $W^{\prime}_{\mathrm{bal}}$ models
for predicting energy expenditure and recovery dynamics during exercise.
Instead of using $CP$ and $W^{\prime}$, they represent energy dynamics as
liquid flow within a system of tanks and pipes. These tanks and pipes are
arranged according to physiological parameters such as maximal oxygen uptake
and estimated phosphocreatine levels of an athlete. The first hydraulic model
was proposed by Margaria, (1976) to provide an intuitive conceptualization of
bioenergetic responses to exercise. Morton, (1986) further elaborated the
model, formalized its dynamics with differential equations, and published it
as the Margaria-Morton (M-M) model. Later, Sundström, (2016) proposed an
extension of the M-M model and named it the Margaria-Morton-Sundström (M-M-S)
model. Compared to $W^{\prime}_{\mathrm{bal}}$ models, hydraulic models can
predict more complex energy expenditure and recovery dynamics and have the
potential to address highlighted shortcomings of $W^{\prime}_{\mathrm{bal}}$
models in recent literature.

A challenge of these hydraulic models is that their parameters require in-
depth knowledge about bioenergetic systems. Indeed, Morton, (2006) concluded
that it remained to be seen to what extent model predictions conform to
reality. Also, the more recently proposed M-M-S model by Sundström, (2016) has
yet to be validated experimentally. Behncke, (1997) applied the M-M model to
world records in competitive running, and while the predictions agreed with
values provided in the literature, he also pointed out situations in which the
naive interpretation of the model would not be justified. Furthermore,
Behncke, (1997) stated that constraints dictated by physiological conditions
made explicit computations with the M-M model “rather cumbersome”.
Collectively, the requirement to set parameters according to physiological
measures impede the application and validation of the M-M and M-M-S models.

To overcome the issues caused by ascribing the model parameters to concrete
bioenergetic measures, we proposed a generalized form of the M-M hydraulic
model in Weigend et al., (2021). Our generalized hydraulic model allows the
fitting of its parameters using an optimization approach that only requires
$CP$ and $W^{\prime}$ as inputs. In this way, our modified model preserves the
flexibility needed to model the observed dynamics without requiring strict
correspondence to parameters pertaining bioenergetics. In a proof-of-concept,
we showed that our fitted hydraulic model could successfully predict both
energy expenditure and recovery kinetics for one example case, the former in
line with predictions of the $CP$ model and the latter in a manner that
matches published observations. While the generalized hydraulic model
demonstrated satisfactory predictivity, it is still unknown as to whether it
can outperform the existing $W^{\prime}_{\mathrm{bal}}$ models.

Therefore, in this work, we compare the prediction quality of our generalized
hydraulic model from Weigend et al., (2021) to that of three
$W^{\prime}_{\mathrm{bal}}$ models. We hypothesized the hydraulic model would
predict the observed recovery ratios compiled from past studies overall more
accurately than the $W^{\prime}_{\mathrm{bal}}$ models. We found that the
hydraulic model outperformed the $W^{\prime}_{\mathrm{bal}}$ models on
objective goodness-of-fit and prediction metrics. We conclude that the
generalized hydraulic model provides a beneficial new perspective on energy
recovery modeling that should be investigated further.

##  2 Material and methods

The _Materials and methods_ are structured in the following way. In Section
2.1, the $W^{\prime}_{\mathrm{bal}}$ and hydraulic models are defined and
their underlying assumptions specified. We then introduce a new
$W^{\prime}_{\mathrm{bal}}$ model that has been fitted to the same recovery
data as the investigated hydraulic model. In Section 2.2, we discuss how we
will objectively compare the $W^{\prime}_{\mathrm{bal}}$ and hydraulic models.
In particular, we propose a procedure to obtain comparable recovery ratios.
The validation data set consists of data compiled from previously published
studies on the recovery from exercise. Section 2.3 lists the data exclusion
criteria and extraction procedures. Finally, in Section 2.4, we describe the
metrics used to assess the model goodness-of-fits and prediction capabilities.

###  2.1 Model definitions

The $CP$, $W^{\prime}_{\mathrm{bal}}$, and the hydraulic models feature
assumptions and parameters that require defining.

####  2.1.1 Energy expenditure and recovery with $W^{\prime}_{\mathrm{bal}}$

The critical power model predicts energy expenditure and is the underlying
model for the $W^{\prime}_{\mathrm{bal}}$ models. The four essential
assumptions of the critical power model are stated as follows (Hill,, 1993;
Morton,, 2006):

  * 1.

An individual’s power output is a function of two energy sources: aerobic
(using oxidative metabolism) and anaerobic (non-oxidative metabolism).

  * 2.

Aerobic energy is unlimited in capacity but its conversion rate into power
output is limited ($CP$).

  * 3.

Anaerobic energy is limited in capacity ($W^{\prime}$) but its conversion rate
is unlimited.

  * 4.

Exhaustion occurs when all of $W^{\prime}$ is depleted.

These assumptions are reflected in Equation 1, in which time to exhaustion
(TTE) is estimated from the available $W^{\prime}$ divided by work above $CP$.
At every time step during which an athlete exercises above $CP$, the product
of the time elapsed and the difference between $CP$ and the actual power
output is subtracted from the energy capacity $W^{\prime}$. Thus, during a
constant power output above $CP$, the critical power model predicts a linear
depletion of $W^{\prime}$. When $W^{\prime}$ is depleted, exhaustion is
reached.

Figure 1: Example energy expenditure and recovery predictions of two models
for intermittent exercise. $W^{\prime}$ expenditure during exercise above $CP$
(the red dashed line in the lower panel) was modeled using the $CP$ model
while $W^{\prime}$ recovery during exercise below $CP$ was modeled using
either the $W^{\prime}_{\mathrm{skib}}$ or $W^{\prime}_{\mathrm{bart}}$.
Differences in predicted recovery kinetics are clearly visible.

As observable in the example in Figure 1, subsequently established
$W^{\prime}_{\mathrm{bal}}$ models combine the assumed linear depletion of
$W^{\prime}$ at power outputs above $CP$ with predictions for $W^{\prime}$
reconstitution during exercise below $CP$. The initial
$W^{\prime}_{\mathrm{bal}}$ model by Skiba et al., (2012) was later updated in
Skiba et al., (2015). Substantial differences between these versions exist,
and as shown by Skiba and Clarke, (2021), the original model by Skiba et al.,
(2012) contradicts the assumption of Equation 1 that $W^{\prime}$ linearly
depletes. As such, we focused on the updated version by Skiba et al., (2015)
and we refer to it henceforth as $W^{\prime}_{\mathrm{bal-ode}}$. We denote
the remaining capacity of $W^{\prime}$ at a discrete time point $t$ during
exercise as $W^{\prime}_{\mathrm{bal-ode}_{t}}$. $P_{t}$ refers to the power
output at a discrete time step $t$. $\Delta t$ is the difference between the
discrete time step $t-1$ and $t$ in seconds. We define the overall
$W^{\prime}_{\mathrm{bal-ode}}$ model as

$W^{\prime}_{\mathrm{bal-ode}_{t}}=\begin{cases}W^{\prime}_{\mathrm{bal-
ode}_{t-1}}-(P_{t}-\lx@glossaries@gls@link{acronym}{cp}{\leavevmode CP})\Delta
t,&\text{for}\;P_{t}\geq\lx@glossaries@gls@link{acronym}{cp}{\leavevmode
CP}\\\ \lx@glossaries@gls@link{acronym}{w^{\prime}}{\leavevmode
W^{\prime}}-(\lx@glossaries@gls@link{acronym}{w^{\prime}}{\leavevmode
W^{\prime}}-W^{\prime}_{\mathrm{bal-ode}_{t-1}})\cdot e^{\frac{-\Delta
t}{\mathcal{T}_{t}}},&\text{for}\;P_{t}<\lx@glossaries@gls@link{acronym}{cp}{\leavevmode
CP}.\\\ \end{cases}$ (2)

During a constant power output above or at $CP$
($P_{t}\geq\lx@glossaries@gls@link{acronym}{cp}{\leavevmode CP}$),
$W^{\prime}_{\mathrm{bal-ode}_{t}}$ decreases linearly as $t$ increases, like
the critical power model predicts. During power outputs below $CP$,
$W^{\prime}_{\mathrm{bal-ode}_{t}}$ increases exponentially with $W^{\prime}$
as its asymptote. $\mathcal{T}_{t}$ affects recovery speed and varies between
distinct $W^{\prime}_{\mathrm{bal-ode}}$ models. At a discrete time step $t$,
the $\mathcal{T}_{t}$ of the $W^{\prime}_{\mathrm{bal-ode}}$ model by Skiba et
al., (2015) ($\mathcal{T}_{\mathrm{skib}_{t}}$) is estimated as

$\mathcal{T}_{\mathrm{skib}_{t}}=\frac{\lx@glossaries@gls@link{acronym}{w^{\prime}}{\leavevmode
W^{\prime}}}{D_{\lx@glossaries@gls@link{acronym}{cp}{\leavevmode CP}_{t}}},$
(3)

where $D_{\lx@glossaries@gls@link{acronym}{cp}{\leavevmode CP}_{t}}$
represents the difference between $P_{t}$ and $CP$. Henceforth, we refer to
Equation 2 with Equation 3 as $W^{\prime}_{\mathrm{skib}}$. Figure 1 depicts
an example for $W^{\prime}_{\mathrm{skib}}$ predictions. For the time steps of
the first 3 minutes $P_{t}$ was below $CP$ and $W^{\prime}_{\mathrm{bal-
ode}_{t}}$, i.e., available $W^{\prime}$ balance, remained at its maximum.
Then, the power output increased above $CP$ for the next three minutes. The
available $W^{\prime}$ balance decreased by
$P_{t}-\lx@glossaries@gls@link{acronym}{cp}{\leavevmode CP}$ per second.
Between 6 and 9 minutes $P_{t}$ dropped below $CP$ again,
$W^{\prime}_{\mathrm{skib}}$ simulated recovery, and $W^{\prime}_{\mathrm{bal-
ode}_{t}}$ rose exponentially with $W^{\prime}$ as its asymptote. Speed of
recovery was affected by $\mathcal{T}_{\mathrm{skib}_{t}}$ from Equation 3,
which took the difference between $P_{t}$ and $CP$ into account. That is
observable by comparing recovery between 6 and 9 minutes to recovery between
12 and 21 minutes. During the second recovery bout $P_{t}$ was lower and
therefore the slope of the exponential recovery was steeper. During the last
three minutes $P_{t}$ was equal to $CP$ and thus $W^{\prime}_{\mathrm{bal-
ode}_{t}}$ did not change. If $W^{\prime}_{\mathrm{bal-ode}_{t}}$ would reach
0, exhaustion would be predicted. In the example in Figure 1 the athlete was
predicted to be close to exhaustion, but some of their energy capacities
remained.

Bartram et al., (2018) investigated the recovery rate of $W^{\prime}$ of elite
cyclists and observed faster recovery rates than Skiba et al., (2015).
Therefore, they proposed another $\mathcal{T}_{t}$ for Equation 2 to predict
quicker recovery ratios. The $\mathcal{T}_{t}$ by Bartram et al., (2018)
($\mathcal{T}_{\mathrm{bart}_{t}}$) was defined as

$\mathcal{T}_{\mathrm{bart}_{t}}=2287.2\cdot{D_{\lx@glossaries@gls@link{acronym}{cp}{\leavevmode
CP}_{t}}}^{-0.688}.$ (4)

Henceforth, we will refer to Equation 2 with Equation 4 as
$W^{\prime}_{\mathrm{bart}}$. Predictions of $W^{\prime}_{\mathrm{bart}}$ are
depicted alongside those of $W^{\prime}_{\mathrm{skib}}$ in the example in
Figure 1. It is observable that $W^{\prime}_{\mathrm{bart}}$ predicted faster
recovery dynamics.

####  2.1.2 The hydraulic model

In Weigend et al., (2021) we expressed our generalized hydraulic tank model
mathematically as a system of discretized differential equations with 8
parameters (A.3 and A.4 in their Appendix). Henceforth, the model will be
referred to as the $\mathrm{hydraulic}_{\mathrm{weig}}$, a schematic of which
is depicted in Figure 2. $\mathrm{Hydraulic}_{\mathrm{weig}}{}{}$ models power
output of an athlete as a function of three interacting energy sources, which
are represented as liquid-containing tanks. As depicted in Figure 2, these
tanks are named the aerobic energy source ($Ae$), anaerobic fast energy source
($AnF$) and anaerobic slow energy source ($AnS$). $Ae$ is assumed to have
infinite volume, which is indicated by the fading color to the left. A pipe
connects $Ae$ to the middle tank and has the maximal flow capacity $m^{Ae}$.
The pipe from the right tank into the middle tank allows flow in both
directions and the maximal flow capacity $m^{AnS}$ from $AnS$ into $AnF$ and
$m^{AnF}$ from $AnF$ into $AnS$. A tap ($p$) is attached to the bottom of the
middle tank $AnF$ and liquid flow from this tap represents energy demand. The
fill levels of tanks and flows through pipes change as liquid flows from the
tap. If the athlete expends energy, the liquid level of $AnF$ drops and
initiates flow from $Ae$ and $AnS$ into $AnF$. When the middle tank is empty,
liquid flow out of $p$ no longer matches the demand and exhaustion is assumed.

In the depicted situation in Figure 2 $p$ was opened and liquid flowed out of
$AnF$. As a result, the liquid level in $AnF$ dropped and thus liquid started
to flow from $Ae$ into $AnF$. The more the fill level of the middle tank $AnF$
dropped, the less liquid pressured against the pipe exit from $Ae$, and the
more the flow from $Ae$ increased. In Figure 2 the flow out of $p$ was so
large, that the fill level of $AnF$ dropped below the top of $AnS$
($h>\theta$) and liquid from $AnS$ started to flow into $AnF$ too. The liquid
volume in $AnF$ is limited such that it can only contribute to flow out of $p$
for a limited time. If the simulated athlete stopped exercise, then their
power output would decrease to 0 and the tap $p$ would close. If the tap was
closed in the depicted situation in Figure 2, liquid from the outer tanks $Ae$
and $AnS$ would refill the middle tank $AnF$ until its fill level rises above
the fill level of $AnS$. Then, liquid from $Ae$ would continue to refill $AnF$
and liquid from $AnF$ would flow into $AnS$. The model would mimic recovery
and $Ae$ would eventually refill both other tanks.

Figure 2: A three-component hydraulic model as defined by Weigend et al.,
(2021). Tanks are conceptualized as the aerobic energy source ($Ae$),
anaerobic fast energy source ($AnF$), and anaerobic slow energy source
($AnS$). $Ae$ is assumed to be infinite in volume, which is indicated by the
fading color to the left. A tap $p$ is attached to the bottom of $AnF$ and
flow from it represents energy demand. Pipes connecting the three tanks have
maximal flow capacities $m^{Ae}$, $m^{AnS}$, and $m^{AnF}$.

Morton, (1986) mathematically expressed liquid flows within this system as
first- and second-order ordinary differential equations and in Weigend et al.,
(2021) we extended these equations so that they apply to all possible
configurations of their generalized model. Due to liquid pressure and flow
dynamics, varying fill levels of the middle and the right tank affect how flow
from the tap $p$ is estimated. With these interactions between three tanks,
$\mathrm{hydraulic}_{\mathrm{weig}}$ is capable of predicting energy
expenditure and recovery as a more complex function than the above introduced
critical power and $W^{\prime}_{\mathrm{bal-ode}}$ models. It features eight
adjustable parameters. Depicted in Figure 2, the parameters $\phi$, $\gamma$,
$\theta$ represent tank positions, $AnF$, $AnS$ tank capacities, and $m^{Ae}$,
$m^{AnF}$, $m^{AnS}$ maximal flow capacities. A configuration $c$ for
$\mathrm{hydraulic}_{\mathrm{weig}}$ entails the positions, sizes and
capacities of each tank and is therefore defined as

$c=[AnF,AnS,m^{Ae},m^{AnS},m^{AnF},\phi,\theta,\gamma].$ (5)

Fitting the $\mathrm{hydraulic}_{\mathrm{weig}}$ to an athlete means finding
the configuration that enables the model output to best reproduce the observed
exercise responses of an athlete. In Weigend et al., (2021) we introduced an
evolutionary computation workflow to derive such configurations. We fitted a
configuration to $CP$ and $W^{\prime}$ measures of an athlete as well as
recovery ratios derived from a publication by Caen et al., (2019). The same
recovery ratios are used for every fitting and thus, to fit
$\mathrm{hydraulic}_{\mathrm{weig}}$ to an athlete, only $CP$ and $W^{\prime}$
of the athlete are required. These are the same measures that are needed to
apply $W^{\prime}_{\mathrm{bal-ode}}$ models and, hence, the required input
measures are the same for all compared $\mathrm{hydraulic}_{\mathrm{weig}}$
and $W^{\prime}_{\mathrm{bal-ode}}$ models. With specific reference to our
fitting method in Weigend et al., (2021), 10 evolutionary fittings to given
$CP$ and $W^{\prime}$ were estimated and the best fitting has been used.

####  2.1.3 An additional $W^{\prime}_{\mathrm{bal}}$ model

The above introduced models $W^{\prime}_{\mathrm{skib}}$,
$W^{\prime}_{\mathrm{bart}}$, and $\mathrm{hydraulic}_{\mathrm{weig}}$ were
each created by fitting to different sets of recovery observations. Objective
metrics to compare model quality, e.g. the Akaike Information Criterion
(Burnham and Anderson,, 2004), require models to be fitted to the same data.
Therefore, to allow a more comprehensive comparison, we added a third
$W^{\prime}_{\mathrm{bal-ode}}$ model with a new
$\mathcal{T}_{\mathrm{weig}_{t}}$. We derived this
$\mathcal{T}_{\mathrm{weig}_{t}}$ with a procedure as close as possible to the
ones of Skiba et al., (2012) and Bartram et al., (2018).

As the first step, a constant value for $\mathcal{T}_{t}$ in Equation 2 was
fitted to each recovery ratio and recovery time combination from Table 1 of
the Appendix of Weigend et al., (2021). For these observations power output
was constant for every discrete time step $t$ during recovery and thus
$\mathcal{T}_{t}$ could be considered as constant with the same value for all
$t$. We used the standard Broyden-Fletcher-Goldfarb-Shanno algorithm
implementation of SciPy (SciPy 1.0 Contributors et al.,, 2020) with 200 as the
initial guess to fit a constant $\mathcal{T}_{t}$ that enabled Equation 2 to
best reproduce the observed recovery ratio. This resulted in twelve pairs of
fitted constant $\mathcal{T}_{t}$s to constant recovery intensities.

As the next step, we then fitted an exponential function to these twelve pairs
using the non-linear least squares implementation of SciPy (SciPy 1.0
Contributors et al.,, 2020). With the recovery intensity as
$D_{\lx@glossaries@gls@link{acronym}{cp}{\leavevmode CP}_{t}}$, the function
was of the form $\mathcal{T}_{t}=a\cdot
e^{D_{\lx@glossaries@gls@link{acronym}{cp}{\leavevmode CP}_{t}}\cdot{b}}+c$.
With the values of Skiba et al., (2012) as the initial guess (546, -0.01,
316), the resulting optimal constants were $a=1274.45,b=-0.0308,$ and
$c=266.65$. Thus, given any
$D_{\lx@glossaries@gls@link{acronym}{cp}{\leavevmode CP}_{t}}$ at a discrete
time step $t$, $\mathcal{T}_{\mathrm{weig}_{t}}$ can be estimated as

$\mathcal{T}_{\mathrm{weig}_{t}}=1274.45\cdot e^{-0.0308\cdot
D_{{\lx@glossaries@gls@link{acronym}{cp}{\leavevmode CP}}_{t}}}+266.65.$ (6)

Unfortunately, this fitted equation failed to satisfactorily fit the data
($R^{2}=0.14$) but we nevertheless used it because it was developed using a
procedure that closely resembled those used to estimate
$\mathcal{T}_{\mathrm{skib}_{t}}$ and $\mathcal{T}_{\mathrm{bart}_{t}}$. The
introduction of $\mathcal{T}_{\mathrm{weig}_{t}}$ is valuable because it
allows the application of the Akaike Information Criterion metric, which
requires compared models to be fitted to the same data points. Henceforth,
Equation 2 with Equation 6 will be referred to as $W^{\prime}_{\mathrm{bal-
ode}}$ model with $\mathcal{T}_{\mathrm{weig}_{t}}$
($W^{\prime}_{\mathrm{weig}}$).

###  2.2 Procedure for computing comparable recovery predictions

WB1 intensity: $P_{\mathrm{work}}$ time: until exhaustion RB intensity:
$P_{\mathrm{rec}}$ time: $T_{\mathrm{rec}}$ WB2 intensity: $P_{\mathrm{work}}$
time: until exhaustion Figure 3: A schematic of the protocol to estimate
recovery ratios. An exhaustive work bout (WB1) at a set intensity
($P_{\mathrm{work}}$) is prescribed. Immediately after exhaustion is reached,
a recovery bout (RB) follows at a lower recovery intensity
($P_{\mathrm{rec}}$) for a set duration ($T_{\mathrm{rec}}$). Then, a second
exhaustive work bout (WB2) is conducted and the ratio of the time to
exhaustion of WB2 to the one of WB1 represents the amount that was recovered
during RB.

We compare the abilities of above defined $W^{\prime}_{\mathrm{bal-ode}}$ and
$\mathrm{hydraulic}_{\mathrm{weig}}$ models to predict “recovery ratios”.
Recovery ratios are computed from exercise protocols involving two exhaustive
work bouts (WB1 and WB2) interspersed with a recovery bout (RB). A schematic
of the protocol is depicted in Figure 3. First, the model simulates exercise
at a fixed work intensity ($P_{\mathrm{work}}$) above $CP$ until exhaustion.
Immediately after exhaustion is reached, exercise intensity switches to a
lower recovery intensity ($P_{\mathrm{rec}}$) below $CP$. After a set time
($T_{\mathrm{rec}}$) at that recovery intensity, a second work bout (WB2)
until exhaustion at $P_{\mathrm{work}}$ is simulated. The time to exhaustion
of WB2 is expected to be shorter than the one of WB1 due to the limited
recovery time in between the work bouts. Because it is assumed that
$W^{\prime}$ is completely depleted at the end of WB1, the ratio of the second
time to exhaustion to the first represents the amount of $W^{\prime}$
recovered. Thus, the time to exhaustion of WB2 divided by the time to
exhaustion of WB1 multiplied by 100 results in a recovery ratio in percent
(%).

This outlined procedure aligns with the assumptions of $CP$ or
$W^{\prime}_{\mathrm{bal-ode}}$ models and enables the direct comparison of
the simulated recovery ratios from each model
$\mathrm{hydraulic}_{\mathrm{weig}}$, $W^{\prime}_{\mathrm{skib}}$,
$W^{\prime}_{\mathrm{bart}}$, $W^{\prime}_{\mathrm{weig}}$, and with the
published data. As an example, the recovery ratio curves in Figure 5 were
obtained by estimating the WB1 $\rightarrow$ RB $\rightarrow$ WB2 protocol for
every recovery duration ($T_{\mathrm{rec}}$) in seconds between 0 and 900
seconds. Published observations by Caen et al., (2021) were added to the plot.

###  2.3 Data extraction for model comparisons

We extracted data from previous studies that investigated energy recovery
dynamics and used to it compare and evaluate recovery ratio predictions of all
models. The studies for comparison were identified from Table 1 of the
comprehensive review by Chorley and Lamb, (2020). From these studies, we
retained those that featured appropriate data, except those that met the
following exclusion criteria:

  * •

Featured a mode of exercise other than cycling. Cycle ergometers measure power
output directly. Power outputs during modes of exercise like running or
swimming are not directly comparable because they are estimated using
different methods or are approximated, e.g, (Morton and Billat,, 2004) focused
only on speed instead of power.

  * •

The observations were made under extreme conditions, e.g., hypoxia or
altitude.

  * •

Insufficient information was reported to simulate the prescribed protocol in
and/or to infer a recovery ratio of $W^{\prime}$ in percent, e.g., the
integral version of the $W^{\prime}_{\mathrm{bal}}$ model by Skiba et al.,
(2012) assumes recovery during high-intensity exercise such that recovery
ratios cannot be straightforwardly inferred.

  * •

The prescribed protocol leaves doubt if reported recovery ratios are
comparable to the “recovery estimation protocol” described earlier, e.g.,
repeated ramp tests until exhaustion, 50% $W^{\prime}$ depletion followed by a
3-min all-out test, or knee-extension maximal voluntary contraction (MVC) test
during recovery.

Five studies were included for comparison, four of which were obtained from
the Chorley and Lamb review: (Bartram et al.,, 2018), (Chidnok et al.,, 2012),
(Ferguson et al.,, 2010), and (Caen et al.,, 2019). After the summary of
Chorley and Lamb, (2020) was published, Caen et al., (2021) published a study
that investigated the $W^{\prime}$ reconstitution dynamics in even more detail
and which was thus added to the list.

The data in the listed studies were presented in diverse ways, such that
modifications were made to some of the data to enable model comparison. The
study by Caen et al., (2019) did not report distinct mean values for every
investigated condition, such that we derived approximate values in Weigend et
al., (2021) to fit our $\mathrm{hydraulic}_{\mathrm{weig}}$ to their
conditions. Hence, the data for comparison are the values from Weigend et al.,
(2021). Further, the study by Bartram et al., (2018) fitted their own
$W^{\prime}_{\mathrm{bart}}$ model, where $\mathcal{T}_{t}$ is defined
according to Equation 4. Therefore, we used $W^{\prime}_{\mathrm{bart}}$ model
predictions for prescribed intensities of Bartram et al., (2018) as the
observations against which the other models were compared.

The study by Chidnok et al., (2012) reported times to exhaustion from their
intermittent exercise protocol instead of recovery ratios. Power output during
recovery was constant in their tests. Therefore, in order to derive recovery
ratio estimations that are comparable with the WB1 $\rightarrow$ RB
$\rightarrow$ WB2 procedure defined above, we fitted a constant value for
$\mathcal{T}_{t}$ of the $W^{\prime}_{\mathrm{bal-ode}}$ model to each of
their prescribed protocols and times to exhaustion. These constant values for
$\mathcal{T}_{t}$ were fitted with the Brent method implementation by SciPy
(SciPy 1.0 Contributors et al.,, 2020) to find a local minimum in the interval
between [100, 1000]. We then used WB1 $\rightarrow$ RB $\rightarrow$ WB2
recovery ratio estimations of $W^{\prime}_{\mathrm{bal-ode}}$ models with
fitted constant $\mathcal{T}_{t}$ as the observations with which to compare
$W^{\prime}_{\mathrm{skib}}$, $W^{\prime}_{\mathrm{bart}}$,
$W^{\prime}_{\mathrm{weig}}$, and $\mathrm{hydraulic}_{\mathrm{weig}}$.

###  2.4 The metrics of goodness of fit

The metrics of goodness of fit used to compare the models were Root Mean
Squared Error (RMSE), Mean Absolute Error (MAE), and the small-sample version
of the Akaike Information Criterion ($\mathrm{AIC}_{\mathrm{c}}$). Chai and
Draxler, (2014) discussed RMSE and MAE as widely adopted metrics for assessing
model prediction capabilities. We compared predictive accuracy by comparing
RMSE and MAE on data to which competing models were not fitted. Lower values
for RMSE and MAE were interpreted as more accurate predictions.

To statistically compare prediction error distributions between models, we
used a bootstrap hypothesis test (Efron and Tibshirani,, 1993; Good,, 2000).
We did so because only small data sets were available and we could not assume
normal distributed prediction errors with equal variances for every compared
model. The null hypothesis of our bootstrap test was that prediction error
distributions of two compared models are the same. Because we used two
prediction error metrics (RMSE and MAE) we investigated the null hypothesis on
both. We used the absolute difference between RMSE and also between MAE of
compared groups as our test statistics. With the null hypothesis that error
distributions are the same, we could bootstrap new samples by randomly
selecting with replacement from all pooled observations. We created a
distribution of test statistics from 1 000 000 bootstrap samples to reliably
approximate the p-value of our observed test statistic at high precision. We
rejected the null hypothesis if the p-value$<.05$.

We also compared models with the $\mathrm{AIC}_{\mathrm{c}}$, which was first
proposed by Sugiura, (1978). The $\mathrm{AIC}_{\mathrm{c}}$ is a model
selection tool to investigate the balance between model complexity and
explanatory capability (Burnham and Anderson,, 2004).
$\mathrm{AIC}_{\mathrm{c}}$ penalizes the number of parameters of the model
and thus provides insight into the balance between model complexity and
goodness of fit. The lower the $\mathrm{AIC}_{\mathrm{c}}$ score, the better
this balance is met. The $\mathrm{AIC}_{\mathrm{c}}$ was calculated as

$\text{AIC}_{c}=n\cdot\text{ln}(\text{MSE})+2k+\frac{2k\cdot(k+1)}{n-k-1},$
(7)

where MSE is the Mean Squared Error, $n$ is the number of data points and $k$
is the number of parameters of the model. Models have to be fitted to and
applied to the same data in order to obtain comparable
$\mathrm{AIC}_{\mathrm{c}}$ scores. Therefore, only
$W^{\prime}_{\mathrm{weig}}$ and $\mathrm{hydraulic}_{\mathrm{weig}}$ were
comparable with this criterion in this work.

Altogether, the hypothesis that the more complex
$\mathrm{hydraulic}_{\mathrm{weig}}$ model fits the data better than the
established $W^{\prime}_{\mathrm{bal-ode}}$ models will be supported if the
overall RMSE, MAE and $\mathrm{AIC}_{\mathrm{c}}$ scores are lower for
$\mathrm{hydraulic}_{\mathrm{weig}}$ than for other models, and if prediction
error distributions are significantly different to those of other models.

##  3 Results

In the following section, we present the extracted data and the prediction
results of $W^{\prime}_{\mathrm{bal-ode}}$ and
$\mathrm{hydraulic}_{\mathrm{weig}}$ models for each listed previous study. We
refer to extracted data from studies by the last name of the first author,
e.g., the extracted data from Bartram et al., (2018) is referred to as
“Bartram data set”. All studies collected their data through performance tests
that required athletes to exercise until volitional exhaustion. Such tests are
affected by circumstances that are hard to measure and control, e.g,
motivation, nutrition, and state-of-mind. Therefore, recovery ratio
observations are noisy and the extracted group averages were accompanied by
large standard deviations. These uncertainties prevented us from drawing
conclusions about model quality on averages of individual data sets and
instead necessitated to perform the comparison between models across all
available data. We begin by presenting the extracted data and model
predictions of individual data sets throughout Sections 3.1, 3.2, 3.3, 3.4 and
3.5 followed by summarizing all prediction errors and resulting RMSE, MAE, and
$\mathrm{AIC}_{\mathrm{c}}$ scores in the final Section 3.6.

###  3.1 Bartram data set

Figure 4: Comparison of model predictions with the defined recovery estimation protocol (WB1 $\rightarrow$ RB $\rightarrow$ WB2). Depicted are the recovery dynamics around 60 seconds at various $D_{CP}$ recovery intensities after a preceding exhaustive exercise at the intensity that is predicted to lead to exhaustion after 100 seconds ($P100$). Chosen intensities and time frames stem from the protocol prescribed by Bartram et al., (2018) and predictions of $W^{\prime}_{\mathrm{bart}}$ were used as the observations with which to compare models. Table 1: The left part of the table summarizes extracted data and conditions from Bartram et al., (2018). The right part of the table displays model predictions. Predictions of $W^{\prime}_{\mathrm{bart}}$ were taken as the observed recovery ratios. Parameters from Bartram et al., (2018) | Observed | Predicted recovery ratio
---|---|---
$CP$ | $W^{\prime}$ | $P_{\mathrm{work}}$ | $P_{\mathrm{rec}}$ | $T_{\mathrm{rec}}$ | recovery ratio | $W^{\prime}_{\mathrm{skib}}$ | $W^{\prime}_{\mathrm{weig}}$ | $\mathrm{hydraulic}_{\mathrm{weig}}$
$393$ | $23\,300$ | $626$ | $393$ | $60$ | $0.0$ | $0.0$ | $0.0$ | $22.7$
$393$ | $23\,300$ | $626$ | $343$ | $60$ | $33.0$ | $12.1$ | $10.6$ | $33.9$
$393$ | $23\,300$ | $626$ | $293$ | $60$ | $47.0$ | $22.8$ | $16.9$ | $44.8$
$393$ | $23\,300$ | $626$ | $243$ | $60$ | $57.0$ | $32.1$ | $19.4$ | $52.7$
$393$ | $23\,300$ | $626$ | $193$ | $60$ | $64.0$ | $40.3$ | $20.0$ | $59.3$

The protocol prescribed by Bartram et al., (2018) consisted of three work
bouts interspersed with two 60-second recovery bouts. The first two work bouts
each lasted for 30 seconds, and the final one until volitional exhaustion.
Work bout exercise intensity ($P_{\mathrm{work}}$) was set to $P100$, i.e.,
the intensity that was predicted to lead to exhaustion after 100 seconds. The
recovery bout intensity ($P_{\mathrm{rec}}$) was set to differences to $CP$
($D_{CP}$) of 200, i.e., $CP$\- 200 watts, or 150, 100, 50, or 0. The group
averaged CP and W’ for the four world-class cyclists featured in Bartram et
al., (2018) were 393 watts and 23,300 joules. Altogether, these input values
resulted in an estimated $P100$ exhaustive intensity of 626 watts and recovery
intensities $D_{CP}$ 0 of 393 watts, $D_{CP}$ 50 of 343 watts, $D_{CP}$ 100 of
293 watts, $D_{CP}$ 150 of 243 watts, and $D_{CP}$ 200 of 193 watts,
respectively.

The resulting recovery predictions of $W^{\prime}_{\mathrm{bal-ode}}$ and
$\mathrm{hydraulic}_{\mathrm{weig}}$ models are summarized in Figure 4 and
Table 1. The $W^{\prime}_{\mathrm{bart}}$ model was not compared because it
was the model that Bartram et al., (2018) fitted to their observations and we
used it to create the observations against which the other models were
compared. The fitted $\mathrm{hydraulic}_{\mathrm{weig}}$ configuration to
$CP$ and $W^{\prime}$ by Bartram et al., (2018) was: $[$23111.91, 65845.28,
391.57, 148.88, 24.15, 0.73, 0.01, 0.24$]$.

Figure 4 and Table 1 show that in all cases except $D_{CP}$0 the recovery
ratios predicted by $\mathrm{hydraulic}_{\mathrm{weig}}$ model were closest to
the ones observed by Bartram et al., (2018), followed by
$W^{\prime}_{\mathrm{skib}}$ and then $W^{\prime}_{\mathrm{weig}}$. On the
contrary, the hydraulic model is the only model to predict recovery at
$D_{CP}$0.

###  3.2 Caen data set

Figure 5: Comparison of model predictions with published observations by Caen et al., (2021). After an exhaustive exercise bout at $P240$, recovery dynamics at an intensity of 161 watts were simulated using the defined recovery estimation protocol of this work (WB1 $\rightarrow$ RB $\rightarrow$ WB2). Published observed recovery ratios by Caen et al., (2021) are depicted in blue. Table 2: The left part of the table summarizes extracted data and conditions from Caen et al., (2021). The right part of the table displays model predictions. Parameters from Caen et al., (2021) | Observed | Predicted recovery ratio
---|---|---
$CP$ | $W^{\prime}$ | $P_{\mathrm{work}}$ | $P_{\mathrm{rec}}$ | $T_{\mathrm{rec}}$ | recovery ratio | $W^{\prime}_{\mathrm{bart}}$ | $W^{\prime}_{\mathrm{skib}}$ | $W^{\prime}_{\mathrm{weig}}$ | $\mathrm{hydraulic}_{\mathrm{weig}}$
$269$ | $19\,200$ | $349$ | $161$ | $30$ | $28.6$ | $28.0$ | $15.5$ | $9.2$ | $26.9$
$269$ | $19\,200$ | $349$ | $161$ | $60$ | $34.8$ | $48.2$ | $28.7$ | $17.5$ | $41.2$
$269$ | $19\,200$ | $349$ | $161$ | $120$ | $44.2$ | $73.2$ | $49.1$ | $31.9$ | $49.8$
$269$ | $19\,200$ | $349$ | $161$ | $180$ | $50.5$ | $86.1$ | $63.7$ | $43.8$ | $52.8$
$269$ | $19\,200$ | $349$ | $161$ | $240$ | $55.1$ | $92.8$ | $74.1$ | $53.6$ | $54.7$
$269$ | $19\,200$ | $349$ | $161$ | $300$ | $56.8$ | $96.3$ | $81.5$ | $61.8$ | $56.3$
$269$ | $19\,200$ | $349$ | $161$ | $600$ | $73.7$ | $99.9$ | $96.6$ | $85.4$ | $64.9$
$269$ | $19\,200$ | $349$ | $161$ | $900$ | $71.3$ | $100.0$ | $99.4$ | $94.4$ | $73.8$

The protocol by Caen et al., (2021) investigated the recovery dynamics
following exhaustive exercise at $P240$ (published average of 349 watts). They
prescribed a recovery intensity of 161 watts on average, which was determined
by selecting $90\%$ of the power at gas exchange threshold (Binder et al.,,
2008) of their participants. The average $CP$ of their participants was 269
watts, and the average $W^{\prime}$ 19 200 joules. The reported observed
recovery ratios were $28.6\pm 8.2\%$ after 30 seconds, $34.8\pm 11.1\%$ after
60 seconds, $44.2\pm 9.7\%$ after 120 seconds, $50.5\pm 12.1\%$ after 180
seconds, $55.1\pm 13.3\%$ after 240 seconds, $56.8\pm 16.4\%$ after 300
seconds, $73.7\pm 19.3\%$ after 600 seconds, and $71.3\pm 20.8\%$ after 900
seconds.

The simulation parameters and results of the defined recovery estimation
protocol are summarized in Figure 5 and Table 2. Fitting
$\mathrm{hydraulic}_{\mathrm{weig}}$ to $CP$ and $W^{\prime}$ group averages
resulted in the configuration $[$17631.06, 46246.13, 267.28, 117.50, 20.09,
0.68, 0.01, 0.29$]$. The recovery ratios predicted by the
$\mathrm{hydraulic}_{\mathrm{weig}}$ model better matched the observed values
compared to all the other models. Nevertheless, some lack of fit for the
$\mathrm{hydraulic}_{\mathrm{weig}}$ model was observed: the model
overpredicted the recovery ratios at early time points and underpredicted
those at longer time points except for the last one.
$W^{\prime}_{\mathrm{skib}}$, $W^{\prime}_{\mathrm{bart}}$, and
$W^{\prime}_{\mathrm{weig}}$ model predictions consistently overestimated
recovery for longer recovery times.

###  3.3 Chidnok data set

Chidnok et al., (2012) prescribed a protocol that alternated between 60-second
work bouts and 30-second recovery bouts until the athlete reached exhaustion.
With their protocol, the work-bout intensity $P_{\mathrm{work}}$ was set to
$P240$. The protocol prescribed four trials each with a different recovery
intensity $P_{\mathrm{rec}}$ (20 watts as the “low” recovery intensity, 95
watts as “medium”, 173 watts as “high”, and 270 watts as the “severe” recovery
intensity). The participants had an average $CP$ of 241 watts and $W^{\prime}$
of 21 100 joules. Their recorded times to exhaustion were $1224\pm 497$
seconds with “low” recovery intensity, $759\pm 243$ seconds with “medium”,
$557\pm 90$ seconds with “high”, and $329\pm 29$ seconds with “severe”.

Figure 6: Predicted recovery dynamics of compared models up to 60 seconds after a preceding exhaustive exercise at $P240$ and at three different recovery intensities (20 watts, 95 watts, and 173 watts). Observations were predicted recovery ratios of $W^{\prime}_{\mathrm{bal-ode}}$ models with a constant $\mathcal{T}$ fitted to reported times to exhaustion by Chidnok et al., (2012). Table 3: The left part of the table summarizes extracted data and conditions from Chidnok et al., (2012). The right part of the table displays model predictions. Parameters from Chidnok et al., (2012) | Observed | Predicted recovery ratio
---|---|---
$CP$ | $W^{\prime}$ | $P_{\mathrm{work}}$ | $P_{\mathrm{rec}}$ | $T_{\mathrm{rec}}$ | recovery ratio | $W^{\prime}_{\mathrm{bart}}$ | $W^{\prime}_{\mathrm{skib}}$ | $W^{\prime}_{\mathrm{weig}}$ | $\mathrm{hydraulic}_{\mathrm{weig}}$
$241$ | $21\,100$ | $329$ | $20$ | $30$ | $16.6$ | $41.6$ | $27.0$ | $10.6$ | $40.6$
$241$ | $21\,100$ | $329$ | $95$ | $30$ | $21.4$ | $33.3$ | $18.8$ | $10.1$ | $30.9$
$241$ | $21\,100$ | $329$ | $173$ | $30$ | $24.4$ | $21.3$ | $9.2$ | $6.8$ | $20.5$

As described in Section 2.3, in order to compare observations of Chidnok et
al., (2012) to WB1 $\rightarrow$ RB $\rightarrow$ WB2 protocol estimations, a
constant value for $\mathcal{T}_{t}$ for the $W^{\prime}_{\mathrm{bal-ode}}$
model was fitted to the protocol by Chidnok et al., (2012) for each of their
recovery conditions. The resulting $\mathcal{T}_{t}$ values were 165.19
seconds for the “low” recovery intensity protocol, a $\mathcal{T}_{t}$ of
124.81 seconds for “medium”, and a $\mathcal{T}_{t}$ of 107.45 seconds for
“high”. The “severe” recovery intensity was left out because 270 watts lies
above the average $CP$ of 241 watts. In this case, no recovery should occur if
the assumptions of $W^{\prime}_{\mathrm{bal-ode}}$ model hold true. Chidnok et
al., (2012) prescribed recovery bouts of 30 seconds and WB1 $\rightarrow$ RB
$\rightarrow$ WB2 protocol estimations with corresponding fitted
$\mathcal{T}_{t}$s and with a $T_{\mathrm{rec}}$ of 30 seconds were 24.6% at
the “low” intensity, 21.7% at “medium”, and 16.7% at the “high” recovery
intensity.

The fitted $\mathrm{hydraulic}_{\mathrm{weig}}$ configuration to $CP$ and
$W^{\prime}$ by Chidnok et al., (2012) was $[$18919.76, 48051.77, 239.55,
115.05, 19.48, 0.68, 0.05, 0.31$]$. Predictions of all models and extracted
conditions for the recovery estimation protocol are summarized in Figure 6 and
Table 3. In the case of the “low” recovery intensity predictions of the
$W^{\prime}_{\mathrm{skib}}$ and $W^{\prime}_{\mathrm{weig}}$ models were the
most accurate. In the case of the “medium” recovery intensity the
$W^{\prime}_{\mathrm{skib}}$ model was the most accurate, and in the remaining
“high” condition $\mathrm{hydraulic}_{\mathrm{weig}}$ and
$W^{\prime}_{\mathrm{bart}}$ model predictions were closest to the data. None
of the models made predictions that were close to all three of the
observations.

###  3.4 Ferguson data set

Figure 7: A comparison of predicted recovery dynamics after an exhaustive
exercise bout at $P360$ and at a recovery intensity of 20 watts. Recovery
ratios are estimated with the (WB1 $\rightarrow$ RB $\rightarrow$ WB2)
protocol, which resembles the prescribed protocol by Ferguson et al., (2010).
Published observations by Ferguson et al., (2010) are depicted in blue.

Ferguson et al., (2010) prescribed a protocol with an initial time to
exhaustion bout at the intensity that was predicted to lead to exhaustion
after 360 seconds ($P360$) followed by a recovery at 20 watts for 2 minutes, 6
minutes, or 15 minutes. After recovery, exercise intensity was then increased
back to one of three possible high-intensity work rates. Thus, each
participant performed nine tests in total with three different constant work
rates after three different recovery times. The $CP$ model was fitted to these
three times to exhaustion after each recovery period to determine changes in
$CP$ and $W^{\prime}$. Ferguson et al., (2010) published their group averages
for $CP$ as 212 watts, $W^{\prime}$ as 21 600 joules, the $P360$ as 269 watts,
and the observed recovery ratios after 2 minutes as ($37\%\pm 5\%$), 6 minutes
($65\%\pm 6\%$), and 15 minutes ($86\%\pm 4\%$).

Table 4: The left part of the table summarizes extracted data and conditions from Ferguson et al., (2010). The right part of the table displays model predictions. Parameters from Ferguson et al., (2010) | Observed | Predicted recovery ratio
---|---|---
$CP$ | $W^{\prime}$ | $P_{\mathrm{work}}$ | $P_{\mathrm{rec}}$ | $T_{\mathrm{rec}}$ | recovery ratio | $W^{\prime}_{\mathrm{bart}}$ | $W^{\prime}_{\mathrm{skib}}$ | $W^{\prime}_{\mathrm{weig}}$ | $\mathrm{hydraulic}_{\mathrm{weig}}$
$212$ | $21\,600$ | $269$ | $20$ | $120$ | $37.0$ | $85.8$ | $65.6$ | $35.9$ | $54.2$
$212$ | $21\,600$ | $269$ | $20$ | $360$ | $65.0$ | $99.7$ | $95.9$ | $73.6$ | $69.4$
$212$ | $21\,600$ | $269$ | $20$ | $900$ | $86.0$ | $100.0$ | $100.0$ | $96.4$ | $98.4$

Extracted parameters for the recovery intensity protocol and model prediction
results are summarized in Figure 7 and Table 3 together with reported means by
Ferguson et al., (2010). The fitted $\mathrm{hydraulic}_{\mathrm{weig}}$
configuration to $CP$ and $W^{\prime}$ group averages by Ferguson et al.,
(2010) was $[$18730.05, 81030.54, 211.56, 94.31, 18.76, 0.63, 0.21, 0.34$]$.
In this setup $W^{\prime}_{\mathrm{weig}}$ was overall closest to published
observations. $\mathrm{Hydraulic}_{\mathrm{weig}}{}{}$ overestimated the
recovery after 120 and after 900 seconds. $W^{\prime}_{\mathrm{skib}}$ and
$W^{\prime}_{\mathrm{bart}}$ overestimated recovery in every instance.

###  3.5 Weigend data set

We derived the values from Table 1 in the Appendix of our Weigend et al.,
(2021) publication from Caen et al., (2019). Reported measures recreate the
depicted means in Figure 3 of the publication by Caen et al., (2019). They
consisted of three recovery ratios for four conditions each: Preceding
exhausting exercise at $P240$ or $P480$ followed by recovery at 33% of $CP$ or
66% of $CP$. The participants of Caen et al., (2019) had an average $CP$ of
248 watts and $W^{\prime}$ of 18 200 joules, which results in a $P240$ of 285
watts, a $P480$ of 323 watts, 33% of $CP$ as 81 watts, and 66% of $CP$ as 163
watts.

Figure 8: Predicted recovery dynamics in comparison to measures that we derived from observations of Caen et al., (2019). We derived three recovery ratios for four conditions each: Preceding exhausting exercise at $P240$ or $P480$ followed by recovery at 33% of $CP$ or 66% of $CP$. Depicted observations are the values from Table 1 in the Appendix of our publication Weigend et al., (2021) and approximate Figure 3 of the publication by Caen et al., (2019). $W^{\prime}_{\mathrm{weig}}$ and $\mathrm{hydraulic}_{\mathrm{weig}}$ were fitted to these observations. Table 5: The left part of the table summarizes extracted data and conditions that Weigend et al., (2021) derived from Caen et al., (2019). The right part of the table displays model predictions. Both $W^{\prime}_{\mathrm{weig}}$ and $\mathrm{hydraulic}_{\mathrm{weig}}$ were fitted to these observations. Their predicted recovery ratios were recorded for the $\mathrm{AIC}_{\mathrm{c}}$ goodness of fit estimation metric in Section 2.4. Parameters from Weigend et al., (2021) | Observed | Predicted recovery ratio
---|---|---
$CP$ | $W^{\prime}$ | $P_{\mathrm{work}}$ | $P_{\mathrm{rec}}$ | $T_{\mathrm{rec}}$ | recovery ratio | $W^{\prime}_{\mathrm{bart}}$ | $W^{\prime}_{\mathrm{skib}}$ | $W^{\prime}_{\mathrm{weig}}$ | $\mathrm{hydraulic}_{\mathrm{weig}}$
$248$ | $18\,200$ | $323$ | $81$ | $120$ | $55.0$ | $83.1$ | $66.7$ | $35.5$ | $58.3$
$248$ | $18\,200$ | $323$ | $81$ | $240$ | $61.0$ | $97.1$ | $89.0$ | $58.3$ | $65.0$
$248$ | $18\,200$ | $323$ | $81$ | $360$ | $70.5$ | $99.5$ | $96.3$ | $73.1$ | $70.1$
$248$ | $18\,200$ | $323$ | $163$ | $120$ | $49.0$ | $67.2$ | $42.9$ | $28.4$ | $46.5$
$248$ | $18\,200$ | $323$ | $163$ | $240$ | $55.0$ | $89.2$ | $67.4$ | $48.7$ | $51.5$
$248$ | $18\,200$ | $323$ | $163$ | $360$ | $58.0$ | $96.5$ | $81.4$ | $63.3$ | $54.2$
$248$ | $18\,200$ | $285$ | $81$ | $120$ | $42.0$ | $83.0$ | $66.8$ | $35.5$ | $46.8$
$248$ | $18\,200$ | $285$ | $81$ | $240$ | $52.0$ | $97.1$ | $89.0$ | $58.3$ | $54.0$
$248$ | $18\,200$ | $285$ | $81$ | $360$ | $59.5$ | $99.5$ | $96.3$ | $73.1$ | $60.4$
$248$ | $18\,200$ | $285$ | $163$ | $120$ | $38.0$ | $67.2$ | $42.9$ | $28.4$ | $38.5$
$248$ | $18\,200$ | $285$ | $163$ | $240$ | $37.5$ | $89.3$ | $67.4$ | $48.7$ | $43.6$
$248$ | $18\,200$ | $285$ | $163$ | $360$ | $50.0$ | $96.5$ | $81.4$ | $63.3$ | $47.4$

Extracted parameters for the recovery ratio estimation protocol and model
predictions are summarized in Figure 8 and Table 5. The best fit
$\mathrm{hydraulic}_{\mathrm{weig}}$ configuration to $CP$ and $W^{\prime}$
group averages was $[$18042.06, 46718.18, 247.4, 106.77, 16.96, 0.72, 0.02,
0.25$]$. As described earlier, the recovery ratio values by Weigend et al.,
(2021) were used to fit the $\mathcal{T}_{\mathrm{weig}_{t}}$ for the
$W^{\prime}_{\mathrm{weig}}$ model and they are used in the evolutionary
fitting process for $\mathrm{hydraulic}_{\mathrm{weig}}$ to fit recovery
dynamics. Therefore, both $W^{\prime}_{\mathrm{weig}}$ and
$\mathrm{hydraulic}_{\mathrm{weig}}$ were not scrutinized for predictive
accuracy on this data set. Their predicted recovery ratios were recorded for
the $\mathrm{AIC}_{\mathrm{c}}$ goodness of fit estimation metric in covered
in the next subsection. Out of the remaining two models predictions of
$W^{\prime}_{\mathrm{skib}}$ were closer to the observations but both
overpredict in nearly all instances.

###  3.6 Summary of metrics of goodness of fit

Table 6: Summary of the model prediction errors and estimated metric scores.
The first two columns summarize prediction errors used to compare predictive
accuracy of $W^{\prime}_{\mathrm{bart}}$ and
$\mathrm{hydraulic}_{\mathrm{weig}}$ via MAE, standard deviation of absolute
errors (SD), and RMSE. Prediction accuracy had to be be assessed using data to
which models were not fitted to, such that we had to exclude the Bartram and
Weigend data sets. The subsequent three columns summarize the prediction
errors used to compare the predictive accuracies of
$W^{\prime}_{\mathrm{skib}}$, $W^{\prime}_{\mathrm{weig}}$ and
$\mathrm{hydraulic}_{\mathrm{weig}}$. Here, the Weigend data set was excluded
because $W^{\prime}_{\mathrm{weig}}$ and $\mathrm{hydraulic}_{\mathrm{weig}}$
were fitted to it. Finally, the $\mathrm{AIC}_{\mathrm{c}}$ metric requires
models to be fitted to and to be evaluated on the same data. Therefore, we
compared $\mathrm{AIC}_{\mathrm{c}}$ scores estimated from prediction errors
of $W^{\prime}_{\mathrm{weig}}$ and $\mathrm{hydraulic}_{\mathrm{weig}}$ on
all data sets. We approximated p-values for absolute differences in MAE and
RMSE with a bootstrap hypothesis tests and considered $\mathrm{p}<.05$ as
significant. For every metric, a lower score means a better result.

| data for prediction scores 1 | data for prediction scores 2 |  data for $\mathrm{AIC}_{\mathrm{c}}$ scores
---|---|---|---
| $W^{\prime}_{\mathrm{bart}}$ | $\mathrm{hydraulic}_{\mathrm{weig}}$ | $W^{\prime}_{\mathrm{skib}}$ | $W^{\prime}_{\mathrm{weig}}$ | $\mathrm{hydraulic}_{\mathrm{weig}}$ | $W^{\prime}_{\mathrm{weig}}$ | $\mathrm{hydraulic}_{\mathrm{weig}}$
Bartram | 0 |  |  | $0.0$ | $0.0$ | $22.7$ | $0.0$ | $22.7$
1 |  |  | $-20.9$ | $-22.4$ | $0.9$ | $-22.4$ | $0.9$
2 |  |  | $-24.2$ | $-30.1$ | $-2.2$ | $-30.1$ | $-2.2$
3 |  |  | $-24.9$ | $-37.6$ | $-4.3$ | $-37.6$ | $-4.3$
4 |  |  | $-23.7$ | $-44.0$ | $-4.7$ | $-44.0$ | $-4.7$
Caen | 0 | $-0.6$ | $-1.7$ | $-13.1$ | $-19.4$ | $-1.7$ | $-19.4$ | $-1.7$
1 | $13.4$ | $6.4$ | $-6.1$ | $-17.3$ | $6.4$ | $-17.3$ | $6.4$
2 | $29.0$ | $5.6$ | $4.9$ | $-12.3$ | $5.6$ | $-12.3$ | $5.6$
3 | $35.6$ | $2.3$ | $13.2$ | $-6.7$ | $2.3$ | $-6.7$ | $2.3$
4 | $37.7$ | $-0.4$ | $19.0$ | $-1.5$ | $-0.4$ | $-1.5$ | $-0.4$
5 | $39.5$ | $-0.5$ | $24.7$ | $5.0$ | $-0.5$ | $5.0$ | $-0.5$
6 | $26.2$ | $-8.8$ | $22.9$ | $11.7$ | $-8.8$ | $11.7$ | $-8.8$
7 | $28.7$ | $2.5$ | $28.1$ | $23.1$ | $2.5$ | $23.1$ | $2.5$
Chid. | 0 | $25.0$ | $24.0$ | $10.4$ | $-6.0$ | $24.0$ | $-6.0$ | $24.0$
1 | $11.9$ | $9.5$ | $-2.6$ | $-11.3$ | $9.5$ | $-11.3$ | $9.5$
2 | $-3.1$ | $-3.9$ | $-15.2$ | $-17.6$ | $-3.9$ | $-17.6$ | $-3.9$
Ferg. | 0 | $48.8$ | $17.2$ | $28.6$ | $-1.1$ | $17.2$ | $-1.1$ | $17.2$
1 | $34.7$ | $4.4$ | $30.9$ | $8.6$ | $4.4$ | $8.6$ | $4.4$
2 | $14.0$ | $12.4$ | $14.0$ | $10.4$ | $12.4$ | $10.4$ | $12.4$
Weigend | 0 |  |  |  |  |  | $-19.5$ | $3.3$
1 |  |  |  |  |  | $-2.7$ | $4.0$
2 |  |  |  |  |  | $2.6$ | $-0.4$
3 |  |  |  |  |  | $-20.6$ | $-2.5$
4 |  |  |  |  |  | $-6.3$ | $-3.5$
5 |  |  |  |  |  | $5.3$ | $-3.8$
6 |  |  |  |  |  | $-6.5$ | $4.8$
7 |  |  |  |  |  | $6.3$ | $2.0$
8 |  |  |  |  |  | $13.6$ | $0.9$
9 |  |  |  |  |  | $-9.6$ | $0.5$
10 |  |  |  |  |  | $11.2$ | $6.1$
11 |  |  |  |  |  | $13.3$ | $-2.6$
MAE |  $24.87$* | $7.11$ |  $17.23$* |  $15.06$* | $7.07$ |  |
$\pm$ SD |  p m $14.35$ |  p m $6.83$ |  p m $9.34$ |  p m $12.19$ |  p m $7.17$ |  |
RMSE |  $28.46$* | $9.69$ |  $19.48$* |  $19.17$* | $9.94$ |  |
AIC |  |  |  |  |  | $181.03$ | $151.85$

  * $*$

significantly different to $\mathrm{hydraulic}_{\mathrm{weig}}$ predictions

Table 6 summarizes the prediction errors of the competing models and resulting
metric scores on our investigated data sets. RMSE and MAE were defined as the
metrics to assess predictive accuracy. Their MAE scores were 24.87 with a
standard deviation of absolute errors (SD) of 14.35 for
$W^{\prime}_{\mathrm{bart}}$ and 7.11(SD=6.83) for
$\mathrm{hydraulic}_{\mathrm{weig}}$ ($\mathrm{p}<.001$ for the difference in
MAE s, bootstrap hypothesis test). The RMSE scores on Caen, Chidnok, and
Ferguson data sets were 28.46 for $W^{\prime}_{\mathrm{bart}}$ and 9.69 for
$\mathrm{hydraulic}_{\mathrm{weig}}$. Also the bootstrap hypothesis test with
the absolute difference in RMSE s as its test statistic resulted in
$\mathrm{p}<.001$.

Both remaining models $W^{\prime}_{\mathrm{skib}}$ and
$W^{\prime}_{\mathrm{weig}}$ could be compared to
$\mathrm{hydraulic}_{\mathrm{weig}}$ on the Bartram, Caen, Chidnok, and
Ferguson data sets. The $\mathrm{hydraulic}_{\mathrm{weig}}$ featured the
lowest MAE with 7.07, the lowest SD with 7.17, and lowest RMSE with 9.94.
$W^{\prime}_{\mathrm{skib}}$ predictions were significantly different to
$\mathrm{hydraulic}_{\mathrm{weig}}$ ($\mathrm{p}<.001$ with the MAE test
statistic and $\mathrm{p}=.001$ with RMSE). $W^{\prime}_{\mathrm{weig}}$
predictions were significantly different to
$\mathrm{hydraulic}_{\mathrm{weig}}$ ($\mathrm{p}=.019$ with the MAE test
statistic and $\mathrm{p}=.031$ with RMSE).

$\mathrm{AIC}_{\mathrm{c}}$ was chosen as the metric to assess which model
provides the best trade-off between predictive capabilities and complexity.
Models must be fitted to and tested on the same data for
$\mathrm{AIC}_{\mathrm{c}}$ scores to be comparable. Hence, as reflected in
the last two columns of Table 6, $W^{\prime}_{\mathrm{weig}}$ and
$\mathrm{hydraulic}_{\mathrm{weig}}$ could be compared on combined data points
of all covered data sets. With a $k$ of 3 for $W^{\prime}_{\mathrm{weig}}$ and
a $k$ of 8 for $\mathrm{hydraulic}_{\mathrm{weig}}$ the resulting scores were
151.85 for $\mathrm{hydraulic}_{\mathrm{weig}}$ and 181.03 for
$W^{\prime}_{\mathrm{weig}}$. The $\mathrm{hydraulic}_{\mathrm{weig}}$
achieved the lower $\mathrm{AIC}_{\mathrm{c}}$ score.

##  4 Discussion

In this study, we compared the prediction capabilities and goodness of fit of
$\mathrm{hydraulic}_{\mathrm{weig}}$ to that of $W^{\prime}_{\mathrm{bal}}$
models. We hypothesized that the hydraulic model would more accuratly predict
observed recovery ratios observed in past studies. Models were compared on
extracted data from five studies and the $\mathrm{hydraulic}_{\mathrm{weig}}$
model outperformed the $W^{\prime}_{\mathrm{skib}}$,
$W^{\prime}_{\mathrm{bart}}$, and $W^{\prime}_{\mathrm{weig}}$ models with
respect to objective RMSE, MAE, and $\mathrm{AIC}_{\mathrm{c}}$ metrics. Our
findings therefore support the hypothesis. We discuss below our results in
more detail and interpret them in context of findings of previous literature.
We present arguments for why the $\mathrm{hydraulic}_{\mathrm{weig}}$
outperformed the $W^{\prime}_{\mathrm{bal-ode}}$ models and we propose
limitations and future work. Finally, we end this section with statements
about significance and implications of our results.

###  4.1 Interpretation and contextualization

We observed that the standard deviations of absolute prediction errors in
Section 3.6 as well as the overall MAE and RMSE were considerably lower for
the $\mathrm{hydraulic}_{\mathrm{weig}}$ than for the
$W^{\prime}_{\mathrm{bal-ode}}$ models. But when averaging the prediction
errors on isolated data sets listed in Table 6,
$\mathrm{hydraulic}_{\mathrm{weig}}$ only made more accurate predictions than
its competitors on the Bartram and Caen data sets. For the Bartram data set,
the MAE of $\mathrm{hydraulic}_{\mathrm{weig}}$ was 6.96, compared to 18.74
for $W^{\prime}_{\mathrm{skib}}$, and 26.82 for $W^{\prime}_{\mathrm{weig}}$
respectively. For the Caen data set, the MAE of
$\mathrm{hydraulic}_{\mathrm{weig}}$ was the lowest with 3.52. On the
remaining Chidnok data set it was $W^{\prime}_{\mathrm{skib}}$ that achieved
the lowest MAE with 9.4 and on the Ferguson data set it was
$W^{\prime}_{\mathrm{weig}}$ with 6.7.

As pointed out by Skiba and Clarke, (2021) and Sreedhara et al., (2019),
$W^{\prime}_{\mathrm{bal-ode}}$ models are meant to be applied to any athlete
on a wide range of possible conditions. A lower MAE score for
$W^{\prime}_{\mathrm{weig}}$ on the Ferguson data set means
$W^{\prime}_{\mathrm{weig}}$ predicted recovery ratios more closely for the
particular group (six recreational active men) under the particular test
conditions that Ferguson tested. However, to determine the usefulness of a
model for predicting performance for high-intensity intermittent exercise in a
more general sense, models have to be evaluated on a multitude of scenarios.
After combining all data sets, $\mathrm{hydraulic}_{\mathrm{weig}}$ achieved
the overall lowest MAE score, which means that
$\mathrm{hydraulic}_{\mathrm{weig}}$ could predict recovery ratios overall
more accurately for a range of groups and settings.

Figure 9: Simulated recovery ratios using the
$\mathrm{hydraulic}_{\mathrm{weig}}$ and $W^{\prime}_{\mathrm{bal-ode}}$
models in response to prior exercise of differing intensities. The plots show
that $W^{\prime}_{\mathrm{bal-ode}}$ models are insensitive to the properties
of prior exhausting exercise, i.e., their predictions were not affected by
$P_{\mathrm{work}}$. In contrast, the $\mathrm{hydraulic}_{\mathrm{weig}}$ was
sensitive to the prior exercise properties. Performance models
$W^{\prime}_{\mathrm{bal-ode}}$ were configured with a $CP$ = 393,
$W^{\prime}$ = 23300 and $\mathrm{hydraulic}_{\mathrm{weig}}$ featured the
configuration $[$23111.91, 65845.28, 391.57, 148.88, 24.15, 0.73, 0.01,
0.24$]$. All simulations differed only in $P_{\mathrm{work}}$, which decreased
from $P100$ to $P480$ in the simulations depicted from left to right.
$P_{\mathrm{work}}$ = $P100$, prescribed by Bartram et al., (2018), was the
highest intensity out of compared studies. They investigated recovery after 60
seconds, therefore the $W^{\prime}_{\mathrm{bart}}$ model prediction after 60
seconds is marked as the observation. $P_{\mathrm{work}}$ = $P480$ was the
lowest prescribed intensity out of compared studies and recovery ratios for
$P_{\mathrm{work}}$ = $P480$, $P_{\mathrm{rec}}$ = $33\%$ of $CP$ from the
Weigend data set were marked as observations on the right.

The less consistent prediction quality across data sets of the
$W^{\prime}_{\mathrm{bal-ode}}$ models agrees with findings by Caen et al.,
(2019), who proposed that the predictive capabilities of
$W^{\prime}_{\mathrm{bal}}$ models may improve with modifications that account
for intensity and duration of prior exhaustive exercise. As an example, out of
all compared studies in this work, Bartram et al., (2018) prescribed the
highest work bout intensity for their experimental setup ($P_{\mathrm{work}}$
= $P100$). Considering the suggestion by Caen et al., (2019) that a shorter
time to exhaustion at a high intensity allows a quicker recovery, it seems
reasonable that the $W^{\prime}_{\mathrm{bart}}$ model estimated the fastest
recovery kinetics out of all recovery models.

Conversely, the Caen et al., (2019) study prescribed the lowest work bout
intensity out of all compared studies ($P_{\mathrm{work}}$ = $P480$). Their
observed recovery ratios are summarized in the Weigend data set and were
slower than the $W^{\prime}_{\mathrm{bart}}$ predictions. This observation
again matches the assumption that a longer exhaustive exercise at a lower
intensity requires a longer recovery.

Despite the differences in observed recovery rates, the
$W^{\prime}_{\mathrm{bal-ode}}$ models allow for only a single recovery rate
no matter the nature of the prior exercise. To illustrate this point, we
conducted simulations to depict the influence of prior exercise intensity on
the recovery ratios predicted by the $W^{\prime}_{\mathrm{bal-ode}}$ and
$\mathrm{hydraulic}_{\mathrm{weig}}$. Figure 9 depicts four simulations. All
simulations shared the same test setup except for differing
$P_{\mathrm{work}}$ intensities. The simulation on the left had
$P_{\mathrm{work}}$ = $P100$ as prescribed by Bartram et al., (2018), the
simulation on the right featured the lowest $P_{\mathrm{work}}$ = $P480$ as
found in the Weigend data set. From left to right, $P_{\mathrm{work}}$ of the
simulations decreased step wise. Bartram et al., (2018) investigated recovery
after 60 seconds and therefore the $W^{\prime}_{\mathrm{bart}}$ prediction
after 60 seconds is marked as the observation on the left. Recovery ratios
with $P_{\mathrm{work}}$ = $P480$ and $P_{\mathrm{rec}}$ = $33\%$ of $CP$ of
the Weigend data set are marked as observations on the right. The recovery
ratios predicted by the $W^{\prime}_{\mathrm{bal-ode}}$ models were the same
for each $P_{\mathrm{work}}$ and their predictions were therefore unable to
fit all observations equally well. In contrast, the hydraulic model could
account for such characteristics.

This result occurred because of the interactions between the three tanks that
$\mathrm{hydraulic}_{\mathrm{weig}}$ uses to model energy recovery. For
example, during high-intensity exercise, the liquid level in $AnF$ would
rapidly decrease and the contribution of $AnS$ would be less than during
lower-intensity exercise when the liquid level in $AnF$ would decrease more
slowly. Differences in fill states of $AnS$ affected recovery estimations and
enabled $\mathrm{hydraulic}_{\mathrm{weig}}$ to predict rapid recovery after
high-intensity exercise and a slower recovery after exercise at a lower
intensity. We suggest that standard deviations of MAE as well as overall MAE
and RMSE scores of $\mathrm{hydraulic}_{\mathrm{weig}}$ model were smaller
than those of $W^{\prime}_{\mathrm{bal-ode}}$ models because the hydraulic
model could account for characteristics of prior exhaustive exercise.

Further, we propose that $\mathrm{hydraulic}_{\mathrm{weig}}$ has also
achieved overall better metric scores because it better captured the bi-
exponential nature of energy recovery, as opposed to the mono-exponential
$W^{\prime}_{\mathrm{bal-ode}}$ models. Indeed, the observed recovery ratios
of the Caen data set increased rapidly from 0 seconds to 120 seconds and then
continued to rise more slowly at longer durations (Figure 5). Caen et al.,
(2021) showed that their observations were well explained with a bi-
exponential model that implements a steeper slope during the beginning of
recovery. Also, the first $W^{\prime}_{\mathrm{bal}}$ paper by Skiba et al.,
(2012) proposed an alternative bi-exponential version of their
$W^{\prime}_{\mathrm{bal}}$ model with two $\mathcal{T}_{t}$s but such bi-
exponential $W^{\prime}_{\mathrm{bal}}$ models have yet to be applied in
practice.

###  4.2 Limitations and future work

The observed improved prediction quality of
$\mathrm{hydraulic}_{\mathrm{weig}}$ comes at the cost of a time-demanding
fitting process. As outlined in Section 2.1.2, our fitting process from
Weigend et al., (2021) requires $CP$ and $W^{\prime}$ of an athlete as inputs
and then obtains fitted $\mathrm{hydraulic}_{\mathrm{weig}}$ configurations
using an evolutionary computation approach. Different $CP$ and $W^{\prime}$
values require a new $\mathrm{hydraulic}_{\mathrm{weig}}$ configuration to be
fitted to them. Additionally, fittings for our comparison took computation
times of 3 hours or more on 7 cores of an
$\text{Intel}^{\text{\textregistered}}$ $\text{Xeon}^{\text{\textregistered}}$
CPU E5-2650 v4 @ 2.20GHz each. On the other hand, obtaining $\mathcal{T}_{t}$
for our $W^{\prime}_{\mathrm{bal-ode}}$ models was solved in milliseconds and
they can be applied to any $CP$ and $W^{\prime}$ combination without fitting
$\mathcal{T}_{t}$ anew. Therefore, the application of
$\mathrm{hydraulic}_{\mathrm{weig}}$ is a time consuming task in comparison to
the application of $W^{\prime}_{\mathrm{bal-ode}}$ models. In order to improve
the feasibility of application, future work must optimize the hydraulic model
fitting process to minimize this limitation.

Further, prediction results show that estimated recovery ratios of the used
WB1 $\rightarrow$ RB $\rightarrow$ WB2 protocol are affected by rounding
errors that arose from the recovery ratios being estimated from simulations in
discrete time steps. For example, in Table 5, it was observed that the
$W^{\prime}_{\mathrm{bal-ode}}$ model predictions between the $P240$ and
$P480$ trials varied slightly in a few cases, even though these the recovery
kinetics should have been the same. The variations were caused by simulation
time steps of size 0.1 that made rounding effects play a bigger role in
shorter work bouts of $P240$. Smaller step sizes would decrease the error, but
to more fully prevent inaccuracies, future work ought to formalize model
simulations in differential equations that don’t require estimations in
discrete time steps.

Predicted recovery ratios of the $\mathrm{hydraulic}_{\mathrm{weig}}$ at
recovery intensities close to $CP$ require further investigation too. The
comparison on the $D_{CP}$0 case of the Bartram data set (See Figure 4)
revealed that $\mathrm{hydraulic}_{\mathrm{weig}}$ predicts a slight recovery
during exercise at $CP$. No liquid from $Ae$ flowed back into the system but
the ongoing flow from $AnS$ to $AnF$ still caused liquid level in $AnF$ to
rise during the recovery bout. Recovery while exercising at $CP$ intensity is
a controversial assumption that is made to an even stronger extent by the
original $W^{\prime}_{\mathrm{bal}}$ model of Skiba et al., (2012). Such
dynamics have to be taken into consideration when the models are used for
predictions and are important directions for future investigation.

Additionally, extracted recovery ratio observations from previous studies come
with associated uncertainties that could not be considered in our MAE, RMSE
and $\mathrm{AIC}_{\mathrm{c}}$ scores. As an example, the standard deviations
of observed recovery ratios by Caen et al., (2021) depicted in Figure 5 are
greater than reported standard deviations by Ferguson et al., (2010) depicted
in Figure 7. One could argue that comparisons to reported means by Ferguson et
al., (2010) therefore provide a better indication of prediction quality.
Unfortunately, we could not incorporate these standard deviations into
goodness-of-fit metrics because of how different recovery ratios were reported
in compared studies. Caen et al., (2021) and Ferguson et al., (2010) reported
averaged observed recovery ratios with standard deviations, for the Bartram
data set we had to use $W^{\prime}_{\mathrm{bart}}$ predictions as
observations to compare to, in Weigend et al., (2021) we derived our values
from Caen et al., (2019) without standard deviations, and for the Chidnok data
set we fitted constant values for $\mathcal{T}_{t}$ for
$W^{\prime}_{\mathrm{bal-ode}}$ models to their reported times to exhaustion
to obtain comparable recovery ratios in percent. We believe these various
formats of observations highlight the need for more and more comparable
studies on energy recovery dynamics.

Larger data sets are vital for more educated investigations of recovery models
and their improvement in future work. We see the combination of data sets in
this work as a step towards this direction. In order to improve and compare
models more holistically, it is important that more comparable studies are
conducted in the future and combined into a larger test bed for performance
models.

###  4.3 Significance and implications

To the best of our knowledge, performance models on energy recovery during
intermittent exercise have yet to be compared in such detail. Our comparison
on data from five studies allowed a more holistic view on recovery dynamics
and confirmed limitations of $W^{\prime}_{\mathrm{bal-ode}}$ models that were
suggested by previous literature. Our results imply that more complex models
like $\mathrm{hydraulic}_{\mathrm{weig}}$ can improve energy recovery
predictions. We propose that further efforts to merge and compare data are
significant steps to bring the research area of energy recovery modeling
forward. We further propose that the predictive capabilities of hydraulic
models look strong and that hydraulic models have to be considered as a
possible future direction to advance energy recovery modeling.

##  5 Conclusion

We conclude that the $\mathrm{hydraulic}_{\mathrm{weig}}$ outperformed
$W^{\prime}_{\mathrm{bal-ode}}$ models when fit to multiple independent data
sets featuring intermittent high-intensity exercise. The predictive accuracy
and goodness of fit was better for the hydraulic model, even for the
$\mathrm{AIC}_{\mathrm{c}}$ metric, which includes a penalty for the number of
model parameters. The hydraulic model is thus likely more generalizable than
$W^{\prime}_{\mathrm{bal-ode}}$ models, which are typically applied within
narrow contexts. Future research should focus on improving the feasibility of
the hydraulic model, because it is computationally more burdensome to use than
$W^{\prime}_{\mathrm{bal-ode}}$ models. Pending such improvements, we foresee
athletes adopting this model to optimize pacing and interval training
workouts. To contribute towards further advancements we publish all material,
extracted data, and simulation scripts here:
https://github.com/faweigend/pypermod.

## Declarations

##### Funding

No funding was received to assist with the preparation of this manuscript.

##### Conflict of interest

The authors have no conflicts of interest to declare that are relevant to the
content of this article.

##### Availability of data and material

All material and extracted data are summarized in Tables 1, 2, 3, 4 and 5 and
are also available here: https://github.com/faweigend/pypermod.

##### Code availability

All code is thoroughly documented and available as source code and as a python
package here: https://github.com/faweigend/pypermod.

## References

  * Bartram et al., (2018) Bartram, J. C., Thewlis, D., Martin, D. T., and Norton, K. I. (2018).  Accuracy of W’ recovery kinetics in high performance cyclists - modeling intermittent work capacity.  International Journal of Sports Physiology and Performance, 13(6):724–728.
  * Behncke, (1997) Behncke, H. (1997).  Optimization models for the force and energy in competitive running.  Journal of Mathematical Biology, 35(4):375–390.
  * Binder et al., (2008) Binder, R. K., Wonisch, M., Corra, U., Cohen-Solal, A., Vanhees, L., Saner, H., and Schmid, J.-P. (2008).  Methodological approach to the first and second lactate threshold in incremental cardiopulmonary exercise testing.  European Journal of Cardiovascular Prevention & Rehabilitation, 15(6):726–734.
  * Burnham and Anderson, (2004) Burnham, K. P. and Anderson, D. R. (2004).  Multimodel Inference: Understanding AIC and BIC in Model Selection.  Sociological Methods & Research, 33(2):261–304.
  * Caen et al., (2021) Caen, K., Bourgois, G., Dauwe, C., Blancquaert, L., Vermeire, K., Lievens, E., Van Dorpe, J., Derave, W., Bourgois, J. G., Pringels, L., and Boone, J. (2021).  W’ Recovery Kinetics following Exhaustion: A Two-Phase Exponential Process Influenced by Aerobic Fitness.  Medicine & Science in Sports & Exercise, Publish Ahead of Print.
  * Caen et al., (2019) Caen, K., Bourgois, J. G., Bourgois, G., Van Der Stede, T., Vermeire, K., and Boone, J. (2019).  The reconstitution of W’ depends on both work and recovery characteristics.  Medicine & Science in Sports & Exercise, 51(8):1745–1751.
  * Chai and Draxler, (2014) Chai, T. and Draxler, R. R. (2014).  Root mean square error (RMSE) or mean absolute error (MAE)? – Arguments against avoiding RMSE in the literature.  Geoscientific Model Development, 7(3):1247–1250.
  * Chidnok et al., (2012) Chidnok, W., Dimenna, F. J., Bailey, S. J., Vanhatalo, A., Morton, R. H., Wilkerson, D. P., and Jones, A. M. (2012).  Exercise tolerance in intermittent cycling: application of the critical power concept.  Medicine & Science in Sports & Exercise, 44(5):966–976.
  * Chorley and Lamb, (2020) Chorley, A. and Lamb, K. L. (2020).  The Application of Critical Power, the Work Capacity above Critical Power (W’), and Its Reconstitution: A Narrative Review of Current Evidence and Implications for Cycling Training Prescription.  Sports, 8(9):123.
  * de Jong et al., (2017) de Jong, J., Fokkink, R., Olsder, G. J., and Schwab, A. (2017).  The individual time trial as an optimal control problem.  Proceedings of the Institution of Mechanical Engineers, Part P: Journal of Sports Engineering and Technology, 231(3):200–206.
  * Efron and Tibshirani, (1993) Efron, B. and Tibshirani, R. (1993).  An introduction to the bootstrap.  Number 57 in Monographs on statistics and applied probability. Chapman & Hall, New York.
  * Ferguson et al., (2010) Ferguson, C., Rossiter, H. B., Whipp, B. J., Cathcart, A. J., Murgatroyd, S. R., and Ward, S. A. (2010).  Effect of recovery duration from prior exhaustive exercise on the parameters of the power-duration relationship.  Journal of Applied Physiology, 108(4):866–874.
  * Good, (2000) Good, P. I. (2000).  Permutation tests: a practical guide to resampling methods for testing hypotheses.  Springer, New York.  OCLC: 681912126.
  * Hill, (1993) Hill, D. W. (1993).  The critical power concept: a review.  Sports Medicine, 16(4):237–254.
  * Hoogkamer et al., (2018) Hoogkamer, W., Snyder, K. L., and Arellano, C. J. (2018).  Modeling the benefits of cooperative drafting: is there an optimal strategy to facilitate a sub-2-hour marathon performance?  Sports Medicine, 48(12):2859–2867.
  * Jones and Vanhatalo, (2017) Jones, A. M. and Vanhatalo, A. (2017).  The ’critical power’ concept: applications to sports performance with a focus on intermittent high-intensity exercise.  Sports Medicine, 47(S1):65–78.
  * Margaria, (1976) Margaria, R. (1976).  Biomechanics and energetics of muscular exercise.  Oxford University Press, Oxford University Press, Walton Street, Oxford, OX2 6DP.
  * Monod and Scherrer, (1965) Monod, H. and Scherrer, J. (1965).  The work capacity of a synergetic muscular group.  Ergonomics, 8(3):329–338.
  * Morton, (1986) Morton, R. H. (1986).  A three component model of human bioenergetics.  Journal of Mathematical Biology, 24(4):451–466.
  * Morton, (2006) Morton, R. H. (2006).  The critical power and related whole-body bioenergetic models.  European Journal of Applied Physiology, 96(4):339–354.
  * Morton and Billat, (2004) Morton, R. H. and Billat, L. V. (2004).  The critical power model for intermittent exercise.  European Journal of Applied Physiology, 91(2-3):303–307.
  * Poole et al., (2016) Poole, D. C., Burnley, M., Vanhatalo, A., Rossiter, H. B., and Jones, A. M. (2016).  Critical power: an important fatigue threshold in exercise physiology.  Medicine & Science in Sports & Exercise, 48(11):2320–2334.
  * SciPy 1.0 Contributors et al., (2020) SciPy 1.0 Contributors, Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S. J., Brett, M., Wilson, J., Millman, K. J., Mayorov, N., Nelson, A. R. J., Jones, E., Kern, R., Larson, E., Carey, C. J., Polat, i., Feng, Y., Moore, E. W., VanderPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E. A., Harris, C. R., Archibald, A. M., Ribeiro, A. H., Pedregosa, F., and van Mulbregt, P. (2020).  SciPy 1.0: fundamental algorithms for scientific computing in Python.  Nature Methods, 17(3):261–272.
  * Skiba et al., (2012) Skiba, P. F., Chidnok, W., Vanhatalo, A., and Jones, A. M. (2012).  Modeling the expenditure and reconstitution of work capacity above critical power.  Medicine & Science in Sports & Exercise, 44(8):1526–1532.
  * Skiba and Clarke, (2021) Skiba, P. F. and Clarke, D. C. (2021).  The W’ Balance Model: Mathematical and Methodological Considerations.  International Journal of Sports Physiology and Performance, 16(11):1561–1572.
  * Skiba et al., (2014) Skiba, P. F., Clarke, D. C., Vanhatalo, A., and Jones, A. M. (2014).  Validation of a novel intermittent W’ model for cycling using field data.  International Journal of Sports Physiology and Performance, 9(6):900–904.
  * Skiba et al., (2015) Skiba, P. F., Fulford, J., Clarke, D. C., Vanhatalo, A., and Jones, A. M. (2015).  Intramuscular determinants of the ability to recover work capacity above critical power.  European Journal of Applied Physiology, 115(4):703–713.
  * Sreedhara et al., (2019) Sreedhara, V. S. M., Mocko, G. M., and Hutchison, R. E. (2019).  A survey of mathematical models of human performance using power and energy.  Sports Medicine - Open, 5(1):54.
  * Sugiura, (1978) Sugiura, N. (1978).  Further analysts of the data by akaike’ s information criterion and the finite corrections: Further analysts of the data by akaike’ s.  Communications in Statistics - Theory and Methods, 7(1):13–26.
  * Sundström, (2016) Sundström, D. (2016).  On a bioenergetic four-compartment model for human exercise.  Sports Engineering, 19(4):251–263.
  * Sundström et al., (2014) Sundström, D., Carlsson, P., and Tinnsten, M. (2014).  Comparing bioenergetic models for the optimisation of pacing strategy in road cycling.  Sports Engineering, 17(4):207–215.
  * Vanhatalo et al., (2011) Vanhatalo, A., Jones, A. M., and Burnley, M. (2011).  Application of critical power in sport.  International Journal of Sports Physiology and Performance, 6(1):128–136.
  * Weigend et al., (2021) Weigend, F. C., Siegler, J., and Obst, O. (2021).  A new pathway to approximate energy expenditure and recovery of an athlete.  In Proceedings of the Genetic and Evolutionary Computation Conference Companion, pages 325–326, Lille France. ACM.
  * Whipp et al., (1982) Whipp, B. J., Huntsman, D. J., Storer, T. W., Lamarra, N., and Wasserman, K. (1982).  A constant which determines the duration of tolerance to high-intensity work.  Federation proceedings, 41(5):1591–1591.