Simplification and construction of ODEs

The general reaction scheme is \begin{align} \color{blue} m + \color{red}s &\mathop{\rightleftharpoons}_{k_{-1,s}}^{k_{1,s}} c_{ms} \mathop{\rightarrow}^{k_{2,s}} (1 - p_s) \color{red} s \\ \color{red} s + \color{green} r &\mathop{\rightleftharpoons}_{k_{-1,r}}^{k_{1,r}} c_{sr} \mathop{\rightarrow}^{k_{2,r}} (1-p_r) \color{green} r \end{align} According to mass action kinetics the temporal behavior of the state variables $(m, s, r, c_{ms}, c_{sr})$ can be described by the following set of ordinary differential equations. Background degradation is assumed to be proportional to concentration and production rates are included as separate terms. \begin{align} \frac{dm}{dt} &= \alpha_m - \beta_m m - k_{1,s}ms + k_{-1,s} c_{ms} \\ \frac{ds}{dt} &= \alpha_s - \beta_s s - k_{1,s}ms + (k_{-1,s}+(1-p_s)k_{2,s}) c_{ms} - k_{1,r}sr + k_{-1,r} c_{sr} \\ \frac{dr}{dt} &= \alpha_r - \beta_r r - k_{1,r} sr + (k_{-1,r} + (1-p_r)k_{2,r}) c_{sr} \\ \frac{d c_{ms}}{dt} &= k_{1,s}ms -(k_{-1,s}+k_{2,s})c_{ms} \\ \frac{d c_{sr}}{dt} &= k_{1,r}sr - (k_{-1,r} + k_{2,r}) c_{sr} \end{align} Assuming $\frac{d c_{ms}}{dt} = 0$ and $\frac{d c_{sr}}{dt} = 0$ it follows that \begin{align} c_{ms} &= \frac{k_{1,s}ms}{k_{-1,s} + k_{2,s}} \\ c_{sr} &= \frac{k_{1,r}sr}{k_{-1,r}+k_{2,r}} \end{align} which is inserted into the ODEs to elliminate $c_{ms}$ and $c_{sr}$ \begin{align} \frac{dm}{dt} &= \alpha_m - \beta_m m - (k_{1,s} - \frac{k_{1,s}k_{-1,s}}{k_{-1,s} + k_{2,s}}) ms \\ \frac{ds}{dt} &= \alpha_s - \beta_s s - \bigg( k_{1,s} - \frac{k_{1,s}k_{-1,s}}{k_{-1,s} + k_{2,s}} - (1-p_s)\frac{k_{1,s}k_{2,s}}{k_{-1,s} + k_{2,s}} \bigg) ms - \bigg( k_{1,r} - \frac{k_{1,r}k_{-1,r}}{k_{-1,r}+k_{2,r}} \bigg) sr \\ \frac{dr}{dt} &= \alpha_r - \beta_r r - \bigg( k_{1,r} - \frac{k_{1,r} k_{-1,r}}{k_{-1,r} + k_{2,r}} - (1- p_r) \frac{k_{1,r} k_{2,r}}{k_{-1,r} + k_{2,r}} \bigg) sr \end{align} For both $ms$ and $sr$ \begin{align} k_1-\frac{k_1 k_{-1}}{k_{-1} + k_2} = \frac{k_1(k_{-1} + k_2)}{k_{-1} + k_2} - \frac{k_1 k_{-1}}{k_{-1} + k_2} = \frac{k_1 k_2}{k_{-1} + k_2} \end{align} and the equations simplify to \begin{align} \label{system} \frac{dm}{dt} &= \alpha_m - \beta_m m - k_s ms \\ \frac{ds}{dt} &= \alpha_s - \beta_s s - p_s k_s ms - k_r sr \\ \frac{dr}{dt} &= \alpha_r - \beta_r r - p_r k_r sr \end{align} where $k_s = \frac{k_{1,s} k_{2,s}}{k_{-1,s} + k_{2,s}}$ and $k_r = \frac{k_{1,r} k_{2,r}}{k_{-1,r} + k_{2,r}}$.

Steady-state analysis

To simplify the notation on steady-state analysis we define a vector of variables ${\bf x} = (m, s, r)$ and a hovering dot as ${\bf \dot{x}} = \frac{d \bf x}{dt}$. Thus the steady-state problem can be compactly specified as ${\bf \dot{x}} = 0$ with some solution $\bf x^*$ satisfying the condition. In this notation the system of differential equations are given by $\dot{\bf x} = {\bf f}({\bf x})$ where ${\bf f}({\bf x}) = (f_1({\bf x}), f_2({\bf x}), f_3({\bf x}))$. For Model 1 the steady state problem \begin{equation} {\bf \dot{x}} = \begin{bmatrix} \alpha_m - \beta_m m - k_s m s \\ \alpha_s - \beta_s - k_r s r \\ \alpha_r - \beta_r r \end{bmatrix} = 0 \end{equation}

is solvable by back substitution and

\begin{equation} {\bf x^*} = \begin{bmatrix} \frac{\alpha_m (\beta_r \beta_s + \alpha_r k_r)}{\beta_r \alpha_s k_s + \beta_m \beta_r \beta_s + \beta_m \alpha_r k_r} \\ \frac{\alpha_s}{\beta_s + \frac{k_r \alpha_r}{\beta_r}} \\ \frac{\alpha_r}{\beta_r} \end{bmatrix} \end{equation}

The Jacobian given by $J_{ij} = \frac{\partial f_i ({\bf x})}{\partial x_j}$ for the simplified system of equations is calculated \begin{equation} \renewcommand{\arraystretch}{1.4} \frac{\partial f_i ({\bf x})}{\partial x_j} = \begin{bmatrix} \frac{\partial f_1}{\partial m} & \frac{\partial f_1}{\partial s} & \frac{\partial f_1}{\partial r} \\ \frac{\partial f_2}{\partial m} & \frac{\partial f_2}{\partial s} & \frac{\partial f_2}{\partial r} \\ \frac{\partial f_3}{\partial m} & \frac{\partial f_3}{\partial s} & \frac{\partial f_3}{\partial r} \end{bmatrix} = \begin{bmatrix} - \beta_m - k_s s & - k_s m & 0 \\ - p_s k_s s & - p_s k_s m - k_r r & - k_r s \\ 0 & - p_r k_r r & - \beta_r - p_r k_r s \end{bmatrix} \end{equation} For Model 1 $p_s = 0$ and $p_r = 0$ the Jacobian is triangular and the eigenvalues are the diagonal entries $\lambda = (- \beta_m - k_s s, \ -k_r r, \ - \beta_r)$. Because all eigenvalues have negative real parts any steady-state is stable.

For model 2 $p_s = 0$ and $p_r = 1$ the eigenvalues of the Jacobian are more complicated \begin{equation} \lambda = \begin{bmatrix} \frac{1}{2} \bigg( -\beta_r - k_r s - k_r r \pm \sqrt{\beta_r^2 + k_r^2 s^2 + 2 k_r^2 r s + k_r^2 r^2 + 2 k_r \beta_r s - 2 k_r \beta_r r} \bigg) \\ - \beta_m - k_s s \end{bmatrix} \end{equation} but fortunately every eigenvalue has negative real parts if \begin{equation} \beta_r + k_r s + k_r r > \sqrt{\beta_r^2 + k_r^2 s^2 + 2 k_r^2 r s + k_r^2 r^2 + 2 k_r \beta_r s - 2 k_r \beta_r r} \end{equation} equivalent to \begin{equation} \beta_r^2 + k_r^2 s^2 + 2 k_r^2 r s + k_r^2 r^2 + 2 k_r \beta_r s + 2 k_r \beta_r r > \beta_r^2 + k_r^2 s^2 + 2 k_r^2 r s + k_r^2 r^2 + 2 k_r \beta_r s - 2 k_r \beta_r r \end{equation} or simply \begin{equation} k_r \beta_r r > - k_r \beta_r r \end{equation} which is true for any steady-state. In conclusion any analytical steady-state is stable.

Parameter estimation

First order background decay of RNA is modeled by assuming that the decay is a stochastic process. \begin{equation} \label{decay} \frac{dN}{dt} = - \beta N \end{equation} where N is the amount of RNA and $\beta$ is a rate constant. The expression is analytically solvable by separation of variables \begin{equation} \int \frac{1}{N} dN = - \beta \int dt \end{equation} equivalent to \begin{align} N(t) &= N_0 e^{-\beta t} \\ \beta &= \frac{ln 2}{t_{1/2}} \end{align} Biological measurements of degradation often assumes first order decay for an experimental setup. One way to experimentally determine the decay rate is to inhibit transcription at $t_0$ and measure the amount of RNA at different time points (Ross, 1995).

Background degradation half-lifes of both chiP and chiX is determined to be $\sim$ 27 $min^{-1}$ (Overgaard, 2009). Hence \begin{equation} \beta_s = 0.0257 min^{-1} \end{equation} Our m consists of chiP fused to lacZ. lacZ mRNA is reported to have a half-life of $1.7 min^{-1}$(Liang, 1999). The actual value of the fusion RNA depends on the particular set of cis-RNA stability determinants but might not be easily inferable from the primary structure. The rate limiting step in mRNA degradation is thought to be endonucleolytic cleavage in the 5' untranslated region (Bechhofer, 1993). In our chiP-lacZ RNA the 5' UTR consists of the 5' UTR of chiP and we might argue that the stability determinants of lacZ have no effect. \begin{equation} \beta_m = 0.0257 min ^{-1} \end{equation} A wide range of sRNA stabilities have been reported (Vogel, 2003) half-lifes ranging from 2-32 min. Based on this limited dataset we could restrict $\beta_r = (0.02 - 0.35) min^{-1}$.

Second order degradation is the coupled RNA degradation and is represented by. \begin{equation} \frac{dm}{dt} = - \beta_m m - k_s m s \end{equation} This representation assumes unchanged background degradation. Separation of variables leads to \begin{equation} \label{second_order} \int \frac{1}{m} dm = - \beta_m \int dt - k_s \int s dt \end{equation} Assuming that s is constant at the conditions at which the data is measured the expression simplifies to \begin{equation} k_s = \frac{1}{[s]} \bigg( \frac{ln 2}{t_{1/2}} - \beta_m \bigg) \end{equation} Overgaard et al measured the s coupled half-life of m to be $\sim 4 min^{-1}$ at $[s] = 180 nmol$, read by eye from fig 1B @MicM(pBAD-micM/pNDM-ybfM). Hence \begin{equation} k_s = 0.000820 \ (nmol \ min)^{-1} \end{equation}

Dimensionless analysis

For model 1 $p_s = 0$ the governing equations are \begin{align} \frac{dm}{dt} &= \alpha_m - \beta_m m - k_s ms \\ \frac{ds}{dt} &= \alpha_s - \beta_s s - k_r sr \\ \frac{dr}{dt} &= \alpha_r - \beta_r r - p_r k_r s r \end{align} We restate the equations into a dimensionless form by rescaling in effect introducing the new variables: $\tau = t \beta_m$, $M = \frac{\beta_m}{\alpha_m}m$, $S = \frac{\beta_s}{\alpha_m}s$, $R = \frac{\beta_r}{\alpha_m}r$. By substitution the equations are algebraically manipulated into \begin{align} \frac{dM}{d \tau} &= 1 - M - \frac{k_s \alpha_m}{\beta_m \beta_s} MS \\ \frac{dS}{d \tau} &= \frac{\beta_s}{\beta_m} \bigg(\frac{\alpha_s}{\alpha_m} - S - \frac{k_r \alpha_m}{\beta_s \beta_r} SR \bigg) \\ \frac{dR}{d \tau} &= \frac{\beta_r}{\beta_m} \bigg(\frac{\alpha_r}{\alpha_m} - R - p_r \frac{k_r \alpha_m}{\beta_s \beta_r} S R\bigg) \end{align} At first this seems like an unnecessary complicated way of rewriting the equations, but the dimensionless parameter groups $p = (\frac{k_s \alpha_m}{\beta_m \beta_s}, \frac{k_r \alpha_m}{\beta_s \beta_r}, \frac{\beta_s}{\beta_m}, \frac{\beta_r}{\beta_m}, \frac{\alpha_s}{\alpha_m}, \frac{\alpha_r}{\alpha_m} )$ has 6 effective parameters fully determining the dynamics and steady-state of the system whereas the original equations had 8 parameters. $p_r$ is asumed to be either one or zero.