Extending invariants in four-variable Michaelis-Menten
Raw Michaelis-Menten Equation Analysis
We perform an example analysis of the Michael-Mentis equations [SS89]:
First we initiate the system from LaTeX and find the maximal scaling matrix such that the system is invariant.
Note that negative subscripts are turned into ‘m’ so that they comply with sympy. The ‘m’s are turned back into
negatives when printint to LaTeX using desr.tex_tools().
>>> import sympy
>>> from desr.matrix_normal_forms import smf
>>> from desr.ode_system import ODESystem
>>> from desr.ode_translation import ODETranslation, scale_action
>>> from desr.tex_tools import expr_to_tex
>>> system_tex = '''\frac{dE}{dt} &= - k_1 E S + k_{-1} C + k_2 C \\\\
... \frac{dS}{dt} &= - k_1 E S + k_{-1} C \\\\
... \frac{dC}{dt} &= k_1 E S - k_{-1} C - k_2 C \\\\
... \frac{dP}{dt} &= k_2 C'''
>>> original_system = ODESystem.from_tex(system_tex)
>>> original_system.variables
(t, C, E, P, S, k_1, k_2, k_m1)
>>> max_scal1 = ODETranslation.from_ode_system(original_system, naming_scheme=('t',['C','E','P','S'], 'c'))
>>> print('Variable order: ', max_scal1.variables_domain)
Variable order: (t, C, E, P, S, k_1, k_2, k_m1)
>>> print(f'Scaling Matrix:\n{max_scal1.scaling_matrix.__repr__()}')
Scaling Matrix:
Matrix([
[1, 0, 0, 0, 0, -1, -1, -1],
[0, 1, 1, 1, 1, -1, 0, 0]])
Inspect the scaling invariants:
>>> print('Invariants: ', max_scal1.invariants())
Invariants: Matrix([[k_m1*t, C*k_1/k_m1, E*k_1/k_m1, P*k_1/k_m1, S*k_1/k_m1, k_2/k_m1]])
Finding the reduced system is also easy.
Since the Hermite multiplier and inverse are compatible with the simplest parameter reduction scheme,
translate() will automatically perform this reduction.
>>> print(f'Reduced system:\n{max_scal1.translate(original_system)}')
Reduced system:
dt/dt = 1
dC/dt = -C*c0 - C + E*S
dE/dt = C*c0 + C - E*S
dP/dt = C*c0
dS/dt = C - E*S
dc0/dt = 0
Changing the variable order
In our previous example, we had \(k_{-1}\) at the end of the variable order, so that the algorithm tries to normalise using \(k_{-1}\).
Instead, we can choose to normalise by \(k_2\), by swapping around the last two variables.
Note that we need to recalculate the ODETranslation instance.
>>> original_system_reorder = original_system.copy()
>>> variable_order = list(original_system.variables)
>>> variable_order[-1], variable_order[-2] = variable_order[-2], variable_order[-1] # Swap the last two variables
>>> original_system_reorder.reorder_variables(variable_order)
>>> original_system_reorder.variables
(t, C, E, P, S, k_1, k_m1, k_2)
>>> max_scal1_reorder = ODETranslation.from_ode_system(original_system_reorder, naming_scheme=('t',['C','E','P','S'], 'c'))
>>> print('Invariants:', ', '.join(map(str, max_scal1_reorder.invariants())))
Invariants: k_2*t, C*k_1/k_2, E*k_1/k_2, P*k_1/k_2, S*k_1/k_2, k_m1/k_2
Now we can reduce to find another, equivalent system.
>>> reduced_system = max_scal1_reorder.translate(original_system_reorder)
>>> reduced_system
dt/dt = 1
dC/dt = -C*c0 - C + E*S
dE/dt = C*c0 + C - E*S
dP/dt = C
dS/dt = C*c0 - E*S
dc0/dt = 0
Extending a choice of invariants
We return to our original variable order: \(t, C, E, P, S, k_1, k_2, k_{-1}\).
Suppose we wish to study the invariants \(\frac{k_1}{k_{2}}C\) and \(\frac{k_1}{k_{-1}}P\). Then we must create a matrix representing these invariants:
We can easily check we have correct matrix:
>>> invariant_choice = sympy.Matrix([[0, 1, 0, 0, 0, 1, -1, 0],
... [0, 0, 0, 1, 0, 1, 0, -1]]).T
>>> scale_action(max_scal1.variables_domain, invariant_choice)
Matrix([[C*k_1/k_2, P*k_1/k_m1]])
Finding a maximal scaling matrix that can be used to rewrite the system in terms of these invariants is also simple.
>>> max_scal2 = max_scal1.extend_from_invariants(invariant_choice=invariant_choice)
>>> max_scal2
A=
Matrix([
[1, 0, 0, 0, 0, -1, -1, -1],
[0, 1, 1, 1, 1, -1, 0, 0]])
V=
Matrix([
[ 0, 0, 0, 0, 1, 0, 0, 0],
[ 0, 0, 1, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 1, 0, 0],
[ 0, 0, 0, 1, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 1, 0],
[ 0, -1, 1, 1, 0, 1, 1, 0],
[ 0, 0, -1, 0, 0, 0, 0, 1],
[-1, 1, 0, -1, 1, -1, -1, -1]])
W=
Matrix([
[1, 0, 0, 0, 0, -1, -1, -1],
[0, 1, 1, 1, 1, -1, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 1, 0, 0, 0, 0],
[1, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 1, 0]])
For Python code that steps through this procedure, see desr.examples.michaelis_menten.
Now, this transformation doesn’t satisfy the conditions of the parameter reduction scheme, so if we try to reduce it
translate() will use the dependent reduction scheme implemented in
translate_dep_var().
>>> max_scal2.invariants()
Matrix([[C*k_1/k_2, P*k_1/k_m1, k_m1*t, E*k_1/k_m1, S*k_1/k_m1, k_2/k_m1]])
>>> max_scal2.translate(original_system)
dt/dt = 1
dx0/dt = 0
dx1/dt = 0
dy0/dt = y0*(-y2*y5 - y2 + y2*y3*y4/(y0*y5))/t
dy1/dt = y0*y2*y5**2/t
dy2/dt = y2/t
dy3/dt = y3*(y0*y2*y5**2/y3 + y0*y2*y5/y3 - y2*y4)/t
dy4/dt = y4*(y0*y2*y5/y4 - y2*y3)/t
dy5/dt = 0
Here, x0 and x1 are auxiliary variables, which can be fixed at any value at all.
(y0, y1, y2, y3, y4) = (C*k_1/k_2, P*k_1/k_m1, k_m1*t, E*k_1/k_m1, S*k_1/k_m1) are our new dependent invariants.
Finally, y5 = k_2/k_m1 is the single parameter of the reduced system.
However, we can see that after performing a permutation of the columns, we can satisfy the parameter reduction scheme. While this isn’t implemented yet, we can do it by hand for the moment. We must apply the cycle \(\begin{pmatrix}0 & 1 & 3 & 2\end{pmatrix}\) to the last \(n-r\) columns.
>>> max_scal3 = max_scal2.herm_mult_n
>>> max_scal3.col_swap(0, 1)
>>> max_scal3.col_swap(0, 3)
>>> max_scal3.col_swap(0, 2)
>>> print(f'Permuted Vn:\n{repr(max_scal3)}')
Permuted Vn:
Matrix([
[1, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0],
[0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 1, 0],
[0, 1, 1, 1, 1, 0],
[0, -1, 0, 0, 0, 1],
[1, 0, -1, -1, -1, -1]])
>>> max_scal3 = sympy.Matrix.hstack(max_scal1.herm_mult_i, max_scal3)
>>> max_scal3 = ODETranslation(max_scal1.scaling_matrix, hermite_multiplier=max_scal3, naming_scheme=('t',['C','E','P','S'], 'c'))
>>> print(max_scal3.translate(original_system).diff_subs({'c0':'alpha'}))
dt/dt = 1
dC/dt = -C*alpha - C + E*S/alpha
dE/dt = C*alpha**2 + C*alpha - E*S
dP/dt = C*alpha**2
dS/dt = C*alpha - E*S
dalpha/dt = 0
So we have found a third different reparametrization of the Michaelis-Menten equations.
Todo
Add a method to ODETranslation that will try and re-order the last \(n-r\) columns so
that the parameter reduction scheme can be applied.