Indexing

The quantities of interest in the mathematical modeling of reaction networks can been indexed in a number of different ways. This is highlighted by the disparity between approaches which primarily emphasize the effect of reactions and those which emphasize the placement of complexes in the reaction graph. For example, notice that ambiguity arises when we attempt to track the complexes in the following reaction network:

${\begin{array}{rclccc}{\mathcal {A}}_{1}&\;\longrightarrow \;&{\mathcal {A}}_{2}&{\mathcal {A}}_{1}\;\;\rightleftarrows \;\;{\mathcal {A}}_{2}\\[0.1cm]{\mathcal {A}}_{2}&\;\longrightarrow \;&{\mathcal {A}}_{1}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\Longleftrightarrow \;\;\;\;\;\;\;\;\;\;&\nwarrow \;\;\;\;\;\swarrow \\[0.1cm]{\mathcal {A}}_{2}&\;\longrightarrow \;&{\mathcal {A}}_{3}+{\mathcal {A}}_{4}&{\mathcal {A}}_{3}+{\mathcal {A}}_{4}.\\[0.1cm]{\mathcal {A}}_{3}+{\mathcal {A}}_{4}&\;\longrightarrow \;&{\mathcal {A}}_{1}&&&\end{array}}$

On the left we have the reaction network represented as a list of elementary reactions whereas on the right we have the reaction graph of the network. Notice that, for instance, the complex corresponding to ${\vec {y}}=[1,0,0,0]^{T}$ can be indexed three times as the reactant complex of the first reaction (${\vec {y}}_{1}$) and the product of the second and fourth reaction (${\vec {y}}_{2}'$ and ${\vec {y}}_{4}'$) or only once as a fixed node (${\vec {y}}_{k}$) with multiple adjacent reactions (one, two, and four). In the second setting, we index the reactions according to the fixed nodes between which they interact so that the first reaction would be corresponded to $y_{k}\rightarrow y_{i}$ and the second and fourth reaction would be corresponded to $y_{j}\rightarrow y_{k}$ for some $i$ and $j$.

Reaction-centered indexing

In reaction-centered indexing, the reactions are numbered sequentially from $1,\ldots ,r,$ where $r$ is the number of reactions. The complexes are indexed according to the reaction they appear in and whether they are reactant or product complex (product complexes are indicated with a prime, so that the product complex of the $i^{th}$ complex is denoted ${\mathcal {C}}_{i}'$). For mass-action systems, the rate constant associated with each reaction is assigned the index corresponding to the reaction, that is to say, we set $k_{i}>0$, $i=1,\ldots ,r$.

Using the reaction-centered indexing scheme, we represent a chemical reaction network as

${\mathcal {R}}_{i}:\;\;\;\;\;\sum _{j=1}^{n}y_{ij}{\mathcal {A}}_{j}\;\;{\stackrel {k_{i}}{\longrightarrow }}\;\;\sum _{j=1}^{n}y_{ij}'{\mathcal {A}}_{j},\;\;\;\;\;{\mbox{for }}i=1,\ldots ,r,$

where $y_{ij}$ is the stoichiometric coefficient of the $j^{th}$ species in the reactant complexes of the $i^{th}$ reaction, and $y_{ij}'$ is the stoichiometric coefficient of the $j^{th}$ species in the product complex of the $i^{th}$ reaction. The complex vectors are given by ${\vec {y}}_{i}=[y_{i1},y_{i2},\ldots ,y_{in}]^{T}$ and ${\vec {y}}_{i}'=[y_{i1}',y_{i2}',\ldots ,y_{in}']^{T}$ for $i=1,\ldots ,r$.

For example, in the above example, we have

${\vec {y}}_{1}=\left[{\begin{array}{c}1\\0\\0\\0\end{array}}\right],\;\;\;{\vec {y}}_{1}'=\left[{\begin{array}{c}0\\1\\0\\0\end{array}}\right],\;\;\;{\vec {y}}_{2}=\left[{\begin{array}{c}0\\1\\0\\0\end{array}}\right],\;\;\;{\vec {y}}_{2}'=\left[{\begin{array}{c}1\\0\\0\\0\end{array}}\right],\;\;\;{\vec {y}}_{3}=\left[{\begin{array}{c}0\\1\\0\\0\end{array}}\right],\;\;\;{\vec {y}}_{3}'=\left[{\begin{array}{c}0\\0\\1\\1\end{array}}\right],\;\;\;{\vec {y}}_{4}=\left[{\begin{array}{c}0\\0\\1\\1\end{array}}\right],\;\;\;{\vec {y}}_{4}'=\left[{\begin{array}{c}1\\0\\0\\0\end{array}}\right].$

Complex-centered indexing

In complex-centered indexing, the complexes are indexed according to their stoichiometrically distinctiveness from $i=1,\ldots ,m$. That is to say, each stoichiometrically distinct complex is assigned only one index, and one corresponding complex vector, regardless of whether it appears in just one or multiple reactions. The reactions are then indexed as ordered pairs depending on which stoichiometric distinct complexes the reaction connections and the direction of the connection. For mass-action systems, the rate constant associated with each reaction is assigned the indices corresponding to the complexes the reaction connects, that is to say, we set $k(i,j)>0$, $i,j=1,\ldots ,m$, if $({\mathcal {C}}_{i},{\mathcal {C}}_{j})\in {\mathcal {R}}$.

Using the complex-centered indexing scheme, we represent a chemical reaction network as

$\sum _{k=1}^{n}y_{ik}{\mathcal {A}}_{j}\;\;{\stackrel {k(i,j)}{\longrightarrow }}\;\;\sum _{k=1}^{n}y_{jk}{\mathcal {A}}_{j},\;\;\;\;\;{\mbox{for }}({\mathcal {C}}_{i},{\mathcal {C}}_{j})\in {\mathcal {R}}$

where $y_{ik}$ is the stoichiometric coefficient of the $k^{th}$ species in the $i^{th}$ complex. The complex vectors are given by ${\vec {y}}_{i}=[y_{i1},y_{i2},\ldots ,y_{in}]^{T}$ for $i=1,\ldots ,m$.

For example, in the above example, if we set ${\mathcal {C}}_{1}={\mathcal {A}}_{1}$, ${\mathcal {C}}_{2}={\mathcal {A}}_{2}$, and ${\mathcal {C}}_{3}={\mathcal {A}}_{3}+{\mathcal {A}}_{4}$, we have

${\vec {y}}_{1}=\left[{\begin{array}{c}1\\0\\0\\0\end{array}}\right],\;\;\;\;\;{\vec {y}}_{2}=\left[{\begin{array}{c}0\\1\\0\\0\end{array}}\right],\;\;\;\;\;{\vec {y}}_{3}=\left[{\begin{array}{c}0\\0\\1\\1\end{array}}\right].$

Conversion from reaction- to complex- centered indexing in mass-action systems

The dynamics of a mass-action system can be formulated in numerous equivalent forms. A primary difference between the representations is a result of the indexing scheme used and, in particular, the use of reaction- versus complex-centered indexing.

In the reaction-centered indexing scheme, it is most natural to represent a mass-action system as

${\frac {d{\vec {c}}}{dt}}=\Gamma {\vec {R}}({\vec {c}})$

where $\Gamma \in \mathbb {Z} ^{n\times r}$ is the stoichiometric matrix and ${\vec {R}}({\vec {c}})$ is the rate vector with entries $R_{i}({\vec {c}})=k_{i}c_{1}^{y_{i1}}\cdots c_{n}^{y_{in}}$.

In order to convert this expression into a form which utilizes complex-centered indexing, we re-index each rate constant $k_{\ell }$, $\ell =1,\ldots ,m$, to be $k(i,j)$ where the reaction from the $i^{th}$ complex ${\mathcal {C}}_{i}$ to the $j^{th}$ complex ${\mathcal {C}}_{j}$ is the $\ell ^{th}$ reaction in the reaction-centered indexing scheme. We now introduce two new matrices, $I_{a}\in \mathbb {Z} ^{r\times m}$ and $I_{k}\in \mathbb {R} _{\geq 0}^{m\times r}$ with entries defined as follows:

$[I_{a}]_{i\ell }=\left\{{\begin{array}{ll}1,\;\;\;\;\;&{\mbox{if }}{\mathcal {C}}_{i}{\mbox{ is the product complex of the }}\ell ^{th}{\mbox{ reaction}}\\-1,\;\;\;\;\;&{\mbox{if }}{\mathcal {C}}_{i}{\mbox{ is the reactant complex of the }}\ell ^{th}{\mbox{ reaction}}\\0,&{\mbox{otherwise}}\end{array}}\right.$

and

$[I_{k}]_{\ell i}=\left\{{\begin{array}{ll}k_{\ell },\;\;\;\;\;&{\mbox{if }}{\mathcal {C}}_{i}{\mbox{ is the reactant complex of the }}\ell ^{th}{\mbox{ reaction}}\\0,&{\mbox{otherwise}}\end{array}}\right.$

for $\ell =1,\ldots ,r$. Recalling the complex matrix $Y\in \mathbb {Z} _{\geq 0}^{n\times m}$ and the mass action vector $\Psi ({\vec {c}})\in \mathbb {R} _{\geq 0}^{m}$, it follows that the matrices $\Gamma$ and ${\vec {R}}({\vec {c}})$ decompose according to $\Gamma =YI_{a}$ and ${\vec {R}}({\vec {c}})=I_{k}\Psi ({\vec {c}})$ so that we can re-write the mass-action system as

${\frac {d{\vec {c}}}{dt}}=YI_{a}I_{k}\Psi ({\vec {c}}).$

In order to further simplify the equations, we notice that the kinetic or Kirchhoff matrix is given by the product of the matrices $I_{a}$ and $I_{k}$ under the re-indexing of the rate constants from the reaction-centered $k_{\ell }$ to the complex-centered $k(i,j)$, i.e. $A_{k}=I_{a}I_{k}\in \mathbb {R} ^{m\times m}$. It follows that

${\frac {d{\vec {c}}}{dt}}=YA_{k}\Psi ({\vec {c}}).$

Example

For the above example, we can use the reaction-centered indexing to write the mass-action system as

${\frac {d{\vec {c}}}{dt}}=\Gamma {\vec {R}}({\vec {c}})=\left[{\begin{array}{cccc}-1&1&0&1\\1&-1&-1&0\\0&0&1&-1\\0&0&1&-1\end{array}}\right]\left[{\begin{array}{c}k_{1}c_{1}\\k_{2}c_{2}\\k_{3}c_{2}\\k_{4}c_{3}c_{4}\end{array}}\right].$

Using the decomposition $\Gamma =YI_{a}$ and ${\vec {R}}({\vec {c}})=I_{k}\Psi ({\vec {c}})$ as above, we have that

${\frac {d{\vec {c}}}{dt}}=YI_{a}I_{k}\Psi ({\vec {c}})=\left[{\begin{array}{ccc}1&0&0\\0&1&0\\0&0&1\\0&0&1\end{array}}\right]\left[{\begin{array}{cccc}-1&1&0&1\\1&-1&-1&0\\0&0&1&-1\end{array}}\right]\left[{\begin{array}{ccc}k_{1}&0&0\\0&k_{2}&0\\0&k_{3}&0\\0&0&k_{4}\end{array}}\right]\left[{\begin{array}{ccc}c_{1}\\c_{2}\\c_{3}c_{4}.\end{array}}\right].$

This can be put into the form

${\frac {d{\vec {c}}}{dt}}=\left[{\begin{array}{ccc}1&0&0\\0&1&0\\0&0&1\\0&0&1\end{array}}\right]\left[{\begin{array}{ccc}-k_{1}&k_{2}&k_{4}\\k_{1}&-k_{2}-k_{3}&0\\0&k_{3}&-k_{4}\end{array}}\right]\left[{\begin{array}{ccc}c_{1}\\c_{2}\\c_{3}c_{4}.\end{array}}\right]$

which after exchanging indices to the complex-centered form yields

${\frac {d{\vec {c}}}{dt}}=YA_{k}\Psi ({\vec {c}})=\left[{\begin{array}{ccc}1&0&0\\0&1&0\\0&0&1\\0&0&1\end{array}}\right]\left[{\begin{array}{ccc}-k(1,2)&k(2,1)&k(3,1)\\k(1,2)&-k(2,1)-k(2,3)&0\\0&k(2,3)&-k(3,1)\end{array}}\right]\left[{\begin{array}{ccc}c_{1}\\c_{2}\\c_{3}c_{4}.\end{array}}\right]$