MetaboliSimv1.0

Interactive simulation & analysis of human skeletal-muscle energy metabolism – Mader (2003), Heck et al. (2022).

Athlete Profile
Auto: Cycling ks4 = 12.322, Running ks4 = 9.858, VO₂base = 5 × BM. Manual: edit in Advanced Settings.
Advanced Settings
Economy Parameters
Auto: ks4 = 12.322 (Cycling) / 9.858 (Running), VO₂base = 5 × BM ml/min. Manual mode: freely editable.
Running Economy
C_run: running economy – lower = more economical (typical 0.24–0.32).
Peak Power (Pmax)
Measured 1 s peak power → calibrates η(0) for dynamic efficiency (Dunst 2023). Enables APR KPI.
Environment
VO₂max −6.3%/1000 m > 580 m · VO₂max −1.5%/°C > 22°C · VLamax +1.0%/°C > 22°C
Starting Conditions
Simulation
Pushes values to Simulation, CP, τ, W′bal & Configuration.
Key Metrics
Analysis Modules

Simulation

Constant-load, interval, ramp & sprint. Full metabolic state: ATP, PCr, pH, lactate, glycogen.

CORE

Steady State

Algebraic SS: VO₂ss, vLass, substrate utilization, fat/CHO oxidation.

ANALYSIS

Metabolic Rates

Oxidative, glycolytic, PCr pathway rates. Michaelis–Menten kinetics.

CORE

Critical Power

P–D models: 2P, 3P, Omni 5P, Bi-Exp (6P). W′ balance & τ recovery analysis.

ANALYSIS

Model Configuration

Initial metabolic state (PCr, La, VO₂) and all 30+ Mader model constants. Customize or reset.

SETTINGS
Scientific Background
Model Theory & References
Mader Model (2003)

Coupled ODE system (5 state variables) for skeletal-muscle energy metabolism. Three ATP-resynthesis pathways: oxidative phosphorylation, anaerobic glycolysis, PCr hydrolysis. Regulation via [ADP], [H+], and glycogen fill level.

dPCr/dt = (Vox + Vgly − Vdemand) − b · vLa,res
dVO2/dt = (VO2,ss − VO2) · kVO2
dLam/dt = vLa(ADP, pH) − vdiff(Lam, Lab)
dLab/dt = vdiff − vox(Lab)
dGly/dt = −vLa · cgly
Glycogen Dynamics (Heck et al. 2022)

Finite muscle glycogen store with intensity-dependent depletion. Modulates VLamax via sigmoidal Hill function (Eq. 34, nH ≈ 3). VO2max declines with falling glycogen (Eq. 35). Enables endurance-limited scenarios below CP.

Eq. 34:  f(r) = rnH / (KS4nH + rnH)
Eq. 35:  VO2maxeff = VO2max · [a + (1−a) · r1/3]
r = Gly / Glyfull
Dynamic Efficiency (Dunst et al. 2023)

Time- and fatigue-dependent mechanical efficiency η(t). ATP-mix weighted: at exercise onset, PCr dominates (ηPCr ≈ 0.10); as aerobic share rises, η → ηox. Steady-state ηox is individualised via ks4 (ΔVO2/ΔP). ηPCr calibrated from Pmax (1 s).

η(t) = ∑(V̇i · ηi) / ∑ V̇i
ηox = 60 / (ks4 · CEfat)   ηgly ≈ 0.15   ηPCr ≈ 0.10
CEfat = 19.584 J/mL O2
W′ Balance Models

Two complementary approaches to track the remaining anaerobic work capacity (W′bal) during variable-power exercise. Integral model (Skiba et al. 2012): convolution of expenditure history with exponential recovery kernel (τW′). Differential model (Skiba et al. 2015): piece-wise ODE — depletion above CP, recovery below CP with DCP-dependent rate. Clarke & Skiba (2021) showed the ODE variant better predicts intermittent exercise.

INT:  W′bal(t) = W′ − ∑ W′exp,i · e−(t−ti)/τ
ODE:  dW′bal/dt = (P − CP)   if P > CP
     dW′bal/dt = (W′ − W′bal) · DCP   if P ≤ CP
Maximum Power Available (Kontro et al. 2026)

3-parameter CP model extends the classical 2P model with Pmax (maximal instantaneous power). MPA declines linearly with W′ expenditure: when fresh, MPA = Pmax; at full depletion, MPA = CP. Task failure occurs when P = MPA, not when W′ = 0. The strain coefficient kstrain amplifies training load as MPA approaches P, enabling energy-system-specific training quantification (SSCP, SSW′, SSPmax).

3P:  tlim = W′/(P−CP) − W′/(Pmax−CP)
MPA = Pmax − (Pmax − CP) · W′exp / W′
kstrain = (Pmax − MPA + CP) / (Pmax − P + CP)
SS = ∑ kstrain · P · (Pmax / CP²) · 100/3600
VO2 Slow Component

Fatigue-driven loss of contractile efficiency above MLSS. Metabolite accumulation ([Pi]↑, [H+]↑, ΔGATP↓) reduces cell excitability in active fibres; additionally, progressive recruitment of less efficient type-II motor units (Jones et al. 2011). Reduces Tlim near PVO2max by 20–40 %.

SC(t) = 1 + kSC · f(P) · tsupra
f(P) = (P − PMLSS) / (PVO2max − PMLSS)
f(P) ∈ [0, 3]

Key References

  • Mader A (2003). Glycolysis and oxidative phosphorylation as a function of cytosolic phosphorylation state and power output of the muscle cell. Eur J Appl Physiol 88(4–5):317–338.
  • Mader A, Heck H (1994). Energiestoffwechselregulation, Erweiterungen des theoretischen Konzepts und seiner Begründungen. In: Mader A, Allmer H (Eds.), Computersimulation. Academia Verlag, Sankt Augustin.
  • Heck H, Bartmus U, Grabow V (2022). Laktat. Springer.
  • Dunst AK, Hesse C, Ueberschär O, Holmberg HC (2023). Sports (Basel) 11(2):29.
  • Monod H, Scherrer J (1965). Ergonomics 8:329–338.
  • Skiba PF, Clarke DC (2021). Int J Sports Physiol Perform 16(11):1561–1572.
  • Skiba PF, Chidnok W, Vanhatalo A, Jones AM (2012). Modeling the expenditure and reconstitution of W′ during intermittent exercise. Med Sci Sports Exerc 44(8):1526–1532.
  • Skiba PF, Fulford J, Clarke DC, Vanhatalo A, Jones AM (2015). Intramuscular determinants of the ability to recover work capacity above critical power. Eur J Appl Physiol 115(4):703–713.
  • Morton RH (1996). A 3-parameter critical power model. Ergonomics 39(4):611–619.
  • Kontro H, Mastracci A, Cheung SS, MacInnis MJ (2026). The three-dimensional impulse-response model: Modeling the training process in accordance with energy system-specific adaptation. PLoS One 21(2):e0341721.
  • Jones AM, Grassi B, Christensen PM, Krustrup P, Bangsbo J, Poole DC (2011). Slow component of VO₂ kinetics: mechanistic bases and practical applications. Med Sci Sports Exerc 43(11):2046–2062.
Simulation Results
Import / Export
Two columns: time (s) and power (W). Header optional. Comma / semicolon / tab.
Summary

                              
Scenarios

Dynamic Metabolic Model

This simulator integrates the nonlinear ODE system proposed by Mader (2003) for cellular energy metabolism, extended by the two-compartment lactate model of Heck, Bartmus & Grabow (2022). The state vector comprises four coupled variables — phosphocreatine [PCr], pulmonary oxygen uptake V̇O2, muscle lactate [Lam], and blood lactate [Lab] — whose dynamics are governed by the creatine-kinase / adenylate-kinase equilibrium (Mader 2003, Eq. 1–6) and Michaelis–Menten-type rate equations for oxidative phosphorylation and anaerobic glycolysis.

State equations (simplified):

d[PCr]/dt = V̇ox + V̇gly − V̇ATP,demand

dV̇O2/dt = kVO₂ · (V̇O2,ss(ADP) − V̇O2)

d[Lam]/dt = v̇La − vox,La − Kdif · ΔLa

d[Lab]/dt = Vrel · (Kdif · ΔLa − vox,b − vres,b)

The phosphorylation state ([ATP], [ADP], [AMP], Pi) is recomputed at every time step from the CHEP equilibrium (creatine kinase KCK = 1.66 · 10⁹, adenylate kinase KAK = 1.05, total adenine nucleotide pool TAN = 8.2 mmol/kg; Mader 2003). [ADP] serves as the central feedback signal driving oxidative phosphorylation (sigmoidal Hill kinetics with nox) and glycolytic flux (via phosphofructokinase activation). Integration uses an explicit fourth-order Runge–Kutta scheme (dt typically 0.1 s).

Two-Compartment Lactate Model

Muscle and blood lactate are coupled by a concentration-dependent diffusion term Kdif(Lab) = Kdif,0 · Lab−K₁, representing monocarboxylate transporter (MCT1/MCT4) saturation kinetics (Heck et al. 2022, Eq. 27). Lactate elimination occurs via oxidation (vox = KLaO₂ · [La]) in both compartments and gluconeogenic resynthesis (vres, ADP-driven, primarily hepatic/renal). The volume ratio Vrel = mactive / (BWeff − mactive) scales inter-compartmental fluxes, where BWeff represents the effective distribution volume (default 45 kg for 75 kg body mass).

pH and Glycolytic Inhibition

Intramuscular pH is derived from lactate accumulation and non-bicarbonate buffer capacity β (slyke; Mader & Heck 1994): pH = pHrest − Δ[Lam] / β. Rising [H⁺] inhibits phosphofructokinase (PFK) and thereby glycolytic rate via:

V̇Lamax,pH = V̇Lamax / (1 + [H⁺]n / Ks3)

(Mader & Heck 1994, Eq. 18). This pH-mediated feedback explains the self-limiting nature of supra-maximal exercise and the emergence of a maximal lactate steady state (MLSS) where lactate production equals elimination.

Glycogen Depletion

Muscle glycogen is tracked as a fifth state variable, depleted by glycolytic flux and partially replenished by exogenous CHO intake. Depletion reduces V̇Lamax via a Hill function: f(Glyc) = GlycnH / (KmnH + GlycnH) (Km = 0.25, nH = 3; Neubig 2021, in Heck et al. 2022, Eq. 34), representing substrate limitation of glycogen phosphorylase. V̇O2max is attenuated via a cube-root relationship (minimum factor ≈ 0.50; Hargreaves & Spriet 2020), reflecting reduced pyruvate flux to the TCA cycle when glycogen stores are critically low.

Dynamic Gross Efficiency

Constant η — A fixed gross efficiency (η ≈ 25 %, corresponding to Ks4 = 12.322; Mader 2003) converts metabolic ATP turnover to mechanical power. This is the default mode and appropriate for most aerobic simulations.

Dynamic η(t) — Following Dunst et al. (2023), instantaneous efficiency is computed as the ATP-flux-weighted mean of pathway-specific efficiencies: ηox ≈ 24.9 % (= 60 / (Ks4 · CEfat) with CEfat = 19.584 J/mL O₂ at RER = 0.70), ηgly ≈ 15 %, and ηpcr ≈ 10–13 %. This yields η ≈ 13 % during all-out sprints (PCr-dominated) and η ≈ 23–25 % at aerobic steady state. When a measured 1 s peak power Pmax is provided, ηpcr is calibrated so the model reproduces the observed sprint power.

Dynamic η(t) + SC — Adds the V̇O₂ slow component. Above MLSS, effective efficiency decreases progressively: SC = 1 + ksc · f(I) · tsupra, where f(I) = clamp((P − PMLSS) / (PV̇O₂max − PMLSS), 0, 3) is the normalised supra-threshold intensity and tsupra the accumulated time above MLSS. This increases ATP cost per watt by ~3–8 %/min in the severe domain, shortening Tlim by ≈ 20–40 % near PV̇O₂max. The mechanism reflects progressive type-II fibre recruitment, Pi-mediated impairment of cross-bridge cycling, and mitochondrial inefficiency (Jones et al. 2011; Grassi et al. 2015). Auto-activates for Sprint and Max Effort modes.

Simulation Modes

Abort on Exhaustion terminates when any metabolic threshold is reached (PCr < 1.0 mmol/kg, pH < 6.4, or ATP < 2.0 mmol/kg), then continues at recovery power until tsimulation. Thresholds represent the physiological limits of the Lohmann reaction, intracellular acidosis tolerance, and adenine nucleotide depletion, respectively.

Energy Limited allows the athlete to continue beyond metabolic thresholds but reduces realisable power to the instantaneously available metabolic capacity (Preal = min(Pdemand, Pmax,metabolic)).

Max Effort computes the maximal sustainable power at each time step from instantaneous V̇max and bWatt, requiring dynamic efficiency and Pmax calibration.

Environmental Conditions

Altitude. V̇O2max decreases linearly by 6.3 % per 1 000 m above 580 m due to reduced arterial O₂ saturation (Wehrlin & Hallén 2006). V̇Lamax is unaffected as glycolysis is O₂-independent.

Temperature. Above the thermoneutral zone (≤ 22 °C), V̇O2max decreases by 1.5 %/°C due to thermoregulatory competition for cardiac output and peripheral vasodilation (Périard, Racinais & Sawka 2015). V̇Lamax increases by 1.0 %/°C due to the Q10 effect on phosphofructokinase and catecholamine-driven glycogenolysis (Febbraio et al. 1994). Altitude and temperature effects on V̇O2max are applied multiplicatively.


References

Mader A (2003) Eur J Appl Physiol 88:317–338 · Mader A, Heck H (1994) BSW 8(2):124–162 · Heck H, Bartmus U, Grabow V (2022) Laktat — Das Buch, Kap. 4, Feldhaus · Dunst AK, Hesse C, Ueberschär O, Holmberg HC (2023) Sports 11(2):29 · Jones AM, Grassi B, Christensen PM et al. (2011) Med Sci Sports Exerc 43:2046–2062 · Grassi B, Rossiter HB, Zoladz JA (2015) J Physiol 593:3413–3434 · Wehrlin JP, Hallén J (2006) Eur J Appl Physiol 96:404–412 · Périard JD, Racinais S, Sawka MN (2015) Compr Physiol 5:255–281 · Febbraio MA et al. (1994) J Appl Physiol 77:2827–2831 · Hargreaves M, Spriet LL (2020) Nat Metab 2:817–828

Energy System Contribution Over Time
Description

This plot shows the relative contribution of the three energy systems (oxidative, glycolytic, and PCr) to overall ATP turnover.

The PCr contribution reflects the instantaneous deficit between ATP demand and aerobic + glycolytic supply. At exercise onset or intensity increases, PCr spikes as it buffers the VO₂ kinetics delay.

Steady-State Analysis

Steady-state model: The algebraic steady-state curves (vLass, vLaox, PD, CLass) describe the 1-compartment muscle model at constant pH = 7.0 (Heck et al. 2022).
MLSS is the power at which PD = vLaox − vLass = 0.

La(b) — simulated: La(b) is computed by running a cumulative step test through the full dynamic 2-compartment model. This integrates the complete ODE system including non-linear diffusion
Kdiff(LaB) = 0.065 · LaB−1.4, pH-dependent glycolysis, and VO₂ kinetics (Mader 2003). The step duration determines how close each step approaches its true steady state.

pH correction: PFK activity falls with rising [H⁺]. vLamaxpH = vLamax / (1 + [H⁺]³ / Ks3) (Mader & Heck 1994, Eq. 18). Effect: vLass is reduced → PD increases, MLSS shifts ~5–10 W higher.

Glycogen depletion: Maximum exercise duration is limited by the ratio of usable glycogen (total × allowed depletion %) to net CHO oxidation rate (steady-state CHO use minus exogenous intake). Background curves show different intake scenarios (0–120 g/h).

Glycogen → VO₂max: VO₂maxeff = VO₂max × [a + (1−a) × (Glyc/Glycfull)¹⁄³]. @ Glyc = 0%: VO₂max drops to a × 100% (≈ 35–60% depending on β-oxidation capacity). Calibrated from McArdle disease data (Lucia et al. 2008, N = 46).

Oxidation and Glycolysis Characteristics
CHEP Equilibrium System and −ΔG_ATP
High-Energy Phosphates vs. Metabolites
Current Parameters
Reference Values

Activation Curves: Shows ADP-dependent activation of oxidative phosphorylation (VO₂ss) and glycolysis (vLass). The vLass curves shift toward higher ADP concentrations at lower pH due to pH inhibition of phosphofructokinase (PFK).

CHEP Equilibrium: The coupled adenylate-creatine phosphate system displays thermodynamically linked concentrations of PCr, ATP, ADP, and AMP, along with the Gibbs free energy of ATP hydrolysis (−ΔG_ATP).

Energy State: Alternative representation with high-energy phosphates (PCr + 2·ATP + ADP) on the X-axis. Zones indicate 'Rest' and 'Contractile Insufficiency' states.

Model Comparison
Simulation Details

                        

                        

Monod & Scherrer (1965) postulated a linear relationship between total work and time to exhaustion:

W_lim = W' + CP × t_lim

CP is the Critical Power — the highest sustainable rate of aerobic energy expenditure. W' is the curvature constant often interpreted as anaerobic work capacity (Whipp et al. 1982).

The 2-Parameter and 3-Parameter models are only valid up to ~30 min (CP domain). Data points beyond 30 min are excluded from these fits but included in the Omni model.

W' Balance Over Time
W'bal-ODE (Differential Model)

Piecewise differential equation:

P ≥ CP: dW'bal/dt = −(P − CP) P < CP: dW'bal/dt = (1 − W'bal/W') · (CP − P)

Recovery is state-dependent: the emptier the tank, the faster the reconstitution. No τ parameter needed.


W'bal-INT (Integral/Convolution Model)

Each supra-CP expenditure is 'discounted' by exponential decay:

W'bal[j] = W' − Σᵢ exp(−(j−i)·Δt/τ) · ΔW'exp[i]

τ = 546 · e^(−0.01·DCP) + 316 (Skiba et al. 2012, group fit)

Reference: Skiba PF, Clarke DC (2021). Int J Sports Physiol Perform, 16:1561–1568.

Editable Parameters


Complete Initial State
Thermodynamics
Phosphate Pools
CHEP Equilibrium
Oxidative Phosphorylation
Glycolysis
Lactate Resynthesis
Lactate Kinetics
Buffering & Thresholds
ATP Demand Mapping

Apply or Reset
Scientific Descriptions of Model Constants