Biochemical systems equation

The biochemical systems equation is a compact equation of nonlinear differential equations for describing a kinetic model for any network of coupled biochemical reactions and transport processes.{{cite journal |last1=Reder |first1=Christine |title=Metabolic control theory: A structural approach |journal=Journal of Theoretical Biology |date=November 1988 |volume=135 |issue=2 |pages=175–201 |doi=10.1016/S0022-5193(88)80073-0|pmid=3267767 |bibcode=1988JThBi.135..175R }}{{cite journal |last1=Hofmeyr |first1=Jan-hendrik S. |title=Metabolic control analysis, in a nutshell |journal=In Proceedings of the 2 Nd International Conference on Systems Biology |date=2001 |pages=291–300 |citeseerx=10.1.1.324.922 }}

The equation is expressed in the following form:

\dfrac{{\bf dx}}{dt} = {\bf N} {\bf v} ({\bf x} (p), p)

The notation for the dependent variable x varies among authors. For example, some authors use s, indicating species. x is used here to match the state space notation used in control theory but either notation is acceptable.

\bf N is the stoichiometry matrix which is an m by n matrix of stoichiometry coefficient. m is the number of species and n the number of biochemical reactions. The notation for \bf N is also variable. In constraint-based modeling the symbol \bf N tends to be used to indicate 'stoichiometry'. However in biochemical dynamic modeling{{cite journal |last1=Stucki |first1=Jörg W. |title=Stability analysis of biochemical systems— A practical guide |journal=Progress in Biophysics and Molecular Biology |date=1979 |volume=33 |issue=2 |pages=99–187 |doi=10.1016/0079-6107(79)90027-0|pmid=674688 }} and sensitivity analysis, \bf N tends to be in more common use to indicate 'number'. In the chemistry domain, the symbol used for the stoichiometry matrix is highly variable though the symbols S and N have been used in the past.{{cite journal |last1=Fjeld |first1=M. |last2=Asbjørnsen |first2=O.A. |last3=Åström |first3=K.J. |title=Reaction invariants and their importance in the analysis of eigenvectors, state observability and controllability of the continuous stirred tank reactor |journal=Chemical Engineering Science |date=September 1974 |volume=29 |issue=9 |pages=1917–1926 |doi=10.1016/0009-2509(74)85009-8|bibcode=1974ChEnS..29.1917F |url=https://lup.lub.lu.se/record/57576b1c-211d-4252-a5c7-9f2e9f9e03b0 }}{{cite journal |last1=Park |first1=David J. M. |title=SMISS, stoichiometric matrix inversion for steady state metabolic networks |journal=Computer Programs in Biomedicine |date=1 September 1975 |volume=5 |issue=1 |pages=46–60 |doi=10.1016/0010-468X(75)90026-4|pmid=1164840 }}

\bf v is an n-dimensional column vector of reaction rates, and p is a p-dimensional column vector of parameters.

Example

Given the biochemical network:

X_o \stackrel{v_1}\longrightarrow\ x_1 \stackrel{v_2}\longrightarrow\ x_2 \stackrel{v_3}\longrightarrow\ x_3 \stackrel{v_4}\longrightarrow\ X_1

where X_o and X_1 are fixed species to ensure the system is open. The system equation can be written as:{{cite journal |last1=Cornish-Bowden |first1=Athel |last2=Hofmeyr |first2=Jan-Hendrik S. |title=The Role of Stoichiometric Analysis in Studies of Metabolism: An Example |journal=Journal of Theoretical Biology |date=May 2002 |volume=216 |issue=2 |pages=179–191 |doi=10.1006/jtbi.2002.2547|pmid=12079370 |bibcode=2002JThBi.216..179C }}

:

\mathbf{N} = \begin{bmatrix}

1 & -1 & \phantom{+}0 & \phantom{+}0 \\

0 & \phantom{+}1 & -1 & \phantom{+}0 \\

0 & \phantom{+}0 & \phantom{+}1 & -1 \\

\end{bmatrix},\

\mathbf{v} = \begin{bmatrix}

v_1 \\

v_2 \\

v_3 \\

v_4 \\

\end{bmatrix}

So that:

\begin{bmatrix}

\dfrac{dx_1}{dt} \\[4pt]

\dfrac{dx_2}{dt} \\[4pt]

\dfrac{dx_3}{dt} \\[4pt]

\dfrac{dx_4}{dt} \\[4pt]

\end{bmatrix}

= \begin{bmatrix}

1 & -1 & \phantom{+}0 & \phantom{+}0 \\

0 & \phantom{+}1 & -1 & \phantom{+}0 \\

0 & \phantom{+}0 & \phantom{+}1 & -1 \\

\end{bmatrix}

\begin{bmatrix}

v_1 \\

v_2 \\

v_3 \\

v_4 \\

\end{bmatrix}

The elements of the rate vector will be rate equations that are functions of one or more species x_i and parameters, p. In the example, these might be simple mass-action rate laws such as v_2 = k_2 x_1 where k_2 is the rate constant parameter. The particular laws chosen will depend on the specific system under study. Assuming mass-action kinetics, the above equation can be written in complete form as:

\begin{bmatrix}

\dfrac{dx_1}{dt} \\[4pt]

\dfrac{dx_2}{dt} \\[4pt]

\dfrac{dx_3}{dt} \\[4pt]

\dfrac{dx_4}{dt} \\[4pt]

\end{bmatrix}

= \begin{bmatrix}

1 & -1 & \phantom{+}0 & \phantom{+}0 \\

0 & \phantom{+}1 & -1 & \phantom{+}0 \\

0 & \phantom{+}0 & \phantom{+}1 & -1 \\

\end{bmatrix}

\begin{bmatrix}

k_1 X_o \\

k_2 x_1 \\

k_3 x_2 \\

k_4 x_3 \\

\end{bmatrix}

Analysis

The system equation can be analyzed by looking at the linear response of the equation around the steady-state with respect to the parameter \bf p .{{cite book |last1=Heinrich |first1=Reinhart |last2=Schuster |first2=Stefan |title=The Regulation of Cellular Systems |date=1996 |publisher=Springer |location=New York, NY |isbn=0412032619}} At steady-state, the system equation is set to zero and given by:

0 = {\bf N} {\bf v} ({\bf x} ({\bf p}), {\bf p})

Differentiating the equation with respect to {\bf p} and rearranging gives:

\dfrac{d{\bf x}}{d{\bf p}} = -\left( {\bf N} \frac{\partial \mathbf{v}}{\partial \mathbf{x}}\right)^{-1} {\bf N} \frac{\partial \mathbf{v}}{\partial \mathbf{p}}

This derivation assumes that the stoichiometry matrix has full rank. If this is not the case, then the inverse won't exist.

=Example=

For example, consider the same problem from the previous section of a linear chain. The matrix \frac{\partial \mathbf{v}}{\partial \mathbf{x}} is the unscaled elasticity matrix:

:

\mathcal{E} =

\begin{bmatrix}

\dfrac{\partial v_1}{\partial x_1} & \cdots & \dfrac{\partial v_1}{\partial x_m} \\ \vdots & \ddots & \vdots \\

\dfrac{\partial v_n}{\partial x_1} & \cdots & \dfrac{\partial v_n}{\partial x_m}

\end{bmatrix}.

In this specific problem there are 3 species (m=3) and 4 reaction steps ( n = 4), the elasticity matrix is therefore a m \times n = 3\ \mbox{by}\ 4 matrix. However, a number of entries in the matrix will be zero. For example \partial v_1/\partial x_3 will be zero since x_3 has no effect on v_1 . The matrix, therefore, will contain the following entries:

:

\mathcal{E} =

\begin{bmatrix}

\dfrac{\partial v_1}{\partial x_1} & 0 & 0 \\

\dfrac{\partial v_2}{\partial x_1} & \dfrac{\partial v_2}{\partial x_2} & 0 \\

0 & \dfrac{\partial v_3}{\partial x_2} & \dfrac{\partial v_3}{\partial x_3} \\

0 & 0 & \dfrac{\partial v_4}{\partial x_3} \\

\end{bmatrix}.

The parameter matrix depends on which parameters are considered. In Metabolic control analysis, a common set of parameters are the enzyme activities. For the sake of argument, we can equate the rate constants with the enzyme activity parameters. We also assume that each enzyme, k_i, only can affect its own step and no other. The matrix \frac{\partial \mathbf{v}}{\partial \mathbf{p}} is the unscaled elasticity matrix with respect to the parameters. Since there are 4 reaction steps and 4 corresponding parameters, the matrix will be a 4 by 4 matrix. Since each parameter only affects one reaction, the matrix will be a diagonal matrix:

:

\mathcal{E} =

\begin{bmatrix}

\dfrac{\partial v_1}{\partial k_1} & 0 & 0 & 0 \\

0 & \dfrac{\partial v_2}{\partial k_2} & 0 & 0 \\

0 & 0 & \dfrac{\partial v_3}{\partial k_3} & 0 \\

0 & 0 & & \dfrac{\partial v_4}{\partial k_4} \\

\end{bmatrix}.

Since there are 3 species and 4 reactions, the resulting matrix \frac{d{\bf x}}{d {\bf p}} will be a 3 by 4 matrix

D = \mathcal{E}^{1}_{1} \mathcal{E}^{2}_{2} (\mathcal{E}^{3}_{3}-\mathcal{E}^{4}_{3})+\mathcal{E}^{1}_{1} \mathcal{E}^{3}_{2} \mathcal{E}^{4}_{3}-\mathcal{E}^{1}_{2} \mathcal{E}^{3}_{2} \mathcal{E}^{4}_{3}

\vphantom{ }

\frac{d{\bf x}}{d {\bf p}} = \frac{1}{D}

\left[

\begin{array}{ll}

\mathcal{E}^{1}_{k_1} (\mathcal{E}^{2}_{2} (\mathcal{E}^{3}_{3}-\mathcal{E}^{4}_{3})+\mathcal{E}^{3}_{2} \mathcal{E}^{4}_{3}) &

-\mathcal{E}^{3}_{2} \mathcal{E}^{4}_{3} \mathcal{E}^{2}_{k_2} \\

\mathcal{E}^{1}_{2} \mathcal{E}^{1}_{k_1} (\mathcal{E}^{3}_{3}-\mathcal{E}^{4}_{3}) &

\mathcal{E}^{1}_{1} \mathcal{E}^{2}_{k_2} (\mathcal{E}^{3}_{3}-\mathcal{E}^{4}_{3}) \\

\mathcal{E}^{1}_{2} \mathcal{E}^{3}_{2} \mathcal{E}^{1}_{k_1} &

\mathcal{E}^{1}_{1} \mathcal{E}^{3}_{2} \mathcal{E}^{2}_{k_2} \\

\end{array}

\right.

\qquad\qquad\qquad\quad

\left.

\begin{array}{ll}

\mathcal{E}^{2}_{2} \mathcal{E}^{4}_{3} \mathcal{E}^{3}_{k_3} &

\mathcal{E}^{2}_{2} \mathcal{E}^{3}_{3} \mathcal{E}^{4}_{k_4} \\

\mathcal{E}^{4}_{3} \mathcal{E}^{3}_{k_3} (\mathcal{E}^{1}_{1}-\mathcal{E}^{1}_{2})

&

\mathcal{E}^{3}_{3} \mathcal{E}^{4}_{k_4} (\mathcal{E}^{1}_{2}-\mathcal{E}^{1}_{1}) \\

\mathcal{E}^{1}_{1} \mathcal{E}^{2}_{2} \mathcal{E}^{3}_{k_3} &

-\mathcal{E}^{4}_{k_4} (\mathcal{E}^{1}_{1} (\mathcal{E}^{2}_{2}-\mathcal{E}^{3}_{2})+\mathcal{E}^{1}_{2} \mathcal{E}^{3}_{2}) \\

\end{array}

\right]

Each expression in the matrix describes how a given parameter influences the steady-state concentration of a given species. Note that this is the unscaled derivative. It is often the case that the derivative is scaled by the parameter and concentration to eliminate units as well as turn the measure into a relative change.

Assumptions

The biochemical systems equation makes two key assumptions:

  1. Species exist in a well-stirred reactor, so there are no spatial gradients.{{cite book |last1=Cowan |first1=Ann E. |last2=Moraru |first2=Ion I. |last3=Schaff |first3=James C. |last4=Slepchenko |first4=Boris M. |last5=Loew |first5=Leslie M. |title=Computational Methods in Cell Biology |chapter=Spatial Modeling of Cell Signaling Networks |date=2012 |volume=110 |pages=195–221 |doi=10.1016/B978-0-12-388403-9.00008-4|pmid=22482950 |pmc=3519356 |isbn=9780123884039 }}{{cite journal |last1=Fell |first1=David A. |title=Theoretical analyses of the functioning of the high- and low-Km cyclic nucleotide phosphodiesterases in the regulation of the concentration of adenosine 3′,5′-cyclic monophosphate in animal cells |journal=Journal of Theoretical Biology |date=May 1980 |volume=84 |issue=2 |pages=361–385 |doi=10.1016/S0022-5193(80)80011-7|pmid=6251314 |bibcode=1980JThBi..84..361F }}{{cite journal |last1=Kholodenko |first1=Boris N. |title=Cell-signalling dynamics in time and space |journal=Nature Reviews Molecular Cell Biology |date=March 2006 |volume=7 |issue=3 |pages=165–176 |doi=10.1038/nrm1838|pmid=16482094 |pmc=1679905 }}
  2. Species concentrations are high enough so that stochastic effects are negligible{{cite journal |last1=Gillespie |first1=Daniel T. |title=Exact stochastic simulation of coupled chemical reactions |journal=The Journal of Physical Chemistry |date=December 1977 |volume=81 |issue=25 |pages=2340–2361 |doi=10.1021/j100540a008|s2cid=2606191 }}{{cite journal |last1=Gillespie |first1=Daniel T. |title=Stochastic Simulation of Chemical Kinetics |journal=Annual Review of Physical Chemistry |date=1 May 2007 |volume=58 |issue=1 |pages=35–55 |doi=10.1146/annurev.physchem.58.032806.104637|pmid=17037977 |bibcode=2007ARPC...58...35G }}{{cite journal |last1=Andrews |first1=Steven S |last2=Bray |first2=Dennis |title=Stochastic simulation of chemical reactions with spatial resolution and single molecule detail |journal=Physical Biology |date=September 2004 |volume=1 |issue=3 |pages=137–151 |doi=10.1088/1478-3967/1/3/001|pmid=16204833 |bibcode=2004PhBio...1..137A |s2cid=16394428 }}

See also

References