Exploiting fluctuations in gene expression to detect causal interactions between genes

Definition 1. A co-transition event is a transition event in the system where two or more components change simultaneously. That is, a reaction where the step size d has more than one vector component.

An example would be a conversion event, where a chemical species $z_m$ converts to another chemical species $z_n$

${\displaystyle {\displaystyle (z_m,z_n)\stackrel{W(z_m)}{\to }(z_m-1,z_n+1)}}$

According to our framework, an arrow would be drawn from $z_m$ to $z_n$ in the network if the above reaction was part of the system. Another example would be two molecules $z_m$, $z_n$ that bind to form a complex $z_l$

${\displaystyle {\displaystyle (z_m,z_n,z_l)\stackrel{W(z_m,z_n)}{\to }(z_m-1,z_n-1,z_l+1).}}$

If the above reaction was part of the system, an arrow would be drawn from $z_m$ to $z_n$, from $z_n$ to $z_m$, and from both $z_m$, $z_n$ to $z_l$.

In order to derive Equation 2, we need to assume that no components in ${\mathbf{z}}_{\text{aff}}^c$ are part of a co-transition event with any variable in ${\mathbf{z}}_{\text{aff}}$. We thus begin by proving the following lemma.

Lemma 1. Let Z_k be a component that is not affected by X or Y. For any network of the class in Appendix 1—figure 21B (Equations 4.3; 4.4), there exists another network with the exact same dynamics in X, Y, and Z_k, but where no components in ${\mathbf{z}}_{\text{aff}}^c$ are part of a co-transition event with any component in ${\mathbf{z}}_{\text{aff}}$.

Proof of Lemma 1. Let the following co-transition reaction be part of the system, where $\{a_k\}$ are variables in ${\mathbf{z}}_{\text{aff}}$ and $\{b_k\}$ are variables in ${\mathbf{z}}_{\text{aff}}^c$

(5.1)

${\displaystyle {\displaystyle (\mathbf{a},\mathbf{b})\stackrel{W(\mathbf{b})}{\to }(\mathbf{a}+{\mathbf{d}}_a,\mathbf{b}+{\mathbf{d}}_b).}}$

By definition of ${\mathbf{z}}_{\text{aff}}$ and ${\mathbf{z}}_{\text{aff}}^c$, the reaction rate $W$ in the above reaction cannot depend on the variables $\{a_k\}$. However, because a is part of ${\mathbf{z}}_{\text{aff}}$, the $\{a_k\}$ variables must be affected by X or Y. Therefore, there must exist one or more reactions, labeled with $i\in \{1,2,\dots \}$, in the system, that change a with reaction rates that depend on the variables affected by X or Y

(5.2)

${\displaystyle {\displaystyle \mathbf{a}\stackrel{W_i\left({\mathbf{z}}_{\text{aff}},{\mathbf{z}}_{\text{aff}}^c\right)}{\to }\mathbf{a}+{\mathbf{d}}_i\text{for}i=1,2,\dots }}$

Note that the above reactions cannot make changes to b; otherwise, the $\{b_k\}$ variables would not be in the variables not affected by X or Y. We now consider an exact copy of the whole system z and the system reactions, with the change that the a variables are decomposed into two sets of mock variables. Specifically, we define the variables $\{a_k^{int}\}$ and $\{a_k^b\}$ that undergo the following reactions

${\displaystyle {\displaystyle \begin{array}{rl}({\mathbf{a}}^b,\mathbf{b})& \stackrel{W(\mathbf{b})}{\to }\left({\mathbf{a}}^b+{\mathbf{d}}_a,\mathbf{b}+{\mathbf{d}}_b\right)\\ {\mathbf{a}}^{int}& \stackrel{W_i\left({\mathbf{z}}_{\text{aff}},{\mathbf{z}}_{\text{aff}}^c\right)}{\to }{\mathbf{a}}^{int}+{\mathbf{d}}_i\text{for}i=1,2,\dots ,\end{array}}}$

which correspond to the reactions in Equations 5.1; 5.2. We replace all the explicit a dependencies in the reaction rates of the system with ${\mathbf{a}}^{int}+{\mathbf{a}}^b$. As a result, all the reactions that depended on a are unchanged, but now we can put the ${\mathbf{a}}^b$ variables in ${\mathbf{z}}_{\text{aff}}^c$, and we can put the ${\mathbf{a}}^{int}$ variables in ${\mathbf{z}}_{\text{aff}}$. We are left with a system where the ${\mathbf{a}}^{int}$ are not part of a co-transition event with variables in ${\mathbf{z}}_{\text{aff}}^c$, and the dynamics of X and Y remain unchanged. Moreover, none of the reaction rates that govern the dynamics of any component in the original ${\mathbf{z}}_{\text{aff}}^c$ have been altered, meaning the dynamics of Z_k remain unchanged. Such a decomposition can be done for any co-transition reaction that involves components from ${\mathbf{z}}_{\text{aff}}$ and ${\mathbf{z}}_{\text{aff}}^c$.

We now let $Z_k$ correspond to another component of interest in the system, as a continuous-time Markov process. It can be any stochastic process that can be measured, like the abundance of a molecular species in the network, the size of the cell, or any parameter that can influence the reaction rates of the system.

Theorem 1. Let $Z_k$ be a component in ${\mathbf{z}}_{\text{aff}}^c$. If the averages $⟨x_t⟩$, $⟨y_t⟩$ and the covariances $\text{Cov}(x_t,z_{k,t})$, $\text{Cov}(x_t,z_{k,t})$ over the ensemble have reached a stationary state (i.e. they have become constant over time), then ${\eta }_{x{z}_k}={\eta }_{y{z}_k}$.

Proof of Theorem 1. We condition on the components not affected by X or Y, ${\mathbf{z}}_{\text{aff}}^c[-\mathrm{\infty },t]$. This corresponds to a hypothetical system where all the variables in ${\mathbf{z}}_{\text{aff}}^c$ become deterministic time-varying signals $\{z_k(t)\}$. The reactions governing X and Y in the conditional system become

(5.3)

${\displaystyle {\displaystyle \begin{array}{cc}x\stackrel{R({\mathbf{z}}_{\text{aff}},t)}{\to }x+1& y\stackrel{\alpha R({\mathbf{z}}_{\text{aff}},t)}{\to }y+1\\ x\stackrel{x\beta (t)}{\to }x-1& y\stackrel{y\beta (t)}{\to }y-1\end{array}}}$

where $R\left({\mathbf{z}}_{\text{aff}},t\right)=R\left({\mathbf{z}}_{\text{aff}},{\mathbf{z}}_{\text{aff}}^c(t)\right)$ and $\beta (t)=\beta \left({\mathbf{z}}_{\text{aff}}^c(t)\right)$ now have an explicit time dependence from the conditioned history ${\mathbf{z}}_{\text{aff}}^c[-\mathrm{\infty },t]$. We let $A$ be the set of all integers k such that the k-th reaction in Equation 4.1 leads to a change in at least one of the components in ${\mathbf{z}}_{\text{aff}}$. In the conditional probability space, the components affected in ${\mathbf{z}}_{\text{aff}}$ follow the following reactions

(5.4)

${\displaystyle {\displaystyle {\mathbf{z}}_{\text{aff}}\stackrel{W_k({\mathbf{z}}_{\text{aff}},t)}{\to }{\mathbf{z}}_{\text{aff}}+{\mathbf{d}}_k\mathrm{\forall }k\in A,}}$

where $W_k({\mathbf{z}}_{\text{aff}},t)=W_k\left({\mathbf{z}}_{\text{aff}},{\mathbf{z}}_{\text{aff}}^c(t)\right)$ now has an explicit time dependence from the conditioned history of the variables not affected by X or Y. Note that if some components in ${\mathbf{z}}_{\text{aff}}^c$ were part of a co-transition event with components in ${\mathbf{z}}_{\text{aff}}$, then Equations 5.3; 5.4 would not hold. This is because conditioning on the history of those extrinsic variables effectively conditions on those birth events in ${\mathbf{z}}_{\text{aff}}$ that are caused by those co-transitions. However, from Lemma 1, we can always work with another network in which there are no such co-transition events, and where the dynamics of X and Y remain unchanged.

This conditional system follows the following master equation

${\displaystyle {\displaystyle {\displaystyle \begin{array}{ll}{\displaystyle \frac{d}{dt}}& P\left(x,y,{\mathbf{z}}_{\mathrm{a}\mathrm{f}\mathrm{f}},t\mid {\mathbf{z}}_{\mathrm{a}\mathrm{f}\mathrm{f}}^c[-\mathrm{\infty },t]\right)=\\ & {\displaystyle \sum _{k\in A}\left[W_k\left({\mathbf{z}}_{\mathrm{a}\mathrm{f}\mathrm{f}}-{\mathbf{d}}_k,t\right)P\left(x,y,{\mathbf{z}}_{\mathrm{a}\mathrm{f}\mathrm{f}}-{\mathbf{d}}_k,t\mid {\mathbf{z}}_{\mathrm{a}\mathrm{f}\mathrm{f}}^c[-\mathrm{\infty },t]\right)-W_k\left({\mathbf{z}}_{\mathrm{a}\mathrm{f}\mathrm{f}},t\right)P\left(x,y,{\mathbf{z}}_{\mathrm{a}\mathrm{f}\mathrm{f}},t\mid {\mathbf{z}}_{\mathrm{a}\mathrm{f}\mathrm{f}}^c[-\mathrm{\infty },t]\right)\right]}\\ & {\displaystyle +R\left(x-1,y,{\mathbf{z}}_{\mathrm{a}\mathrm{f}\mathrm{f}},t\right)P\left(x-1,y,{\mathbf{z}}_{\mathrm{a}\mathrm{f}\mathrm{f}},t\mid {\mathbf{z}}_{\mathrm{a}\mathrm{f}\mathrm{f}}^c[-\mathrm{\infty },t]\right)-R\left(x,y,{\mathbf{z}}_{\mathrm{a}\mathrm{f}\mathrm{f}},t\right)P\left(x,y,{\mathbf{z}}_{\mathrm{a}\mathrm{f}\mathrm{f}},t\mid {\mathbf{z}}_{\mathrm{a}\mathrm{f}\mathrm{f}}^c[-\mathrm{\infty },t]\right)}\\ & {\displaystyle +\alpha R\left(x,y-1,{\mathbf{z}}_{\mathrm{a}\mathrm{f}\mathrm{f}},t\right)P\left(x,y-1,{\mathbf{z}}_{\mathrm{a}\mathrm{f}\mathrm{f}},t\mid {\mathbf{z}}_{\mathrm{a}\mathrm{f}\mathrm{f}}^c[-\mathrm{\infty },t]\right)-\alpha R\left(x,y,{\mathbf{z}}_{\mathrm{a}\mathrm{f}\mathrm{f}},t\right)P\left(x,y,{\mathbf{z}}_{\mathrm{a}\mathrm{f}\mathrm{f}};t\mid {\mathbf{z}}_{\mathrm{a}\mathrm{f}\mathrm{f}}^c[-\mathrm{\infty },t]\right)}\\ & {\displaystyle +(x+1)\beta (t)P\left(x+1,y,{\mathbf{z}}_{\mathrm{a}\mathrm{f}\mathrm{f}},t\mid {\mathbf{z}}_{\mathrm{a}\mathrm{f}\mathrm{f}}^c[-\mathrm{\infty },t]\right)-x\beta (t)P\left(x,y,{\mathbf{z}}_{\mathrm{a}\mathrm{f}\mathrm{f}},t\mid {\mathbf{z}}_{\mathrm{a}\mathrm{f}\mathrm{f}}^c[-\mathrm{\infty },t]\right)}\\ & +(y+1)\beta (t)P\left(x,y+1,{\mathbf{z}}_{\mathrm{a}\mathrm{f}\mathrm{f}},t\mid {\mathbf{z}}_{\mathrm{a}\mathrm{f}\mathrm{f}}^c[-\mathrm{\infty },t]\right)-y\beta (t)P\left(x,y,{\mathbf{z}}_{\mathrm{a}\mathrm{f}\mathrm{f}},t\mid {\mathbf{z}}_{\mathrm{a}\mathrm{f}\mathrm{f}}^c[-\mathrm{\infty },t]\right).\end{array}}}}$

We consider the averages of X and Y conditioned on the upstream history, $\overline{x}(t)=E\left[x_t|{\mathbf{z}}_{\text{aff}}^c[-\mathrm{\infty },t]\right]$ and $\overline{y}(t)=E\left[y_t|{\mathbf{z}}_{\text{aff}}^c[-\mathrm{\infty },t]\right]$, where $x_t$ and $y_t$ are the X and Y abundances at time t. From the above master equation, the time-evolution for these first moments can be derived (Joly-Smith et al., 2021; Hilfinger and Paulsson, 2011)

(5.5)

${\displaystyle {\displaystyle \frac{d\overline{x}}{dt}=\overline{R}(t)-\overline{x}\overline{\beta }(t)\mathrm{\&}\frac{d\overline{y}}{dt}=\alpha \overline{R}(t)-\overline{y}\overline{\beta }(t),}}$

where $\overline{R}(t)=E\left[R(x_t,y_t,{\mathbf{z}}_{\text{aff}},t)|{\mathbf{z}}_{\text{ aff}}^c[-\mathrm{\infty },t]\right]$ and $\overline{\beta }(t)=E\left[\beta ({\mathbf{z}}_{\text{aff}}^c)|{\mathbf{z}}_{\text{aff }}^c[-\mathrm{\infty },t]\right]$ are the average production and degradation rates conditioned on the history of the variables not affected by X or Y. Note that $\overline{\beta }(t)=E\left[\beta ({\mathbf{z}}_{\text{aff}}^c)|{\mathbf{z}}_{\text{aff }}^c[-\mathrm{\infty },t]\right]=\beta ({\mathbf{z}}_{\text{aff}}^c(t))=\beta (t)$, because the time trajectory of ${\mathbf{z}}_{\text{aff}}^c$ is set through the conditioning on the upstream history. We can then take the expectation of $\overline{x}$ over all possible histories of ${\mathbf{z}}_{\text{aff}}^c$ to get

(5.6)

${\displaystyle {\displaystyle E[\overline{x}(t){]}_{\text{histories}}=E{\left[E\left[x_t\mid {\mathbf{z}}_{\text{aff}}^c[-\mathrm{\infty },t]\right]\right]}_{\text{histories}}=⟨x_t⟩,}}$

which follows from the law of total expectation. We now let $Z_k$ correspond to any component in the network that is not affected by X or Y. It can be a molecular abundance, concentration, or another stochastic cellular variable like the growth rate of the cell. It follows that

(5.7)

${\displaystyle {\displaystyle \begin{array}{rl}E{\left[\overline{x}(t)z_k(t)\right]}_{\text{histories}}& =E{\left[E\left[x_t\mid {\mathbf{z}}_{\text{aff }}^c[-\mathrm{\infty },t]\right]\cdot E\left[z_{k,t}\mid {\mathbf{z}}_{\text{aff }}^c[-\mathrm{\infty },t]\right]\right]}_{\text{histories }}\\ & =E{\left[E\left[x_t{z}_{k,t}\mid {\mathbf{z}}_{\text{aff }}^c[-\mathrm{\infty },t]\right]\right]}_{\text{histories }}=⟨x_t{z}_{k,t}⟩\end{array},}}$

where the second step comes from the fact that conditioning on the history of ${\mathbf{z}}_{\text{aff}}^c$ effectively also conditions on the history of $Z_k$ with $z_k(t)=E\left[z_{k,t}|{\mathbf{z}}_{\text{aff}}^c[-\mathrm{\infty },t]\right]$ (so x and $z_k$ are independent when conditioning on the ${\mathbf{z}}_{\text{aff}}^c$ history), the last step follows from the law of total expectation, and $z_{k,t}$ is the measured amount of $Z_k$ at time t. From Equations 5.6; 5.7, it follows that

(5.8)

${\displaystyle \begin{array}{c}\text{Cov}(x_t,z_{k,t})=\text{Cov}(\overline{x}(t),{\overline{z}}_k(t)),\end{array}}$

where ${\overline{z}}_k(t)=E\left[z_{k,t}|{\mathbf{z}}_{\text{aff}}^c[-\mathrm{\infty },t]\right]=z_k(t)$, where the last step comes from the fact that the time trajectory of $z_k$ is set through the conditioning of the upstream histories. Intuitively, Equation 5.8 says that when X does not affect $Z_k$, the stochastic fluctuations of X average out when taking the covariance between X and $Z_k$. Strikingly, this is independent of any type of feedback that X may impose through interactions in the cloud of components $\mathbf{z}(t)$. As a result, the same should hold for Y, and since $\overline{y}(t)$ is governed by the same differential equation as $\overline{x}$, the intrinsic fluctuations that differentiate X and Y will average out when taking the covariances with $Z_k$.

That is, dividing the right equation in Equation 5.5 with α, we write the general solution for $\overline{x}(t)$ and $\overline{y}(t)/\alpha$, and find

(5.9)

${\displaystyle {\displaystyle \overline{x}(t)=\overline{x}(0)e^{-{\int }_0^t\beta (u)du}+{\int }_0^t{e}^{-\left({\int }_0^t\beta (u)du-{\int }_0^{t^{\prime }}\beta (v)dv\right)}R(t^{\prime })d{t}^{\prime }}}$

(5.10)

${\displaystyle {\displaystyle \overline{y}(t)/\alpha =\overline{y}(0)e^{-{\int }_0^t\beta (u)du}/\alpha +{\int }_0^t{e}^{-\left({\int }_0^t\beta (u)du-{\int }_0^{t^{\prime }}\beta (v)dv\right)}R(t^{\prime })d{t}^{\prime }}}$

Taking the average over all histories and subtracting the equations we have

(5.11)

${\displaystyle {\displaystyle \begin{array}{rl}E[\overline{x}(t){]}_{\text{histories}}-E[\overline{y}(t){]}_{\text{histories}}& {\displaystyle =E{\left[\left(\overline{x}(0)-\overline{y}(0)/\alpha \right)e^{-{\int }_0^t\beta (u)du}\right]}_{\text{histories}}}\\ \Rightarrow ⟨x_t⟩-⟨y_t⟩/\alpha & {\displaystyle =E{\left[\left(\overline{x}(0)-\overline{y}(0)/\alpha \right)e^{-{\int }_0^t\beta (u)du}\right]}_{\text{histories}},}\end{array}}}$

where in the second step we used Equation 5.6 which also holds for $y$ by symmetry. We now invoke the requirement that the averages $⟨x_t⟩$ and $⟨y_t⟩$ are stationary, meaning they are constant over time. In that case, $⟨x_t⟩-⟨y_t⟩/\alpha$ is constant over time, and so

(5.12)

${\displaystyle {\displaystyle ⟨x_t⟩-⟨y_t⟩/\alpha =\underset{t\to \mathrm{\infty }}{lim}\left(⟨x_t⟩-⟨y_t⟩/\alpha \right).}}$

Substituting Equation 5.11, we thus have

(5.13)

${\displaystyle {\displaystyle ⟨x_t⟩-⟨y_t⟩/\alpha =E{\left[\left(\overline{x}(0)-\overline{y}(0)/\alpha \right)e^{-{\int }_0^{\mathrm{\infty }}\beta (u)du}\right]}_{\text{histories}}.}}$

Now, we must have ${\int }_0^{\mathrm{\infty }}\beta (u)du\to \mathrm{\infty }$ when $⟨\beta {⟩}_t=\underset{T\to \mathrm{\infty }}{lim}\frac{1}{T}{\int }_0^T\beta (u)du>0$. Therefore, the right-hand side of Equation 5.13 becomes 0:

(5.14)

${\displaystyle {\displaystyle ⟨x_t⟩-⟨y_t⟩/\alpha =E{\left[\left(\overline{x}(0)-\overline{y}(0)/\alpha \right)e^{-{\int }_0^{\mathrm{\infty }}\beta (u)du}\right]}_{\text{histories}}=0,}}$

Similarly, we now multiply Equations 5.10; 5.9 with $z_k(t)$, average over all histories, and subtract to obtain

(5.15)

${\displaystyle {\displaystyle \begin{array}{rl}llE[\overline{x}(t)z_k(t){]}_{\text{histories}}-E[\overline{y}(t)z_k(t){]}_{\text{histories}}& {\displaystyle =E{\left[\left(\overline{x}(0)-\overline{y}(0)/\alpha \right)z_k(t)e^{-{\int }_0^t\beta (u)du}\right]}_{\text{histories}}}\\ \Rightarrow ⟨x_t{z}_{k,t}⟩-⟨y_t{z}_{k,t}⟩/\alpha & {\displaystyle =E{\left[\left(\overline{x}(0)-\overline{y}(0)/\alpha \right)z_k(t)e^{-{\int }_0^t\beta (u)du}\right]}_{\text{histories}}=0,}\end{array}}}$

where in the second step we used Equation 5.7 which also holds for y by symmetry. The last step follows when $⟨x_t{z}_{k,t}⟩$ and $⟨y_t{z}_{k,t}⟩$ have reached stationarity and are constant over time, along with $⟨\beta {⟩}_t=\underset{T\to \mathrm{\infty }}{lim}\frac{1}{T}{\int }_0^T\beta (u)du>0$.

It then follows from Equations 5.14; 5.15 that

(5.16)

${\displaystyle {\displaystyle ⟨x⟩=⟨y⟩/\alpha \mathrm{\&}\text{Cov}(x,z_k)=\text{Cov}(y,z_k)/\alpha .}}$

Dividing $\text{Cov}(x,z_k)$ by $⟨x⟩⟨z⟩$, and $\text{Cov}(y,z_k)/\alpha$ by $⟨y⟩⟨z⟩/\alpha$, we find ${\eta }_{x{z}_k}={\eta }_{y{z}_k}$.

If the reporter Y is engineered to be passive (i.e., it does not affect components in the network), then a violation of Equation 2 would imply that X affects $Z_k$. Otherwise, such a violation would imply that X or Y affect $Z_k$.

Thus far, we assumed that all the components not affected by X or Y are part of a continuous-time Markov chain. This was in order to make a rigorous definition of causal interaction in our framework as a path in the topology of the transition rates. Alternatively, if we relax the requirement that the components not affected by X or Y be Markov chains (i.e. they can be a set of arbitrary stochastic processes), we can operationally define ‘no causal interaction from X or Y’ to mean that we can condition on the history of those stochastic processes and write down Equations 5.3; 5.4. We can then operationally define any violation of Equation 2 as a ‘causal interaction from X or Y’ to a stochastic process $Z_k$.