Reaction network realizations of rational biochemical systems and their structural properties

Abstract

In this paper, a frequently used representation of mass-action type reaction networks is extended to a more general system class where the reaction rates are in rational function form. An algorithm is given to compute a possible reaction graph from the kinetic differential equations. However, this structure is generally non-unique, as it is illustrated through the phenomenon of dynamical equivalence, when different reaction network structures correspond to exactly the same dynamics. It is shown that under some technical assumptions, the so-called dense realization containing the maximal number of reactions, forms a super-structure in the sense that the reaction graph of any dynamically equivalent reaction network is the sub-graph of the dense realization. Additionally, optimization based methods are given to find dynamically equivalent realizations with preferred properties, such as dense realizations or sparse realizations. The introduced concepts are illustrated by examples.

Appendices

Appendix 1: Nomenclature

Notation	Description
$A_\mathrm {k}\in \mathbb {R}^{m\times m}$	Kirchhoff matrix containing the reaction rate coefficients
$\tilde{A}_\mathrm {k}\in \mathbb {R}^{m\times \kappa }$	Kirchhoff matrix containing the principal reaction rate coefficients
$\mathcal {C} = \{C_1,\dots C_m\}$	Set of complexes
$d_1, \dots d_m$	Number of different kinetics from complexes $C_1,\dots C_m$, respectively
$D:x\in \overline{\mathbb {R}}^n_+ \rightarrow {\mathbb {R}}_+$	Positive multivariate polynomial with leading 1
$\varphi :x\in \overline{\mathbb {R}}^n_+ \rightarrow \overline{\mathbb {R}}^\kappa _+$	Vector of kinetics
$g_{il}$	The $l$th kinetics of the reaction in complex $C_i$
$k_{ijl}$	Reaction rate coefficient of the reaction $C_i \rightarrow C_j$ of the $l$th kinetics
$\kappa$	Number of different kinetics in the network
$E_d$	Edges of the weighted, directed reaction graph
$\eta ^{(i)}$	Complex composition vector of complex $C_i$
$m$	Number of complexes
$M = YA_\mathrm {k}$	Monomial coefficient matrix
$\tilde{M} = Y\tilde{A}_\mathrm {k}$	Kinetic coefficient matrix
$\mu _{ij}$	Stoichiometric coefficient of specie $X_i$ in (product) complex $C_j$
$n$	Number of species
$\nu _{ij}$	Stoichiometric coefficient of specie $X_i$ in (source) complex $C_j$
$P(x) \in \mathbb {R}^{\kappa \times m}$	Rate weighting matrix
$\Psi (x) \in \overline{\mathbb {R}}^n_+ \rightarrow \overline{\mathbb {R}}^m_+$	Monomial vector function
$r(x)\in \overline{\mathbb {R}}^n_+ \rightarrow \overline{\mathbb {R}}_+$	Reaction rate function
$\mathbb {R}^n$	The space of n-dimensional real vectors
$\overline{\mathbb {R}}^n_+ = [0,\infty )^n$	n-dimensional non-negative orthant
${\mathbb {R}}^n_+ = (0,\infty )^n$	n-dimensional positive orthant
$\mathbb {R}^{m\times n}$	The space of m $\times$ n dimensional real matrices
$\mathcal {R}$	Set of reactions
$\rho ^{(l,k)} = \eta ^{(l)}-\eta ^{(k)}$	Reaction vector corresponding to the reaction $C_k \rightarrow C_l$
$\mathcal {S}$	Set of species
$S$	Stoichiometric subspace
$V_d$	Vertices of the weighted, directed reaction graph
$X_1,\,\dots X_n$	Species
$x\in \mathbb {R}^n$	Species concentration
$Y$	Complex composition (stoichiometric) matrix
$\tilde{Y}$	Truncated complex composition (stoichiometric) matrix

Appendix 2: Further examples and decompositions of biochemical reaction rate functions

In this section we show for a set of biochemical reaction rate functions—mostly from the well-known modeling software COPASI [16], how to formulate them using elementary biochemical reaction rate functions according to (4). We have to note that not every reaction rate function can be formulated as (4).

Whenever a reaction rate function is reversible, two irreversible elementary reaction rate functions are formulated, $r_f$ for the forward and $r_r$ for the reverse reaction. In what follows $S$, $P$, $I$ and $A$ stand for the concentration of substrates, products, inhibitors and activators, respectively. The notations used for the constant parameters of the reaction rate functions are adapted from the COPASI software.

• Mass Action (reversible)

  $$\begin{aligned} k_1\prod _iS_j - k_2\prod _jP_j ~~\Longrightarrow ~~ r_f = k_1\prod _iS_j; ~~r_r = k_2\prod _iP_j, \end{aligned}$$

• Michaelis-Menten (reversible)

  $$\begin{aligned} \frac{V_f\frac{S}{K_{ms}}-V_r\frac{P}{K_{mp}} }{1 +\frac{S}{K_{ms}}+\frac{P}{K_{mp}}} ~~\Longrightarrow ~~ r_f = k_1\frac{S}{1 +k_3S+k_4P}; ~~ r_b = k_2\frac{P}{1 +k_3S+k_4P} \end{aligned}$$

  where $k_1 =\frac{V_f}{K_{ms}}$, $k_2=\frac{V_r}{K_{mp}}$, $k_3 = 1/K_{ms}$ and $k_4 = 1/K_{mp}$.

• Hill Cooperativity (irreversible)

  $$\begin{aligned} \frac{V S^h}{K^h + S^h} = \frac{V \left( \frac{S}{K}\right) ^h}{ 1 + \left( \frac{S}{K}\right) ^h} ~~\Longrightarrow ~~ r = k_1\frac{S^h}{1 +k_2S^h}, \end{aligned}$$

  where $k_1 = V/K^h$ and $k_2 = 1/K^h$.

• Ordered Bi Uni

  $$\begin{aligned}&\frac{V_f\left( S_aS_b-\frac{P}{K_{eq}}\right) }{ S_aS_b + K_{ma}S_b +K_{mb}S_a + \frac{V_f}{V_r K_{eq}}\left( K_{mp} + P\left( 1+\frac{S_a}{K_{ia}}\right) \right) }\Longrightarrow ~~ \\&r_f = k_f\frac{S_aS_b}{ 1+ k_1S_aS_b + k_3S_b +k_4S_a + k_5P + k_6PS_a};\\&r_b = k_r\frac{P}{ 1+ k_1S_aS_b + k_3S_b +k_4S_a + k_5P + k_6PS_a}, \end{aligned}$$

  where $k_1 = \frac{V_rK_{eq}}{V_fK_{mp}}$, $k_f = V_fk_1$, $k_r = V_fk_1/K_{eq}$, $k_3 = K_{ma}k_1$, $k_4 = K_{mb}k_1$, $k_5 = k_1/K_{mp}$ and $k_6 = k_1/(K_{mp}K_{ia})$.

• Allosteric inhibition (reversible)

  $$\begin{aligned}&\frac{V_f\frac{S}{K_{ms}}-V_r\frac{P}{K_{mp}}}{1 + \frac{S}{K_{ms}} + \frac{P}{K_{mp}} + \left( \frac{I}{K_{i}}\right) ^n} \Longrightarrow ~~ \\&r_f = k_f \frac{S}{1+k_1S + k_2P+k_3I^n} \\&r_r = k_r \frac{P}{1+k_1S + k_2P+k_3I^n}, \end{aligned}$$

  where $k_f = V_f/K_{ms}$, $k_r = V_r/K_{mp}$, $k_1 = 1/K_{ms}$, $k_2 = 1/K_{mp}$ and $k_3 = 1/K_{i}^n$.

• Mixed Inhibition (irreversible)

  $$\begin{aligned} \frac{V\frac{S}{K_{m}}}{1+\frac{I}{k_{is}}+ \frac{S}{K_{m}} + \frac{S}{K_{m}}\frac{I}{k_{ic}}}\Longrightarrow ~~ r = k_1 \frac{S}{1+k_2I + k_3S + k_4SI} , \end{aligned}$$

  where $k_1 = V/K_m$, $k_2 = 1/k_{is}$, $k_3 = 1/K_m$ and $k_4 = 1/(K_mk_{ic})$.

• Catalytic Activation (irreversible)

  $$\begin{aligned} \frac{V_{max}\frac{S}{K_{mS}} \frac{A}{K_{mA}} }{1 + \frac{S}{K_{mS}} + \frac{A}{K_{mA}}+\frac{S}{K_{mS}}\frac{A}{K_{mA}}} \Longrightarrow ~~ r = k\frac{SA}{1+k_1S + k_2A + k_3SA}, \end{aligned}$$

  where $k = V_{max}/(K_{mS}K_{mA})$, $k_1 = 1/K_{mS}$, $k_2 = 1/K_{mA}$ and , $k_3 = 1/(K_{mS}K_{mA})$.

• Substrate inhibition (irreversible)

  $$\begin{aligned} \frac{V\frac{S}{K_{m}}}{1+\frac{S}{K_{m}}+\left( \frac{S}{K_{si}}\right) ^2}\Longrightarrow ~~ r = k \frac{S}{1+k_1S + k_2S^2}, \end{aligned}$$

  where $k = V/K_{m}$, $k_1 = 1/K_{m}$ and $k_2 = 1/K_{si}^2$.

• Substrate activation (irreversible)

  $$\begin{aligned} \frac{V\left( \frac{S}{K_{ms}}\right) ^2}{1+\frac{S}{K_{sc}}+\frac{S}{K_{sa}}+\left( \frac{S}{K_{sa}}\right) ^2}\Longrightarrow ~~ r = k \frac{S^2}{1+k_1S + k_2S^2}, \end{aligned}$$

  where $k = V/K_{ms}^2$, $k_1 = \frac{K_{sa} + K_{sc}}{K_{sa}K_{sc}}$ and $k_2 = 1/K_{sa}^2$.

Appendix 3: Non-negativity of the solutions of bio-CRNs

In order to prove the non-negativity of the solution, one need to check the Lipschitz condition and the essential non-negativity of the right hand side of (7). It is easy to see, that the right hand sides of the ODEs are continuously differentiable, therefore they are locally Lipschitz.

To show the essential non-negativity of the right hand side functions, insert (11) into (7), for $p = 1,\,\dots n$, then the $p$-th equation reads as

$$\begin{aligned} f_p(x) = \sum _{l = 1}^m Y_{pl} \cdot \sum _{j = 1}^\kappa \tilde{A}_{\mathrm {k},lj}\cdot \varphi _j(x) \end{aligned}$$

(34)

Rewriting the sum over all the $\kappa$ kinetics into two sums: over the reactant complexes and over the kinetics in each of these complexes, one arrives to

$$\begin{aligned} f_p(x) = \sum _{l = 1}^m Y_{pl} \sum _{j = 1}^m \sum _{i=1}^{d_i} \tilde{A}_{\mathrm {k},l\,z_j+i}\cdot \varphi _{z_j+i}(x), \end{aligned}$$

(35)

where $z_j = \sum _{k=1}^{j-1}d_k$. Using (11) and (4)

$$\begin{aligned} f_p(x) = \sum _{l = 1}^m Y_{pl} \sum _{j = 1}^m \sum _{i=1}^{d_i} \tilde{A}_{\mathrm {k},l\,z_j+i}\cdot g_{ji}(x) = \sum _{l = 1}^m Y_{pl} \sum _{j = 1}^m \sum _{i=1}^{d_i} \tilde{A}_{\mathrm {k},l\,z_j+i}\cdot \frac{\prod _{o=1}^{n}x_o^{Y_{oj}} }{ D_{ji}(x) }. \end{aligned}$$

(36)

From the definition of $\tilde{A}_k$ we know that the coefficients $\tilde{A}_{\mathrm {k},l\,z_j+i}$ are negative when $l=j$ and non-negative otherwise. So we decompose the summation over $j$ into the two cases

$$\begin{aligned} f_p(x)= & {} \sum _{l = 1}^m Y_{pl} \sum _{j = 1, j\ne l }^m \sum _{i=1}^{d_i} \tilde{A}_{\mathrm {k},l\,z_j+i}\cdot \frac{\prod _{o=1}^{n}x_o^{Y_o,j} }{ D_{ji}(x) }\nonumber \\&- \sum _{l = 1}^m Y_{pl} \sum _{i=1}^{d_i} |\tilde{A}_{\mathrm {k},l\,z_j+i}|\cdot \frac{\prod _{o=1}^{n}x_o^{Y_{ol}} }{ D_{li}(x) }. \end{aligned}$$

(37)

Notice that the first term is always non-negative and the second term contains the factor $Y_{pl}x_p^{Y_{pl}}$. If $Y_{pl}=0$, then $\lim _{x_p\rightarrow 0}\frac{0\cdot x_p^{0}}{D(x)} = 0$. If $Y_{pl}>0$, since the denominator term (6) cannot approach zero, $\lim _{x_p\rightarrow 0}\frac{Y_{pl}x_p^{Y_{pl}}}{D(x)} = 0$ and $f_p$ is indeed essentially non-negative.

Appendix 4: Kinetic realizability

1.1 Dynamic equations with reaction vector formalism

The balance equations can also be written using the reaction rate functions and the reaction vectors. Let us denote the reaction rate corresponding to the reaction $(C_i,C_j,G_l)\in \mathcal R$ by $r_{ijl}$ and the set of all reaction rates by $\nabla = \{r_{ijl}\}$ for $i=1,\dots m; \,i=1,\dots m; \,l=1,\dots d_i$. Then, the balance equation reads as

(38)

where the reaction vector $\eta ^i$ is the $i$-th column of the complex composition matrix $Y$. Let $\{\omega _i\}$ be the standard basis on $\mathbb {R}^m$, thus $\eta ^{(i)} = Y\omega _i$. Inserting above yields

Rewrite the sum for each source complex, product complex and each kinetics

$$\begin{aligned} \dot{x} = Y\sum _{i = 1}^m \sum _{j = 1}^m\sum _{l=1}^{d_i} r_{ijl}(x)\left( \omega _j - \omega _i\right) . \end{aligned}$$

Here we used the convention, that if the reaction $(C_i,C_j,G_{l})\not \in \mathcal {R}$, then the corresponding principal reaction rate coefficient $k_{ijl}=0$, thus the reaction rate $r_{ijl}(x) = 0$. Let also use the definition of the bio-chemical reaction rate functions from (4)

$$\begin{aligned} \dot{x} = Y\sum _{i = 1}^m\sum _{j = 1}^m \sum _{l=1}^{d_i} k_{ijl}\frac{x^{\eta ^{(i)}}}{D_{il}(x)} \left( \omega _j - \omega _i\right) \end{aligned}$$

(39)

1.2 The proof of kinetic realizability

Consider the following autonomous ordinary differential equations

$$\begin{aligned} \dot{x}(t) = f(x(t)),\quad x(0)=x_0,\quad t \in [0, t_f), \end{aligned}$$

where $f(x):\mathbb {R}^{n}\rightarrow \mathbb {R}^{n}$ and $x \in \mathbb {R}^n$. Assume that the right hand side function $f$ is composed by the linear combination of biochemical reaction rate functions (4). Then, there exists a bio-CRN with $n$ species included if and only if for each $o = 1,\,\dots n$ the $f_o([x_1, x_2,\dots , x_{o-1},0,x_{o+1} \dots x_n])$ is a non-negative linear combination of the biochemical reaction rate functions.

Proof

(Sufficiency) In the $i$th equation each term in $f_i(x)$ has the form

$$\begin{aligned} a_i\frac{x_1^{p_1}x_2^{p_2}\dots x_i^{p_i}\dots x_n^{p_n}}{D(x)}, \end{aligned}$$

(40)

where the exponents $p_j\ge 0$ for $j = 1\dots n$ and specially $p_i > 0$, or

$$\begin{aligned} b_i\frac{x_1^{q_1}x_2^{q_2}\dots x_{i-1}^{q_{i-i}}x_{i+1}^{q_{i+i}}\dots x_n^{q_n}}{D(x)}, \end{aligned}$$

(41)

where the $b_i>0$ and the $x_i$ is missing from the nominator, i.e. $q_i \equiv 0$. Let

$$\begin{aligned} \hbox {sign}(x) = {\left\{ \begin{array}{ll} 1 &{} \hbox { if } x > 0\\ 0 &{} \hbox { if } x = 0\\ -1 &{} \hbox { if } x < 0 \end{array}\right. }. \end{aligned}$$

The following reaction

results the term (40) in $f_i$, but does not contribute to any other $f_j$, $j = 1, \dots n$, $j\ne i$, since the stoichiometric coefficients for any other $X_j$ are the same in the reactant and in the product side. Further, the reaction

contributes (41) in $f_i$, but nothing for the other $f_j$.

Since all terms in $f$ can be realized by a biochemical reaction with either the biochemical rate $|a_i|\dfrac{\sum _{j=1}^n x_j^{p_j}}{D(x)}$ or $b_i\dfrac{\sum _{j=1}^n x_j^{q_j}}{D(x)}$, we proved the sufficiency.

Necessity Let $o = 1,\,\dots n$ and recall (38)

Inserting the reaction rate functions (4) yields

Let $x_o = 0$. If $\eta ^{(i)}_o>0$, then the monomial $x^{\eta ^{(i)}}$ is zero and therefore, the corresponding terms in the sum disappear. On the other hand, if $\eta ^{(i)}_o = 0$, then the $0^0$ is the problem, which is treated as

$$\begin{aligned} \frac{k_{ijl}}{D_{il}(x)}x^{\eta ^{(i)}}\left( \eta ^{(j)}_o - \eta ^{(i)}_o\right)&= \lim _{x_o \rightarrow 0^+} \frac{k_{ijl}}{D_{il}(x)} x^{\eta ^{(i)}_1}x^{\eta ^{(i)}_2}\dots x^{\eta ^{(i)}_{o-1}} x_o^0 x^{\eta ^{(i)}_{o+1}} \dots x^{\eta ^{(i)}_{n}}\eta ^{(j)}_o \\&= \frac{k_{ijl}}{D_{il}([x_1, \dots , x_{o-1},0, x_{o-1}\dots x_n ])} x^{\eta ^{(i)}_1}x^{\eta ^{(i)}_2}\\&\qquad \dots x^{\eta ^{(i)}_{o-1}} 1 x^{\eta ^{(i)}_{o+1}} \dots x^{\eta ^{(i)}_{n}} \eta ^{(j)}_o. \end{aligned}$$

Let $\overline{ \mathcal {R}} = \{\hbox { terms in which }\eta ^{(i)}_o = 0\}$ and so

$$\begin{aligned} f_o([x_1, x_2,\dots , x_{o-1},0,x_{o+1} \dots x_n])=\sum _{\overline{\mathcal {R}}} \frac{\eta ^{(j)}_ok_{ijl}}{D_{il}(x)} x^{\eta ^{(i)}_1}x^{\eta ^{(i)}_2}\dots x^{\eta ^{(i)}_{o-1}} x^{\eta ^{(i)}_{o+1}} \dots x^{\eta ^{(i)}_{n}}, \end{aligned}$$

which is indeed a non-negative linear combination of elementary biochemical reaction rate functions. $\square$

Appendix 5: Properties of the realizations

In this section we prove the three claims from Sect. 4.3.

Indirect proof of P1 and P2 Let $M = Y\tilde{A}_k$ and assume that $(Y,\tilde{A}_\mathrm {k},P)$ is a dense realization of $\varSigma$ and thus according to (13), $\tilde{A}_\mathrm {k}$ has the most number of positive entries among the possible solutions of $M = Y\tilde{A}_\mathrm {k}$. Further assume that $\tilde{A}_\mathrm {k}'$ is also a valid modified Kirchhoff matrix solution of $M = Y\tilde{A}_\mathrm {k}'$, but there is $i,j,l$, $i \ne j$ such that $\tilde{A}_{\mathrm {k},i,z_j+l} = 0$ but $\tilde{A}_{\mathrm {k},i,z_j+l}' > 0$. Then, it follows from (23) that $(Y,\tilde{A}_\mathrm {k}'',P)$ with $\tilde{A}_\mathrm {k}'' = \frac{1}{2}\tilde{A}_\mathrm {k}' + \frac{1}{2}\tilde{A}_\mathrm {k}$ is also a valid dynamically equivalent realization of $\varSigma$, but $\tilde{A}_\mathrm {k}''$ has more positive elements than $\tilde{A}_\mathrm {k}$, which is a contradiction.

Proof of P3 ($\Rightarrow$) If the sparse and the dense realizations are structurally identical, then all the realizations are structurally the same, since the dense realisation is structurally unique. Thus any realisation is structurally unique. ($\Leftarrow$) If the structure of the realization is unique then the dense and the sparse realizations are trivially identical.