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. The equation is expressed in the following form:
d x d t = N v ( x ( p ) , p ) {\displaystyle {\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.
N {\displaystyle {\bf {N}}} is the stoichiometry matrix which is an m {\displaystyle m} by n {\displaystyle n} matrix of stoichiometry coefficient. m {\displaystyle m} is the number of species and n {\displaystyle n} the number of biochemical reactions. The notation for N {\displaystyle {\bf {N}}} is also variable. In constraint-based modeling the symbol N {\displaystyle {\bf {N}}} tends to be used to indicate 'stoichiometry'. However in biochemical dynamic modeling and sensitivity analysis, N {\displaystyle {\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.
v {\displaystyle {\bf {v}}} is an n-dimensional column vector of reaction rates, and p {\displaystyle p} is a p-dimensional column vector of parameters.
Example Given the biochemical network:
X o ⟶ v 1 x 1 ⟶ v 2 x 2 ⟶ v 3 x 3 ⟶ v 4 X 1 {\displaystyle 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 {\displaystyle X_{o}} and X 1 {\displaystyle X_{1}} are fixed species to ensure the system is open. The system equation can be written as:
N = [ 1 − 1 + 0 + 0 0 + 1 − 1 + 0 0 + 0 + 1 − 1 ] , {\displaystyle \mathbf {N} ={\begin{bmatrix}1&-1&{\phantom {+}}0&{\phantom {+}}0\\0&{\phantom {+}}1&-1&{\phantom {+}}0\\0&{\phantom {+}}0&{\phantom {+}}1&-1\\\end{bmatrix}},\ } v = [ v 1 v 2 v 3 v 4 ] {\displaystyle \mathbf {v} ={\begin{bmatrix}v_{1}\\v_{2}\\v_{3}\\v_{4}\\\end{bmatrix}}}
So that:
[ d x 1 d t d x 2 d t d x 3 d t d x 4 d t ] = [ 1 − 1 + 0 + 0 0 + 1 − 1 + 0 0 + 0 + 1 − 1 ] {\displaystyle {\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}}} [ v 1 v 2 v 3 v 4 ] {\displaystyle {\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 {\displaystyle x_{i}} and parameters, p. In the example, these might be simple mass-action rate laws such as v 2 = k 2 x 1 {\displaystyle v_{2}=k_{2}x_{1}} where k 2 {\displaystyle 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:
[ d x 1 d t d x 2 d t d x 3 d t d x 4 d t ] = [ 1 − 1 + 0 + 0 0 + 1 − 1 + 0 0 + 0 + 1 − 1 ] {\displaystyle {\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}}} [ k 1 X o k 2 x 1 k 3 x 2 k 4 x 3 ] {\displaystyle {\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 p {\displaystyle {\bf {p}}} . At steady-state, the system equation is set to zero and given by:
0 = N v ( x ( p ) , p ) {\displaystyle 0={\bf {N}}{\bf {v}}({\bf {x}}({\bf {p}}),{\bf {p}})}
Differentiating the equation with respect to p {\displaystyle {\bf {p}}} and rearranging gives:
d x d p = − ( N ∂ v ∂ x ) − 1 N ∂ v ∂ p {\displaystyle {\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 ∂ v ∂ x {\displaystyle {\frac {\partial \mathbf {v} }{\partial \mathbf {x} }}} is the unscaled elasticity matrix:
E = [ ∂ v 1 ∂ x 1 ⋯ ∂ v 1 ∂ x m ⋮ ⋱ ⋮ ∂ v n ∂ x 1 ⋯ ∂ v n ∂ x m ] . {\displaystyle {\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 {\displaystyle m=3} ) and 4 reaction steps ( n = 4 {\displaystyle n=4} ), the elasticity matrix is therefore a m × n = 3 by 4 {\displaystyle m\times n=3\ {\mbox{by}}\ 4} matrix. However, a number of entries in the matrix will be zero. For example ∂ v 1 / ∂ x 3 {\displaystyle \partial v_{1}/\partial x_{3}} will be zero since x 3 {\displaystyle x_{3}} has no effect on v 1 {\displaystyle v_{1}} . The matrix, therefore, will contain the following entries:
E = [ ∂ v 1 ∂ x 1 0 0 ∂ v 2 ∂ x 1 ∂ v 2 ∂ x 2 0 0 ∂ v 3 ∂ x 2 ∂ v 3 ∂ x 3 0 0 ∂ v 4 ∂ x 3 ] . {\displaystyle {\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 {\displaystyle k_{i}} , only can affect its own step and no other. The matrix ∂ v ∂ p {\displaystyle {\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:
E = [ ∂ v 1 ∂ k 1 0 0 0 0 ∂ v 2 ∂ k 2 0 0 0 0 ∂ v 3 ∂ k 3 0 0 0 ∂ v 4 ∂ k 4 ] . {\displaystyle {\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 d x d p {\displaystyle {\frac {d{\bf {x}}}{d{\bf {p}}}}} will be a 3 by 4 matrix
D = E 1 1 E 2 2 ( E 3 3 − E 3 4 ) + E 1 1 E 2 3 E 3 4 − E 2 1 E 2 3 E 3 4 {\displaystyle D={\mathcal {E}}_{1}^{1}{\mathcal {E}}_{2}^{2}({\mathcal {E}}_{3}^{3}-{\mathcal {E}}_{3}^{4})+{\mathcal {E}}_{1}^{1}{\mathcal {E}}_{2}^{3}{\mathcal {E}}_{3}^{4}-{\mathcal {E}}_{2}^{1}{\mathcal {E}}_{2}^{3}{\mathcal {E}}_{3}^{4}}
{\displaystyle {\vphantom {}}}
d