Michaelis-Menten with initial conditions and a constant

We derive the analysis of the Michaelis-Menten equations found in [SS89]: in a systematic manner.

Michaelis-Menten model, two variables

There are two chemical species \(E\) and \(S\), which reversibly combine to form \(ES\), which decomposes into \(E\) and \(P\).

\[ E+S\leftrightarrow ES \rightarrow E+P \]

After using the law of mass action, we are able to translate to the ODE system

\begin{align} \frac{ds}{dt} &= - k_1 e_0 s + k_1 c s + k_{-1} c \\ \frac{dc}{dt} &= k_1 e_0 s - k_1 c s - k_{-1} c - k_2 c \\ s(0) &= s_0. \end{align}

Model reduction in desr

First we must create the system in Python. We start with the same system, represented in a string in TeX.

>>> system_tex = '''\frac{ds}{dt} &= - k_1 e_0 s + k_1 c s + k_{-1} c \\\\
...          \frac{dc}{dt} &= k_1 e_0 s - k_1 c s - k_{-1} c - k_2 c'''
>>> system = ODESystem.from_tex(system_tex)

Add the constraint with the Michaelis-Menten constant \(K_m = k_2 + k_{-1}\).

>>> system.add_constraint('K_m', '(k_m1 + k_2) / k_1')

To match the analysis in Murray, We need to add in our initial condition for \(s\).

>>> system.update_initial_conditions({'s': 's_0'})
>>> system.exponent_matrix()
Matrix([
[1, 1, 1,  1, 1, 1,  1,  0,  0,  0],
[0, 0, 0, -1, 0, 1,  1,  0,  0,  0],
[1, 0, 0,  1, 0, 0, -1,  1,  0,  0],
[0, 0, 0,  1, 1, 0,  0,  0,  0,  0],
[1, 0, 0,  1, 1, 1,  0,  0,  1,  0],
[0, 1, 0,  0, 0, 0,  0,  0,  0,  1],
[0, 0, 1,  0, 0, 0,  1,  0, -1, -1],
[0, 0, 0,  0, 0, 0,  0,  0,  1,  0],
[0, 0, 0,  0, 0, 0,  0, -1,  0,  0]])

Choose the variable order.

>>> system.reorder_variables(['t', 's', 'c', 'K_m', 'k_m1', 'k_2', 'k_1', 'e_0', 's_0'])

Compute the reduced system by making an ODETranslation from the original system,

>>> translation = ODETranslation.from_ode_system(system, naming_scheme=('t',['s','c'], 'c'))

Optionally, observe some properties. First, the scaling matrix \(A\):

>>> translation.scaling_matrix
Matrix([
[1, 0, 0, 0, -1, -1, -1, 0, 0],
[0, 1, 1, 1,  0,  0, -1, 1, 1]])

The corresponding scaling invariants:

>>> translation.invariants()
Matrix([[k_1*s_0*t, s/s_0, c/s_0, K_m/s_0, k_m1/(k_1*s_0), k_2/(k_1*s_0), e_0/s_0]])

These don’t match those from Murray, but they are very close. We’ll recover it next.

Just for curiosity, here’s what we get if we translate the system at this point:

>>> translation.translate(system)
dt/dt = 1
dc/dt = -c*c1 - c*c2 - c*s + c3*s
ds/dt = c*c1 + c*s - c3*s
dc1/dt = 0
dc2/dt = 0
dc3/dt = 0
dc0/dt = 0
s(0) = 1
c0 == c1 + c2

Column operations to match the canonical system

Returning back to the original system with initial conditions,

>>> system
dt/dt = 1
ds/dt = c*k_1*s + c*k_m1 - e_0*k_1*s
dc/dt = -c*k_1*s - c*k_2 - c*k_m1 + e_0*k_1*s
dK_m/dt = 0
dk_m1/dt = 0
dk_2/dt = 0
dk_1/dt = 0
de_0/dt = 0
ds_0/dt = 0
s(0) = s_0
K_m == (k_2 + k_m1)/k_1
>>> translation = ODETranslation.from_ode_system(system, naming_scheme=('tau',['u','v'], 'c'))

Do some column operations on the Hermite multiplier

>>> translation.multiplier_add_columns(2, 8, 1)  # Scale time by e_0 not s_0
>>> translation.multiplier_add_columns(4, 8, -1)  # Scale c by e_0

Optionally, see the invariants and substitutions

>>> translation.invariants()
Matrix([[e_0*k_1*t, s/s_0, c/e_0, K_m/s_0, k_m1/(k_1*s_0), k_2/(k_1*s_0), e_0/s_0]])

We have substitutions to make. Note that many get normalized to \(1\).

>>> translation.translate_parameter_substitutions(system)
{t: tau/c3, s: u, c: c3*v, K_m: c0, k_m1: c1, k_2: c2, k_1: 1, e_0: c3, s_0: 1}

Translate

>>> reduced_system = translation.translate(system)
>>> reduced_system
dtau/dtau = 1
du/dtau = c1*v + u*v - u
dv/dtau = -c1*v/c3 - c2*v/c3 - u*v/c3 + u/c3
dc1/dtau = 0
dc2/dtau = 0
dc3/dtau = 0
dc0/dtau = 0
u(0) = 1
c0 == c1 + c2

Now we have a system, but the names of symbols don’t match. We’ll do substitutions.

>>> to_sub =  {'c0': 'K', 'c1': 'K-L', 'c2':'L', 'c3':'epsilon'}
>>> reduced_system.diff_subs(to_sub, subs_constraints=True)
dtau/dtau = 1
du/dtau = K*v - L*v + u*v - u
dv/dtau = -K*v/epsilon - u*v/epsilon + u/epsilon
dK/dtau = 0
dL/dtau = 0
depsilon/dtau = 0
u(0) = 1

The only step that desr cannot do is multiply \(\frac{dv}{d\tau}\) by \(\epsilon\). That has to be done by hand. And finally we have the result from [Mur07] (6.13).

\begin{align*} \frac{du}{d\tau} & = -u + (u+K-\lambda)v & u(0)&=1 \\ \epsilon \frac{dv}{d\tau} &= u - (u+K)v & v(0)&=0. \end{align*}