System of \(r\) reactions involving \(p\) different substances, that can be described by \(j=1,\ldots,r\) reaction equations:

\[ \begin{aligned} k_{1j}R_1+\ldots+k_{pj}R_p \stackrel{\theta_j}{\longrightarrow} s_{1j}R_1+\ldots+s_{pj}R_p. \label{eqn:simple-chemical-reaction} \end{aligned} \]

The terms on the left-hand side are the reactants, whereas the terms on the right-hand side are the products.

We will denote with \(Y\) the overall state of the system, in particular \[ Y_{ti} = \mbox{number of molecules of substrate $i$ at time $t$.}\]

So, if we start the system at time \(0\) with \(Y_0\) molecules and if after time \(t\), reaction \(j\) is the first reaction to be triggered, then the new state \(Y_t\) for the \(p\) substrates can be calculated as follows, \[ Y_{ti} = Y_{0i} + s_{ij} - k_{ij}. \] The difference \(v_{ij}= s_{ij} - k_{ij}\) is called the net-change. For a set of \(r\) reactions and \(p\) species, the molecular transfer from reactant to product species is according to a net change of \(\mathbf{V}=\mathbf{S}-\mathbf{K}\), where \(\mathbf{S}\) denotes the \(p\times r\) dimensional matrix of stoichiometric coefficients of the products, \(\mathbf{K}\) is the \(p\times r\) dimensional matrix of stoichiometric coefficients of the reactants, and \(\mathbf{V}\) is called the \(p\times r\) dimensional or matrix.

If each single combination has rate \(\theta_j\), then the overall rate for reaction \(j\) is: \[ h_j(\mathbf{Y}_0,\boldsymbol\theta) = \theta_j \prod_{i=1}^p {Y_{0i} \choose k_{ij}} , \label{eq:sderate} \] where \({x \choose y}=0\) for \(x<y\). The waiting time for reaction \(j\) will therefore be \[ T_j \sim \mbox{Exp}(h_j(\mathbf{Y}_{0},\boldsymbol\theta)) . \]

Transcription factors binding to DNA system

Two different proteins, \(P\) and \(Q\), can bind to DNA. Using the notation \(P\cdot DNA\) to denote \(P\) bound to DNA, etc, the chemical equations for the process are

\[\begin{aligned} P+DNA \quad &\stackrel{r_{1}}{\longrightarrow} P \cdot DNA \\ Q+DNA \quad & \stackrel{r_{2}}{\longrightarrow} Q \cdot DNA \\ P+Q \cdot DNA \quad & \stackrel{r_3}{\longrightarrow} P \cdot Q \cdot D N A \\ Q+P \cdot DNA \quad & \stackrel{r _4}{\longrightarrow} P \cdot Q \cdot DNA \\ P \cdot DNA \quad & \stackrel{r_5}{\longrightarrow} P+DNA \\ Q \cdot DNA \quad & \stackrel{r_6}{\longrightarrow} Q + DNA \\ P \cdot Q \cdot DNA \quad & \stackrel{r_{7}}{\longrightarrow} P + Q \cdot DNA \\ P \cdot Q \cdot DNA \quad & \stackrel{r_{8}}{\longrightarrow} Q + P \cdot DNA \end{aligned}\]

We consider the following ODE system

\[ \begin{aligned} r_{1}: & {\text{Vegetation}} & \rightarrow & {\text{Litter}}\\ r_{2}: & {\text{Vegetation}} & \rightarrow & {\text{Heterotrophs}}\\ r_{3}: & {\text{Heterotrophs}} & \rightarrow & {\text{Litter}}\\ r_{4}: & {\text{Soil}} & \rightarrow & {\text{Vegetation}} \\ r_{5}: & {\text{Litter}} & \rightarrow & {\text{Soil}}\\ r_{6}: & {\text{Soil}} & \rightarrow & {\text{Ground water}}\\ r_{7}: & {\text{Grounf water}} & \rightarrow & {\text{Surface water}}\\ r_{8}: & {\text{Surface water}} & \rightarrow & \emptyset \\ r_{9}: & {\text{Soil}} & \rightarrow & {\text{Surface water}}\\ r_{10}: & {\text{Soil}} & \rightarrow & \emptyset \end{aligned} \]

consider also the external influence

\[ \begin{aligned} q_{1}: & \emptyset & \rightarrow & {\text{Litter}}\\ q_{2}: & \emptyset & \rightarrow & {\text{Vegetation}}\\ q_{3}: & \emptyset & \rightarrow & {\text{Surface water}} \end{aligned} \]

The stoichiometric matrix is an important matrix to study the dynamics of a reaction system. For example, the mean change of the state \(\mathbf{Y}_t\) at time \(t\) is given by \[ E [\mathbf{Y}_{t+dt} - \mathbf{Y}_t|\mathbf{Y}_t] = \mathbf{Vh}(\mathbf{Y}_t,\boldsymbol\theta)dt. \label{eq:sdemean} \] Therefore, the system is in mean if \[ \mathbf{V h}(\mathbf{\mu},\boldsymbol\theta) = 0.\] It is possible to find the mean stationary states \(\mathbf{\mu}\) by solving the above equation. In particular, one needs to find the null space of \(\mathbf{V}\) and then match it with \(\mathbf{h}(\mathbf{\mu},\boldsymbol\theta)\).

We are interested in how a quasi-reaction system develops when started at \(\mathbf{Y}(0)\). Defining the mean function \(\boldsymbol{\mu_Y}(t) = E[\boldsymbol{Y}_t|\boldsymbol{Y}_0]\), it is clear that \(\boldsymbol{\mu_Y}(t)\) is the solution of the ordinary differential equation \[ \begin{align} \frac{d\boldsymbol{\mu_Y}(t)}{dt} &=& \bf{V h}(\mu_Y (t),\boldsymbol\theta) \nonumber \\ &=& \underbrace{{\bf V} \mbox{diag}\{\boldsymbol\theta\}}_{\boldsymbol{P_\theta}}\bf{H}(\mu_Y(t)). \end{align} \] Given the polynomial nature of \(\boldsymbol{H}\), where \[H_j(\boldsymbol{\mu_Y}) = \prod_{i=1}^p {\mu_{Yi} \choose k_{ij}}\]

We consider \(n\) observations at locations \(t_1\leq t_2\leq \ldots \leq t_n\) and that each observation is given as \[ \mathbf{Y}_{t_i} = \boldsymbol{\mu_Y}(t_i) + \epsilon_i, \] with \(i=1,\ldots,n\) and \(V(\epsilon_i)=\sigma^2\).

The idea is then to find the best possible \(\boldsymbol{\theta}\), so that one obtains \[ \begin{align} \frac{d\boldsymbol{\hat{\mu}_Y}(t)}{dt} &\approx& \mathbf{V h}(\boldsymbol{\hat{\mu}_Y}(t),\boldsymbol{\theta}) \nonumber \\ &=& \underbrace{\mathbf{V} \mbox{diag}\{\mathbf{H}(\boldsymbol{\hat\mu_Y}(t))\}}_{\mathbf{g}(t)} \boldsymbol{\theta} \end{align} \]

which is linear in \(\boldsymbol\theta\).

By integrating both sides, we obtain \[ \boldsymbol{\hat{\mu}_Y}(t) \approx \boldsymbol{\mu_0}+\underbrace{\int_{t_1}^t \mathbf{g}(s)ds}_{\mathbf{G}(t)} \boldsymbol\theta, \] where \(\boldsymbol{\mu_0}=\boldsymbol{\mu_Y}(t_1)\) is the initial value of the ODE. The best is defined as the one that minimizes the integrated difference between the left and right-hand sides. \[ (\boldsymbol{\hat{\theta}},{\boldsymbol{\hat\mu_0}}) = \arg\min_{\boldsymbol{\theta},\boldsymbol{\mu_0}} \int_{t_1}^{t_n} ||\boldsymbol{\hat{\mu}_Y}(t) - \boldsymbol{\mu_0} - \boldsymbol{G}(t)\boldsymbol{\theta}||^2dt \]

Thus,

\[ (\boldsymbol{\hat{\theta}},{\boldsymbol{\hat\mu_0}}) = \arg\min_{\boldsymbol{\theta},\boldsymbol{\mu_0}} \begin{bmatrix}\boldsymbol{\theta} \\ \boldsymbol{\mu_0}\end{bmatrix}^T \underbrace{\begin{bmatrix} \int_{t_1}^{t_n} \mathbf{G}^T(t)\mathbf{G}(t)dt & \int_{t_1}^{t_n} \mathbf{G}^T(t)dt \\ \int_{t_1}^{t_n} \mathbf{G}(t)dt & (t_n-t_1) \mathbf{I}_p \end{bmatrix}}_{\mathbf{D}} \begin{bmatrix} \boldsymbol{\theta} \\ \boldsymbol{\mu_0} \end{bmatrix} -2 \underbrace{\begin{bmatrix} \int_{t_1}^{t_n} \boldsymbol{\hat{\mu}_Y}(t)^T\mathbf{G}(t)dt \\ \int_{t_1}^{t_n} \boldsymbol{\hat{\mu}_Y}(t)^Tdt \end{bmatrix}}_{\mathbf{d}} \begin{bmatrix} \boldsymbol{\theta} \\ \boldsymbol{\mu_0} \end{bmatrix} \] In principle, one has, again, an explicit least squares estimate for \(\boldsymbol\theta\), \[ \begin{bmatrix} \boldsymbol{\hat\theta} \\ \boldsymbol{\hat\mu_0} \end{bmatrix} = \boldsymbol D^{-1}\boldsymbol d. \]