Consider the toy example set of reactions
\begin{align*}
S_0 &\xrightarrow{k_1} S_1 + S_2\\
S_2 &\xrightarrow{k_2}S_3 + S_4\\
S_1 + S_3 &\xrightarrow{k_3} S_5
\end{align*}
We have (noisy) data on the concentrations of the species as a function of time. We want to infer the rates $k_1$ to $k_3$.
Let’s write the derivatives:
\begin{align*}
\dot S_0 &=- k_1 S_0\\
\dot S_1 &= k_1 S_0 -k_3 S_1 S_3\\
\dot S_2 &= k_1 S_0 – k_2 S_2\\
\dot S_3 &= k_2 S_2 – k_3 S_1 S_3\\
\dot S_4 &= k_2 S_2\\
\dot S_5 &= k_3 S_1S_3
\end{align*}
In the general case, we will have a sequence of equations $$ L_i \xrightarrow{k_i} R_i.$$
A particular species diminishes whenever it’s on the left-hand side, by the product of the species on that side. Let’s define $$q_j = \prod_{i \in L_j} S_i.$$ If we let $A_{ij}$ be a binary variable that says if species $i$ is in the left-hand side of equation $j$, we get a contribution of $-A_{ij} q_j k_j .$
Similarly, the species will increase whenever it’s on the right-hand side, by the product of the species on the left hand side. If we let $B_{ij}$ be a binary variable that indicates if species $i$ is in the right-hand side of equation $j$, we get a contribution of $B_{ij} q_j k_j$. If we combine these, we get $$ \dot S_i = \sum_j (B_{ij} – A_{ij}) q_j k_j.$$ This is a linear system of equations in the unknown rates. If we define $ \mathbf Q(\mathbf s) = \diag{q_j}$, we can write the above as one matrix equation $$ \dot {\mathbf s} = (\mathbf B – \mathbf A) \mathbf Q(\mathbf s) \mathbf k = \mathbf H(\mathbf s) \mathbf k.$$
Inferring the rates
We have (noisy) measurements of the species concentrations at each point in time, so we have both $\dot {\mathbf s}$ and $\mathbf H (\mathbf s)$ at each time point. Therefore we can define a loss using the squared error at each time point, and add those up, with a regularizer on $\mathbf k$:
$$ L(\mathbf k) = {1 \over 2} \sum_t \|\dot{\mathbf s}_t – \mathbf{H}_t \mathbf{k} \|_2^2 + {1 \over 2} \alpha \|\mathbf k\|_2^2,$$ where $\mathbf H_t$ is shorthand for $\mathbf H(\mathbf s_t).$
The gradient is
$$ \nabla_{\mathbf k} L= \alpha \mathbf k-\sum_t \mathbf H_t^T(\mathbf s_t – \mathbf H_t \mathbf k) .$$Setting this to zero, we get
$$ (\alpha \mathbf I + \sum_t \mathbf{H}_t^T \mathbf{H}_t) \mathbf k = \sum_t \mathbf H_t^T \dot {\mathbf s_t}.$$ Defining $$ \mathbf G = (\alpha \mathbf I + \sum_t \mathbf{H}_t^T \mathbf{H}_t), \quad \mathbf y = \sum_t \mathbf H_t^T \dot {\mathbf s_t},$$ we have $$ \mathbf { G k} = \mathbf y.$$ Solving for $\mathbf k$ we get
$$ \mathbf k = \mathbf G^{-1} \mathbf y.$$
Trying it out
We try this out in https://github.com/stootoon/inferring-reaction-rates/blob/main/demo.ipynb, demonstrating that we can infer the correct rates in the noise-free setting, and can get close under moderate observation noise if we use regularization.
$$\blacksquare$$
Leave a Reply