An algorithmic method to simplify the mathematical modeling of a reaction network using conservation laws

Question

Consider the reaction scheme:

$$\mathrm{S + E}\overset{k_1}{\to}\mathrm{C_1}, \qquad \mathrm{C_1}\overset{k_2}{\to}\mathrm{E + P}, \qquad \mathrm{S + C_1}\overset{k_3}{\underset{k_4}{\rightleftharpoons}}\mathrm{C_2}$$

where $\mathrm{S}$ is the substrate, $\mathrm{E}$ is the enzyme, $\mathrm{P}$ is the product, $\mathrm{C_1}$ and $\mathrm{C_2}$ are enzyme substrate complexes. Let $[\mathrm{S}] = s$, $[\mathrm{E}] = e$, $[\mathrm{C_1}] = c_1$, $[\mathrm{C_2}] = c_2$ and $[\mathrm{P}] = p$ be the concentrations of each respective chemical. I have used the law of mass reaction to convert this reaction into a system of differential equations:

$$ \left\{ \begin{align} \frac{ds}{dt} &= -k_1se-k_3sc_1+k_4c_2\\ \frac{de}{dt} &= -k_1se+k_2c_1\\ \frac{dc_1}{dt} &= k_1se-k_2c_1-k_3sc_1+k_4c_2\\ \frac{dc_2}{dt} &= k_3sc_1-k_4c_2\\ \frac{dp}{dt} &= k_2c_1 \end{align} \right. $$

I am asked to find a conservation equation and use it to simplify this system however I haven't been taught as to what a conservation equation is. Can anyone push me in the correct direction? (and check the above differential equations are correct).

Answer 1

Note: This answer is edited, don't forget to read the edit part as well.

---

I checked, your equations are correct for the system using mass-action kinetics.

Write a matrix, rows corresponded to equations and columns related to your complexes, $\{se,c1,ep,sc1,c2\}$, then find it's row kernel space, column kernel of its transpose. Choose a basis for that and write equations arrived by the inner product of them and the vector $[s,e,c1,c2,p]$. And replace them by one of the equations related to the species presented to the conservation equations. This way you reduce number of non-linear equations of your system and this is the reason it is helpful.

Doing this with Maple can be done as below;

with(LinearAlgebra):
M := Matrix([[-k[1],-k[3],k[4],0],[-k[1],0,0,k[2]],[k[1],-k[3],k[4],-k[2]],
[0,k[3],-k[4],0],[0,0,0,k[2]]]):
NullSpace(Transpose(M));

But if you pay attention you can avoid writing reaction constants as the result matrix has scalar times of columns of the previous one so it won't effect on its row kernel.

with(LinearAlgebra):
M := Matrix([[-1,-1,1,0],[-1,0,0,1],[1,-1,1,-1],[0,1,-1,0],[0,0,0,1]]):
NullSpace(Transpose(M));

Answer of both are same and is $$\{\begin{bmatrix} 0\\ 1\\ 1\\ 1\\ 0\end{bmatrix},\begin{bmatrix} 1\\ -2\\ -1\\ 0\\ 1\end{bmatrix}\}$$ Usually try to have positive combinations for your conservation laws so for example replace your second vector with sum of 2 times of the first with the second, the result is still a basis for row kernel.

So a maximal set of conservation laws for your reaction network is as below (this is not unique so I used maximal not maximum). $$\{e+c_1+c_2=T_1,s+c_1+2c_2+p=T_2\}$$ Where $T_1$ and $T_2$ are constants and you can find them using initial values of your variables in your differential system, it means $$\begin{array}{l} T_1=e(0)+c_1(0)+c_2(0)\\ T_2=s(0)+c_1(0)+2c_2(0)+p(0)\end{array}$$ You can replace those two conservation laws with equations of species participating in them so for example I replace the first one with $\dot{c_1}=k_1se-k_3sc_1+k_4c_2-k_2c_1$ and the second one with $\dot{s}=-k_1se-k_3sc_1+k_4c_2$. It's a good idea to replace them with those who are not linear already and contains more components. So the simplified differential equation is as below. $$\left\{\begin{array}{l} \dot{e}=-k_1se+k_2c_1\\ \dot{c_2}=k_3sc_1-k_4c_2\\ \dot{p}=k_2c_1\\ e+c_1+c_2=T_1\\ s+c_1+2c_2+p=T_2\end{array}\right.$$

---

Edit: In general we must consider the stoichiometric matrix, not the one introduced above with monomials of the reactants corresponding to the colums. It was ok for the network in this post, because each reaction has a seperate reactant complex. However, if we had something like $\require{mhchem} S+E\ce{<=>[k_1][k_2]}Y\ce{->[k_3]}S^*+E$, then we have terms such as $k_2y$ and $k_3y$ in our steady-state equations which can't be cancelled by real (non-parametric) linear summations even they have the same monomial. For a general case, the former answer must be corrected in this way, that the matrix that we want its row kernel is in fact the Stoichiometric matrix which is of the size number of species times number of reactions. The $(i,j)$-th entry contains the difference of the stochiometric coefficients of the $i$-th species in the $j$-th reaction. In your case, ordering species as $S,E,C_1,C_2,P$ and ordering the reactions the same way you indexed their rates, the stoichiometric matrix is the following. $$\begin{bmatrix}-1 & 0 & -1 & 1\\ -1 & 1 & 0 & 0\\ 1 & -1 & -1 & 1\\ 0 & 0 & 1 & 1\\ 0 & 1 & 0 & 0\end{bmatrix}$$ 4 and half year ago I gave the code to compute the row-kernel with Maple, so this time let's do it with Mathematica. The simple short code is the following.

matrix = {{-1, 0, -1, 1}, {-1, 1, 0, 0}, {1, -1, -1, 1}, {0, 0, 1, -1}, {0, 1, 0, 0}};
NullSpace[Transpose[matrix]]

To see things in a bit better display (one possible way to do it, you may find better ways);

matrix = {{-1, 0, -1, 1}, {-1, 1, 0, 0}, {1, -1, -1, 1}, {0, 0, 1, -1}, {0, 1, 0, 0}};
MatrixForm[matrix]
base = NullSpace[Transpose[matrix]];
baseNiceDisplay = {};
For[i = 1, i <= Length[base], i++,
  AppendTo[baseNiceDisplay, MatrixForm[base[[i]]]];
  ];
baseNiceDisplay

Anyway, you will get the following from Mathematica. $$\lbrace\begin{bmatrix}1\\ -2\\ -1\\ 0\\ 1\end{bmatrix},\begin{bmatrix}0\\ 1\\ 1\\ 1\\ 0\end{bmatrix} \rbrace$$ The same as the previous attempt. And the rest of the story is the same.