Derivatives at Tiny States: Michaelis-Menten

The problem in one sentence

Two ways of writing the same rate law can give the same value but a different derivative in floating point: the one with a state in the denominator of a denominator breaks down when that state gets tiny.

Michaelis-Menten

The Michaelis-Menten rate of an enzyme with substrate concentration S:

\[v = \frac{v_{max} \, S}{K_m + S}\]

Dividing numerator and denominator by S gives a form often found in model code:

\[v = \frac{v_{max}}{1 + K_m / S}\]

Mathematically both are the same function, and its derivative

\[\frac{dv}{dS} = \frac{v_{max} \, K_m}{(K_m + S)^2}\]

is finite for every S >= 0.

What goes wrong at tiny S

In ODE and DAE models concentrations can become very small, e.g. S = 1e-200, or exactly 0. For the second form:

• value: Km / S = 5e199, so v = vmax / (1 + 5e199), a tiny finite number. Nothing unusual.
• derivative (forward or reverse automatic differentiation applies the chain rule to every operation):
  • d(Km / S) = -Km / S^2, and S^2 = 1e-400 underflows to 0, giving -Inf
  • d(1 / D) = -dD / D^2, and D^2 = 2.5e399 overflows to Inf, giving Inf / Inf = NaN

At S = 0 exactly, Km / 0 = Inf, the value is still vmax / Inf = 0, and the derivative is again NaN.

Difference quotients do not see this, because they only use values. It shows up as soon as a Jacobian comes from automatic differentiation, e.g. a jac function for ida() / cvode() built with jacobian().

What ast2ast reports

f <- function(S) {
  argtypes(S |> type(vec(double)))
  rate <- function(s) {
    argtypes(s |> type(vec(double)))
    returntype(vec(double))
    vmax <- 2
    Km <- 0.5
    return(vmax / (1 + Km / s))
  }
  return(jacobian(rate, S))
}
fcpp <- ast2ast::translate(f)
fcpp(c(1e-200, 0.3))

Error: In 'return(vmax/(1 + Km/s))': jacobian: d f[1] / d x[1] is not finite,
although f[1] is. The first non-finite derivative appears in
'return(vmax/(1 + Km/s))'. Typical cause: a division by a (nearly) zero state,
e.g. vmax / (1 + K / c); the same value as vmax * c / (c + K) has a finite
derivative.

The check runs when a Jacobian contains a non-finite entry while the function value is finite: one extra sweep finds the first line where a derivative stops being finite. If the function value itself is NaN, the Jacobian is NaN too and no error is raised, as in R.

The fix: multiply through

Multiply numerator and denominator by the state that sits in the inner denominator:

return(vmax * s / (Km + s))

The value is the same, and S now only appears multiplied, so nothing becomes huge: at S = 1e-200 and at S = 0 the derivative is just vmax / Km.

     [,1]   [,2]
[1,]    4 0.0000
[2,]    0 1.5625

Rule of thumb

If a state sits in the denominator of a denominator, multiply numerator and denominator by that state.

With several substrates, multiply by all of them:

# vmax / (1 + K1 / a + K2 / b)
vmax * a * b / (a * b + K1 * b + K2 * a)

# vmax / (1 + K1 / a + K2 / b + K1 * K2 / (a * b))   (random-order bi-substrate)
vmax * a * b / (a * b + K1 * b + K2 * a + K1 * K2)

The multiplied-through form is usually the textbook form anyway; the form with the double fraction tends to come from code that was written for values only.