Module 4 65–80 min

Dynamic systems

In Module 3, every model was a fixed formula: give it a duration or a time since the start of exercise, and it returns a value. That works as long as the input does one simple thing, like a single step in power. Real training is not like that: power goes up and down, sessions follow each other, and what the body does now depends on what it did before. This module is about models that can follow that history.


From static to dynamic

A static model has no memory: its output depends only on the current value of its input, through a fixed formula. A dynamic model has memory: its output depends on the history of its input, because it carries a state that evolves over time (described by differential equations). The static VO₂ curveM3, for example, is a fixed formula of the time since a single step in power; it cannot respond to a power that keeps changing.

Comes back inModule 3

The chart below shows the difference on real data from an incremental exercise test, where the load goes up in steps of a few minutes each.

Loading fat_max.csv data...

Top panel: a static model. Each dot is the VO₂ measured at the end of a step, when it has settled, plotted against the load of that step. A straight line fits these points well:

VO2SS=VO2RES+ϵ⋅PVO_{2SS} = VO_{2RES} + \epsilon \cdot P

where ϵ\epsilon (the slope) is the oxygen cost of power, in mL/min of O₂ per watt, and VO2RESVO_{2RES} (the intercept) is the VO₂ extrapolated to zero load. This is a perfectly good static model: tell it the load, and it tells you the VO₂ the body will settle at.

Bottom panel: what really happens over time. The gray lines are the measured VO₂ (faint: raw; solid: an 11-point moving average). The teal line is what the static model predicts at each moment from the load at that moment. Hover over a dot in the top panel and an orange line marks its prediction in the time series.

Look at what happens at each step. The static model jumps to the new value immediately; the measured VO₂ does not. It gets there, but it takes a minute or two. The static model knows where VO₂ will settle, but not how and when it gets there. To describe the transition, the model needs memory: the VO₂ of the next second depends on the VO₂ of this second, not only on the load. That is what a dynamic model adds.


Differential equations

Dynamic models are typically represented using differential equations. In mathematics, a differential equation relates one or more unknown functions to their derivatives. A derivative quantifies how a function changes in response to changes in its input.

In our case, we are interested in how a variable evolves over time. This means we are specifically concerned with the time derivative of a variable. We write the derivative term as:

ddty(t)\frac{d}{dt}y(t)

This is read as: "the derivative of y with respect to time (t)." Do not worry too much about the notation or mathematics just yet. The important takeaway is that differential equations, those involving derivatives, help us describe how systems change over time, and thus are key to building dynamic models.

For example, a simple model such as:

ddty(t)=3\frac{d}{dt}y(t)=3

describes a system where a property yy increases at a constant rate of 3.

One of the most important differential equations in sport science, though, is:

τddty(t)+y(t)=x(t)\tau\frac{d}{dt}y(t)+y(t)=x(t)

This equation describes a system that reacts exponentially to changes in the input x(t)x(t). The parameter τ\tau is the time constantM3 of the system, which we already encountered in Module 3! This is one of the most versatile equations in your modeling toolbox 🛠️, believe me. If you get the behavior of this equation, you can model a wide range of systems, such as the oxygen consumption dynamics during generic exercise, the supercompensation principle and the Banister model, exponentially weighted moving averages, a resistance-capacitance circuit, and on and on and on (we will put them side by side in One equation, many systems).

In this module we look at two families of dynamic models. First-order models, like the equation above, will help you understand the importance of the time characteristics of a system. Second-order models will help you understand the importance of inertia, stiffness, and damping.


First-order dynamic systems

The VO₂ kinetics model (dynamic version)

It is highly recommended to review the static version of the VO₂ modelM3 before approaching this dynamic version.

The simplest dynamic formulation of the VO₂ kinetics model is:

τddtVO2(t)+VO2(t)=VO2SS(t)\tau \frac{d}{dt} VO_2(t) + VO_2(t) = VO_{2SS}(t)

This equation describes how the actual oxygen uptake VO2(t)VO_2(t) responds over time to a steady-state target VO2SS(t)VO_{2SS}(t).

The steady-state oxygen uptake is the static model from the beginning of this module:

VO2SS(t)=P(t)⋅ϵ+VO2RESVO_{2SS}(t) = P(t) \cdot \epsilon + VO_{2RES}

Where:

  • P(t)P(t) is the power output at time tt
  • ϵ\epsilon is the oxygen cost per unit of power (e.g., in mL·W⁻¹·min⁻¹)
  • VO2RESVO_{2RES} is the resting or baseline oxygen consumption

Substituting into the first equation, we get:

τddtVO2(t)+VO2(t)=P(t)⋅ϵ+VO2RES\tau \frac{d}{dt} VO_2(t) + VO_2(t) = P(t) \cdot \epsilon + VO_{2RES}

Here, you can recognize the derivative operator (rate of change) and the time constant τ\tau, which governs how fast VO₂ adapts.

To implement this model in a spreadsheet or software, we approximate the derivative using a finite difference method:

ddtVO2(t)≈VO2,k+1−VO2,ktk+1−tk\frac{d}{dt} VO_2(t) \approx \frac{VO_{2,k+1} - VO_{2,k}}{t_{k+1} - t_k}

Substituting into the dynamic equation:

τ⋅VO2,k+1−VO2,ktk+1−tk+VO2,k=Pk⋅ϵ+VO2RES\tau \cdot \frac{VO_{2,k+1} - VO_{2,k}}{t_{k+1} - t_k} + VO_{2,k} = P_k \cdot \epsilon + VO_{2RES}

Now solve for VO2,k+1VO_{2,k+1}:

VO2,k+1=VO2,k+tk+1−tkτ(Pk⋅ϵ+VO2RES−VO2,k)VO_{2,k+1} = VO_{2,k} + \frac{t_{k+1} - t_k}{\tau} \left( P_k \cdot \epsilon + VO_{2RES} - VO_{2,k} \right)

This recursive equation lets you compute oxygen uptake step by step from power data, and so simulate the transient response of oxygen consumption during exercise. It's especially useful for interval training or transitions between work rates. Once you have implemented it in a spreadsheet, you can visualize how VO₂ responds to different power profiles. In the chart below, pick an exercise protocol and change τ and ε to see how the VO₂ response follows its steady-state target.

Resting VO₂, VO₂RES (mL/min): 350

Select Exercise Protocol:

VO₂ Kinetics

Power Profile

Notice that VO₂ always lags behind its steady-state target: the larger τ, the longer the lag.

Comes back inW′ balance · Interactive Wbal simulation

Supercompensation: the Banister model

The supercompensation principle describes how the body adapts to training stimuli over time.
After a training session, performance capacity temporarily decreases due to fatigue.
During recovery, capacity returns to its baseline and may then exceed it. This is supercompensation: an improved capacity compared to pre-training levels.

If a new training stimulus is applied at the right time, the adaptation curve can progressively shift upward, leading to long-term performance improvements.
Conversely, if training is too frequent (insufficient recovery) or too infrequent (loss of adaptation), performance gains are reduced or absent.

This is a dynamic problem by nature: how you perform today depends on all the training of the previous weeks, not only on today's session. Mathematically, it is often described using the Banister model (also known as the fitness-fatigue model).
In this model:

  • Fitness represents the positive adaptations to training.
  • Fatigue represents the negative, transient effects of training load.
  • Performance is modeled as fitness minus fatigue, each with its own time constant and sensitivity.

For a given training load (TL), we generically refer to CTL as the chronic training load (or fitness, i.e., the positive and long-lasting effects of training) and to ATL as the acute training load (or fatigue, i.e., the short-term effects of training that are negative for performance). With TSB we refer to the training stress balance (or form, i.e., the balance between the positive and negative effects of training on the performance potential).

The model is two first-order equations, one for fitness and one for fatigue, combined into form:

dCTL(t)dt=TL(t)−CTL(t)τ1\frac{dCTL(t)}{dt} = TL(t) - \frac{CTL(t)}{\tau_1}
dATL(t)dt=TL(t)−ATL(t)τ2\frac{dATL(t)}{dt} = TL(t) - \frac{ATL(t)}{\tau_2}
TSB(t)=k1⋅CTL(t)−k2⋅ATL(t)TSB(t)= k_1\cdot CTL(t) - k_2\cdot ATL(t)

Each day's load adds to both fitness and fatigue, and both then decay, fitness slowly (time constant τ1\tau_1) and fatigue quickly (τ2\tau_2). This is Banister's original formulation, and it is the one the chart below uses. Training software such as TrainingPeaks uses a variant, τ1dCTL(t)dt+CTL(t)=TL(t)\tau_1 \frac{dCTL(t)}{dt} + CTL(t) = TL(t), in which CTL and ATL become moving averages of the load and TSB is simply CTL − ATL.

To replicate it in your spreadsheet, use the same finite-difference recipe as for VO₂, one day at a time:

CTLk+1=CTLk+(tk+1−tk)(TLk−CTLkτ1)CTL_{k+1} = CTL_k + (t_{k+1} - t_k)\left(TL_k - \frac{CTL_k}{\tau_1}\right)

and the same for ATL with τ2\tau_2. This way you can feed the model any sequence of training loads.

The parameters τ1\tau_1, τ2\tau_2, k1k_1, and k2k_2 are estimated empirically and depend on the athlete and on the software: TrainingPeaks, for example, uses τ1=42\tau_1=42 days and τ2=7\tau_2=7 days.

In the chart below, set the daily training load with the sliders at the bottom and see how CTL, ATL, and TSB respond. The chart uses τ1=27\tau_1=27 days, τ2=10\tau_2=10 days, k1=1k_1=1, and k2=2k_2=2.

📈 Training Load & Form

📊 Daily Training Load


🔨 Set Daily Training Load

01
02
03
04
05
06
07
08
09
010
011
012
013
014
015
016
017
018
019
020
021
022
023
024
025
026
027
028
029
030

Comes back inModule 6 · 🔍 More examples of spurious correlations

One equation, many systems

Completely different systems often follow identical mathematical rules. The first-order differential equation of this module appears throughout science and engineering:

τddty(t)+y(t)=x(t)\tau \frac{d}{dt} y(t) + y(t) = x(t)

Here are four applications of this same equation, each from a different field:

1. VO₂ kinetics (exercise physiology)

τddtVO2(t)+VO2(t)=VO2SS(t)\tau \frac{d}{dt} VO_2(t) + VO_2(t) = VO_{2SS}(t)

Here, oxygen uptake responds to changes in steady-state demand with time constant τ\tau.

2. Exponentially weighted moving average (data analysis)

τddtxˉ(t)+xˉ(t)=x(t)\tau \frac{d}{dt} \bar{x}(t) + \bar{x}(t) = x(t)

Here, a smoothed average xˉ(t)\bar{x}(t) tracks a noisy input signal x(t)x(t) with smoothing parameter τ\tau. This is the continuous-time version of the exponentially weighted moving average, used in fields from stock market analysis to GPS tracking.

3. Banister fitness-fatigue model (training science)

τ1ddtCTL(t)+CTL(t)=TL(t)\tau_1 \frac{d}{dt} CTL(t) + CTL(t) = TL(t)
τ2ddtATL(t)+ATL(t)=TL(t)\tau_2 \frac{d}{dt} ATL(t) + ATL(t) = TL(t)

Here, fitness (CTL) and fatigue (ATL) respond to the same training load with different time constants. This is the moving-average variant used by training software; Banister's original form, used in the chart above, multiplies the load by τ\tau, which changes the size of each response but not its shape.

4. RC circuit (electrical engineering)

RCddtVC(t)+VC(t)=Vin(t)RC \frac{d}{dt} V_C(t) + V_C(t) = V_{in}(t)

Here, capacitor voltage VC(t)V_C(t) responds to input voltage Vin(t)V_{in}(t) with time constant τ=RC\tau = RC (resistance × capacitance).

The universal pattern. In every case, we have:

  • A state variable that changes over time: VO2VO_2, moving average, fitness/fatigue, or voltage
  • An input that drives the change: power demand, data signal, training load, or applied voltage
  • A time constant that determines response speed: metabolic τ\tau, smoothing τ\tau, adaptation τ\tau, or electrical τ\tau

The mathematical structure is identical; only the physical interpretation changes. A sports scientist optimizing training periodization uses the same differential equation as an electrical engineer designing a filter circuit.

Why does this matter for sports scientists?

  • In the form used by training software, CTL and ATL are exponentially weighted moving averages of the training load, one with a long time constant and one with a short one: the same equation, read twice.
  • The recursive (step-by-step) form you wrote for VO₂ works unchanged for fitness, fatigue, a smoothed power signal, or the voltage of a capacitor: only the names and the value of τ change.
  • VO₂ smooths a changing power signal the way an RC circuit smooths a changing voltage: changes much faster than τ barely show up in the output. In the model, a sprint much shorter than τ hardly moves VO₂.
  • Recognizing the same equation in another field lets you borrow that field's methods and solutions.

If you can read τddty+y=x\tau\frac{d}{dt}y + y = x in one of these fields, you can read it in all of them.

Comes back inLoop gain and oscillatory breathing · What comes nextWorkout Reserve · The exponentially weighted mean (EWM) frameworkModule 6 · Looking back


Second-order dynamic systems

The mass-spring-damper system

One of the main reasons we introduce second-order dynamic models in sports science is that they derive naturally from Newton's second law of motion, the foundation of classical mechanics. You may already be familiar with the general form of this law:

m⋅a(t)=∑i=1nFi(t)m \cdot a(t) = \sum_{i=1}^n F_i(t)

This tells us that an object with mass mm will experience an acceleration a(t)a(t) that is proportional to the sum of the n external forces F acting on it. Since acceleration is the second derivative of position, we can rewrite this more formally as:

ddtx(t)=v(t)\frac{d}{dt}x(t) = v(t)
d2dt2x(t)=a(t)\frac{d^2}{dt^2}x(t) = a(t)

So Newton's law becomes:

m⋅d2dt2x(t)=∑i=1nFi(t)m \cdot \frac{d^2}{dt^2}x(t) = \sum_{i=1}^n F_i(t)

In many physical systems, especially those involving movement or impact, we observe three key mechanical properties: inertia, stiffness, and damping. These are particularly relevant in sports settings, for example in the dynamics of skis, rackets, bicycles, and other equipment.

Each of these effects can be captured by a parameter in the classic mass-spring-damper system:

  • Inertia is represented by the mass m.
  • Stiffness corresponds to a restoring force proportional to the displacement x.
  • Damping corresponds to a resistive force proportional to velocity v.

For a mass connected to the ground through a spring and a damper, a typical equation of motion is:

m⋅d2dt2x(t)+C⋅ddtx(t)+K⋅x(t)=F(t)m \cdot \frac{d^2}{dt^2}x(t) + C \cdot \frac{d}{dt}x(t) + K \cdot x(t) = F(t)

Here, F(t) is a single generic external force applied to the system.
This type of equation is useful for modeling human biomechanics (for example, in running) or, as in this case, a mountain bike suspension system. In that system, changes in ground elevation drive corresponding changes in the rider's vertical position.

Let's define:

  • x1x_1: the vertical position of the rider (the mass),
  • x2x_2: the position of the ground (known and imposed),
  • m: the mass of the rider,
  • K: stiffness of the suspension,
  • C: damping coefficient,
  • L₀: the equilibrium length of the suspension.

Modeling a mountain bike suspension

The forces acting on the rider include:

  • Spring force: depends on the displacement between rider and ground,
  • Damping force: depends on the relative velocity between rider and ground,
  • Gravity: acts downward on the rider's mass, with gg the gravitational acceleration (9.81 m/s²).

We write the equation of motion as follows, replacing the second derivative of position with the first derivative of velocity (velocity is the derivative of position, so d2dt2x1=ddtv1\frac{d^2}{dt^2}x_1 = \frac{d}{dt}v_1):

m⋅ddtv1(t)+C⋅(v1(t)−v2(t))+K⋅((x1(t)−L0)−x2(t))=−m⋅gm \cdot \frac{d}{dt}v_1(t) + C \cdot (v_1(t)-v_2(t)) + K\cdot ((x_1(t)-L_0)-x_2(t)) = -m\cdot g

Rearranging for acceleration:

ddtv1(t)=1m(C⋅(v2(t)−v1(t))+K⋅(x2(t)−(x1(t)−L0))−m⋅g)\frac{d}{dt}v_1(t) = \frac{1}{m}(C \cdot (v_2(t)-v_1(t)) + K \cdot (x_2(t) - (x_1(t)-L_0)) - m\cdot g)

To use this model in a spreadsheet or code, we discretize the time derivative with the same finite difference approximation we used for VO₂:

ddtv(t)≈vk+1−vktk+1−tk\frac{d}{dt}v(t) \approx \frac{v_{k+1} - v_k}{t_{k+1} - t_k}

Substituting this into the equation:

v1,k+1−v1,ktk+1−tk=1m(C⋅(v2,k−v1,k)+K⋅(x2,k−(x1,k−L0))−m⋅g)\frac{v_{1,k+1} - v_{1,k}}{t_{k+1} - t_k} = \frac{1}{m}(C \cdot (v_{2,k} - v_{1,k}) + K \cdot (x_{2,k} - (x_{1,k} - L_0)) - m \cdot g)

And solving for v1,k+1v_{1,k+1}:

v1,k+1=(C⋅(v2,k−v1,k)+K⋅(x2,k−(x1,k−L0))−mg)⋅Δtm+v1,kv_{1,k+1} = \left( C \cdot (v_{2,k} - v_{1,k}) + K \cdot \left(x_{2,k} - (x_{1,k} - L_0)\right) - mg \right) \cdot \frac{\Delta t}{m} + v_{1,k}

where Δt=tk+1−tk\Delta t = t_{k+1} - t_k is the time step.

Since x2x_2 is imposed by the terrain, we can compute v2v_2 directly from it. Then, starting from an initial condition for x1x_1 and v1v_1 at time step k = 0, we can simulate the system step by step.

Try it yourself!
Use the sliders to adjust the parameters and see how the model responds to changes in mass, stiffness, damping, or terrain input.

Parameter selection

Fixed: g = 9.81 m/s², Δt = 10 ms, L₀ = 0.8 m

Vertical displacements

Velocities


Equation of motion for locomotor sports

Equations of motion, the second-order differential equations introduced above, form the basis for modeling the dynamics of many locomotor sports such as cycling, cross-country skiing, and running. In Module 1 we used a simplified versionM1 that kept only gravity; here we add every force.

In cycling, the primary forces acting on the system (cyclist and bicycle) are:

  • FP(t)F_P(t): Propulsive force generated by the cyclist
  • FW(t)F_W(t): Aerodynamic drag
  • FR(t)F_R(t): Rolling resistance
  • FG(t)F_G(t): Gravitational force due to slope

The dynamic model becomes:

m⋅d2dt2x(t)=FP(t)−FW(t)−FR(t)−FG(t)m \cdot \frac{d^2}{dt^2}x(t) = F_P(t) - F_W(t) - F_R(t) - F_G(t)

In cycling, we usually monitor effort and performance through power, not force. Power is defined as the product of force and velocity:

P(t)=F(t)⋅v(t)P(t) = F(t) \cdot v(t)

Writing d2dt2x\frac{d^2}{dt^2}x as ddtv\frac{d}{dt}v (velocity is the derivative of position) and multiplying both sides of the equation of motion by v(t)v(t) gives:

v(t)⋅m⋅ddtv(t)=v(t)⋅FP(t)−v(t)⋅FW(t)−v(t)⋅FR(t)−v(t)⋅FG(t)v(t) \cdot m \cdot \frac{d}{dt}v(t) = v(t) \cdot F_P(t) - v(t) \cdot F_W(t) - v(t) \cdot F_R(t) - v(t) \cdot F_G(t)

This equation expresses how the power produced by the cyclist is used to accelerate the system and overcome resistive forces.

We now replace each force on the right-hand side with its physical expression:

v(t)⋅m⋅ddtv(t)=P(t)−v(t)⋅12ρ⋅CdA⋅v2(t)−v(t)⋅mgcos⁡(β)⋅Crr−v(t)⋅mgsin⁡(β)v(t) \cdot m \cdot \frac{d}{dt}v(t) = P(t) - v(t) \cdot \frac{1}{2} \rho \cdot C_dA \cdot v^2(t) - v(t) \cdot m g \cos(\beta) \cdot C_{rr} - v(t) \cdot m g \sin(\beta)

Where:

  • mm is the total mass of cyclist and bicycle
  • ρ\rho is the air density
  • CdAC_dA is the aerodynamic drag coefficient times frontal area
  • CrrC_{rr} is the coefficient of rolling resistance
  • β\beta is the road gradient (in radians)
  • gg is gravitational acceleration (9.81 m/s²)

To use this model in a spreadsheet or in code, we discretize the time derivative using a finite difference approximation:

ddtv(t)≈vk+1−vktk+1−tk\frac{d}{dt}v(t) \approx \frac{v_{k+1} - v_k}{t_{k+1} - t_k}

Substituting into the equation:

vk⋅m⋅vk+1−vktk+1−tk=Pk−vk⋅12ρ⋅CdA⋅vk2−vk⋅mgcos⁡(β)⋅Crr−vk⋅mgsin⁡(β)v_k \cdot m \cdot \frac{v_{k+1} - v_k}{t_{k+1} - t_k} = P_k - v_k \cdot \frac{1}{2} \rho \cdot C_dA \cdot v_k^2 - v_k \cdot m g \cos(\beta) \cdot C_{rr} - v_k \cdot m g \sin(\beta)

Solving for vk+1v_{k+1}:

vk+1=vk+tk+1−tkm⋅vk[Pk−12ρ⋅CdA⋅vk3−mgcos⁡(β)⋅Crr⋅vk−mgsin⁡(β)⋅vk]v_{k+1} = v_k + \frac{t_{k+1} - t_k}{m \cdot v_k} \left[ P_k - \frac{1}{2} \rho \cdot C_dA \cdot v_k^3 - m g \cos(\beta) \cdot C_{rr} \cdot v_k - m g \sin(\beta) \cdot v_k \right]

And of course we do the same to find the distance:

xk+1=xk+(tk+1−tk)⋅vkx_{k+1} = x_k + (t_{k+1} - t_k) \cdot v_k

These recursive equations allow you to compute the velocity and distance at the next time step, based on the power input and environmental parameters.

Because it uses power instead of force, this model can be fed directly with power data and solved step by step in a spreadsheet or in code. The chart below does exactly that: set the power output, mass, air density, CdA, Crr and road gradient, and see how velocity and distance evolve over time.

🚴 Bicycle Locomotion Simulation

Fixed: g = 9.81 m/s², initial distance = 0 m, simulation duration = 60 s, time step Δt = 0.1 s

Velocity vs. Time

Distance vs. Time

With constant power, velocity moves toward a steady value: the speed at which the cyclist's power exactly balances drag, rolling resistance, and gravity, so that ddtv=0\frac{d}{dt}v = 0.


Wrap-up

Well done, you have completed the fourth module of this course.