In this post we show that, given a set of first order reactions with unknown rates, inferring the reaction rates from a dataset of instantaneous species concentrations is a linear regression problem. We solve it for some toy examples.<\/p>\n\n\n\n
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$. <\/p>\n\n\n\n
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*}<\/p>\n\n\n\n
In the general case, we will have a sequence of equations $$ L_i \\xrightarrow{k_i} R_i.$$<\/p>\n\n\n\n
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 .$<\/p>\n\n\n\n
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.$$<\/p>\n\n\n\n
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).$<\/p>\n\n\n\n
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.$$<\/p>\n\n\n\n