Path mutual information for a class of biochemical reaction networks

TL;DR

This paper introduces a method to estimate the mutual information between the complete trajectories of two interacting molecular species in biochemical reaction networks, addressing mathematical challenges in stochastic systems.

Contribution

It presents a novel approach for calculating path mutual information in biochemical networks, enabling better understanding of cellular information processing.

Findings

Method successfully applied to three case studies

Provides a quantitative measure of interdependence in biochemical systems

Addresses mathematical challenges in stochastic reaction network analysis

Abstract

Living cells encode and transmit information in the temporal dynamics of biochemical components. Gaining a detailed understanding of the input-output relationship in biological systems therefore requires quantitative measures that capture the interdependence between complete time trajectories of biochemical components. Mutual information provides such a measure but its calculation in the context of stochastic reaction networks is associated with mathematical challenges. Here we show how to estimate the mutual information between complete paths of two molecular species that interact with each other through biochemical reactions. We demonstrate our approach using three simple case studies.

Figures

Equations

$$\sum _{i=1}^D{\alpha }_{k,i}{\text{SS}}^{(i)}\rightharpoonup \sum _{i=1}^D{\beta }_{k,i}{\text{SS}}^{(i)}$$

$$h_k({\text{SS}}_t)=c_k\prod _{i=1}^D\left(\genfrac{}{}{0ex}{}{{\text{SS}}_t^{(i)}}{{\alpha }_{k,i}}\right),$$

$${\text{SS}}_t={\text{SS}}_0+\sum _{k=1}^K{N}_k\left({\int }_0^t{h}_k({\text{SS}}_s)\mathrm{d}s\right){\nu }_k,$$

$$\mathbb{I}(A_0^t,B_0^t)=\mathbb{E}\text{bqty}\log \frac{\mathrm{d}{\mathbb{P}}^{AB}}{\mathrm{d}({\mathbb{P}}^A\times {\mathbb{P}}^B)}.$$

$$A_t$$

$$B_t$$

$$A_t$$

$$B_t$$

$$\begin{split}\mathrm{d}\pi^{B}(a,t)&=\sum_{k\in\mathcal{R}_{A}}\Big{[}h_{k}(a-\nu_{k}^{A},B_{t})\pi^{B}(a-\nu_{k}^{A},t)\\ &\quad\quad\quad\quad\quad-h_{k}(a,B_{t})\pi^{B}(a,t)\Big{]}\mathrm{d}t\\ &\quad-\sum_{k\in\mathcal{R}_{B|A}}\pqty{h_{k}(a,B_{t})-\lambda_{k}^{A}(t)}\pi^{B}(a,t)\mathrm{d}t\\ &\quad+\sum_{k\in\mathcal{R}_{B|A}}\frac{h_{k}(a,B_{t})-\lambda_{k}^{A}(t)}{\lambda_{k}^{A}(t)}\pi^{B}(a,t)\mathrm{d}N_{k}(t)\end{split}$$

$$\frac{\mathrm{d}{\mathbb{Q}}^X}{\mathrm{d}{\hat{\mathbb{Q}}}^X}=\frac{{\prod }_{k=1}^K{\prod }_{j=1}^{N_k(t)}{\eta }_k(T_{k,j}^{-})e^{-{\int }_0^t{\eta }_k(s)\mathrm{d}s}}{{\prod }_{k=1}^K{\prod }_{j=1}^{N_k(t)}{\hat{\eta }}_k(T_{k,j}^{-})e^{-{\int }_0^t{\hat{\eta }}_k(s)\mathrm{d}s}},$$

$$\begin{split}\frac{\mathrm{d}\mathbb{P}^{AB}}{\mathrm{d}(\mathbb{P}^{A}\times\mathbb{P}^{B})}&=\prod_{k=1}^{K}\prod_{j=1}^{N_{k}(t)}\lambda_{k}^{AB}(T_{k,j}^{-})e^{-\int_{0}^{t}\lambda_{k}^{AB}(s)\mathrm{d}s}\\ &\quad\times\Bigg{(}\prod_{k\in\mathcal{R}_{A}}\prod_{j=1}^{N_{k}(t)}\lambda_{k}^{A}(T_{k,j}^{-})e^{-\int_{0}^{t}\lambda_{k}^{A}(s)\mathrm{d}s}\\ &\quad\times\prod_{k\in\mathcal{R}_{B}}\prod_{j=1}^{N_{k}(t)}\lambda_{k}^{B}(T_{k,j}^{-})e^{-\int_{0}^{t}\lambda_{k}^{B}(s)\mathrm{d}s}\Bigg{)}^{-1}.\end{split}$$

$$\begin{array}{rl}& \log \frac{\mathrm{d}{\mathbb{P}}^{AB}}{\mathrm{d}({\mathbb{P}}^A\times {\mathbb{P}}^B)}=\\ & \sum _{k=1}^K\sum _{j=1}^{N_k(t)}\log {\lambda }_k^{AB}(T_{k,j}^{-})-{\int }_0^t{\lambda }_k^{AB}(s)\mathrm{d}s\\ & -\sum _{k\in {\mathcal{R}}_A}\sum _{j=1}^{N_k(t)}\log {\lambda }_k^A(T_{k,j}^{-})+{\int }_0^t{\lambda }_k^A(s)\mathrm{d}s\\ & -\sum _{k\in {\mathcal{R}}_B}\sum _{j=1}^{N_k(t)}\log {\lambda }_k^B(T_{k,j}^{-})+{\int }_0^t{\lambda }_k^B(s)\mathrm{d}s.\end{array}$$

$$\mathbb{I}(A_0^t,B_0^t)=\mathbb{E}\text{bqty}\log \frac{\mathrm{d}{\mathbb{P}}^{AB}}{\mathrm{d}({\mathbb{P}}^A\times {\mathbb{P}}^B)}$$

$$\displaystyle=\mathbb{E}\Big{[}\sum_{k\in\mathcal{R}_{A|B}}\int_{0}^{t}\pqty{\log\lambda_{k}^{AB}(s^{-})-\log\lambda_{k}^{A}(s^{-})}\mathrm{d}N_{k}(s)$$

$$-{\int }_0^t\text{pqty}{\lambda }_k^{AB}(s)-{\lambda }_k^A(s)\mathrm{d}s$$

$$+\sum _{k\in {\mathcal{R}}_{B|A}}{\int }_0^t\text{pqty}\log {\lambda }_k^{AB}(s^{-})-\log {\lambda }_k^B(s^{-})\mathrm{d}{N}_k(s)$$

$$\displaystyle\quad\quad\quad\quad\quad\quad\quad-\int_{0}^{t}\pqty{\lambda_{k}^{AB}(s)-\lambda_{k}^{B}(s)}\mathrm{d}s\Big{]},$$

$$\mathbb{I}(A_0^t,B_0^t)=$$

$$-\sum _{k\in {\mathcal{R}}_{A|B}}{\int }_0^t\mathbb{E}\text{bqty}{\lambda }_k^A(s)\log {\lambda }_k^A(s)\mathrm{d}s$$

$$-\sum _{k\in {\mathcal{R}}_{B|A}}{\int }_0^t\mathbb{E}\text{bqty}{\lambda }_k^B(s)\log {\lambda }_k^B(s)\mathrm{d}s,$$

$$\#1:\varnothing \stackrel{{\gamma }_A}{\rightharpoonup }\mathrm{A}$$

$$\#2:\mathrm{A}\stackrel{{\delta }_A}{\rightharpoonup }\varnothing$$

$$\displaystyle\mathbb{I}({A}_{0}^{t},{B}_{0}^{t})=\mathbb{E}\Big{[}\int_{0}^{t}\pqty{\log\pqty{\gamma_{B}A_{s}}-\log(\gamma_{B}\hat{A}_{1}(s))}\mathrm{d}N_{3}(s)$$

$$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,-\gamma_{B}\int_{0}^{t}\pqty{A_{s}-\hat{A}_{1}(s)}\mathrm{d}s\Big{]},$$

$$\mathrm{d}{\pi }^B(a,t)$$

$$+{\delta }_A\text{bqty}(a+1){\pi }^B(a+1,t)-a{\pi }^B(a,t)\mathrm{d}t$$

$$-{\gamma }_B\text{bqty}a-{\hat{A}}_1(t){\pi }^B(a,t)\mathrm{d}t$$

$$+\frac{a-{\hat{A}}_1(t)}{{\hat{A}}_1(t)}{\pi }^B(a,t)\mathrm{d}{N}_3(t).$$

$$\mathrm{d}{\hat{A}}_1=$$

$$+\frac{{\hat{A}}_2(t)-{\hat{A}}_1^2(t)}{{\hat{A}}_1(t)}\mathrm{d}{N}_3(t).$$

$$\mathrm{d}{\hat{A}}_2=$$

$$\displaystyle-2\gamma_{B}\frac{\hat{A}_{2}(t)}{\hat{A}_{1}(t)}\pqty{\hat{A}_{2}(t)-\hat{A}_{1}^{2}(t)}\Big{]}\mathrm{d}t$$

$$+2\frac{{\hat{A}}_2(t)}{{\hat{A}}_1^2(t)}\text{pqty}{\hat{A}}_2(t)-{\hat{A}}_1^2(t)\mathrm{d}{N}_3(t).$$

$$\#1:\mathrm{A}\stackrel{{\gamma }_A}{\rightharpoonup }2\mathrm{A}$$

$$\#2:\mathrm{A}+\mathrm{B}\stackrel{{\delta }_A}{\rightharpoonup }\mathrm{B}$$

Full text

Path mutual information for a class of biochemical reaction networks

Lorenzo Duso

Max Planck Institute of Molecular Cell Biology and Genetics, 01307 Dresden, Germany

Center for Systems Biology Dresden, 01307 Dresden, Germany

Christoph Zechner

Max Planck Institute of Molecular Cell Biology and Genetics, 01307 Dresden, Germany

Center for Systems Biology Dresden, 01307 Dresden, Germany

Correspondence to: zechner@mpi-cbg.de

**Abstract – Living cells encode and transmit information in the temporal dynamics of biochemical components. Gaining a detailed understanding of the input-output relationship in biological systems therefore requires quantitative measures that capture the interdependence between complete time trajectories of biochemical components. Mutual information provides such a measure but its calculation in the context of stochastic reaction networks is associated with mathematical challenges. Here we show how to estimate the mutual information between complete paths of two molecular species that interact with each other through biochemical reactions. We demonstrate our approach using three simple case studies. **

1 Introduction

The ability to continuously sense and respond to a wide spectrum of environmental signals and stresses is a hallmark of living systems. At the same time, biochemical processes inside cells are often significantly affected by random fluctuations [1, 2], which stands in contrast with the remarkable robustness of biological systems. Understanding how living cells reliably transmit and process information is a fundamental problem in biology, which gained wide attention in the past.

The combination of probabilistic methods and information theory [3] provides a rich framework to gain a quantitative understanding of signal processing in biological systems. Indeed, the use of information theoretical quantities such as the Shannon entropy and mutual information (MI) has led to insights into diverse biological systems [4, 5], ranging from neural networks [6, 7], to biomolecular signalling and decision making systems [8, 9, 10, 11, 12, 13]. The mutual information and related quantities such as the channel capacity have proven particularly useful in the context of biology, since they capture how information is transmitted between certain inputs (e.g., inducer molecules) and outputs (e.g., downstream targets) through cascades of biochemical reactions.

Experimental studies have demonstrated that living cells frequently encode information in the temporal dynamics of biochemical processes [14, 15]. In this case, the mutual information calculated at individual time-points would significantly underestimate the information that is transmitted across a system. A key advantage of information theoretic quantities is that they can be defined in very general terms such that they apply also to time-spanning objects. Extensions of the mutual information to trajectories have been proposed [16, 17, 18, 19] as well as the closely related transfer entropy [20, 21, 22]. However, the analysis of path-related information measures is mathematically more demanding, especially in the context of continuous-time jump processes such as frequently encountered in biochemical systems. While the information theoretical foundation has been established for such processes [23, 16, 17], the explicit calculation of the path mutual information for concrete systems remains a challenge. Previous studies have addressed this problem either numerically using particle filters [17, 19] or by approximating the original jump process by a continuous Gaussian process [18].

In this work we present an efficient strategy to calculate the path mutual information between two chemical species that are coupled through biochemical reactions. The method does not rely on continuous approximations of the underlying process and therefore correctly takes into account the discrete nature of biochemical networks. Our theoretical results build upon the recently proposed marginal process framework [24, 25], which allows us to estimate the path mutual information by combining conventional stochastic simulations with the solution of a filtering problem.

The structure of the remaining paper is as follows. We review some basic definitions and concepts of stochastic chemical kinetics in Section 2.1. In Section 2.2 we formally state the problem of calculating path mutual information for a class of biochemical networks and in 2.3 we show how it can be addressed using the marginal process framework. We provide in Section 2.4 explicit expressions for the Radon-Nikodym derivative between the joint and marginal path measures which then allows us to state our main result in Section 2.5. In Section 2.6 we briefly relate the obtained expression to previous theoretical results from the literature. Finally, we apply the method to study information transmission in three simple case studies (Section 3).

2 Theoretical Results

2.1 Stochastic chemical kinetics

We consider a stochastic reaction network $\mathcal{R}_{\SS}$ that describes the time evolution of $D$ chemical species $\mathrm{\SS}^{(1)},\ldots,\mathrm{\SS}^{(D)}$ and $K$ reaction channels. The network can be defined by a set of stoichiometric equations

[TABLE]

for $k=1,\ldots,K$. The non-negative integer numbers $\alpha_{{k},{i}}$ and $\beta_{{k},{i}}$ are respectively the reactant and product stoichiometric coefficient of species $i$ in the reaction $k$. We introduce the stochastic process $\SS=(\SS_{t})_{t\geq 0}$, which takes values in $\mathbb{N}_{0}^{D}$ and is subject to the reaction dynamics in (1). The state vector $\SS_{t}$ collects the discrete copy numbers of each species at time $t$ and it changes by $\nu_{k}=\beta_{k}-\alpha_{k}$ when reaction $k$ occurs. Each reaction channel $k$ is equipped with a propensity function $h_{k}(\SS_{t})$, which sets the firing rate of this reaction. Throughout this work, we consider mass-action kinetics such that

[TABLE]

with $c_{k}$ as the stochastic rate constant of reaction $k$. In this setting, $\SS$ admits a continuous-time Markov chain (CTMC), whose time evolution satisfies an integral equation of the form

[TABLE]

where $N_{k}$ denotes a unit Poisson process counting the occurrences of reaction $k$ until time $t$ [26]. We define by ${\SS}_{0}^{t}$ a complete path of $\SS$ collecting all information about the types and firing times of the reactions that happen between time zero and $t$. Moreover, let $\mathbb{P}^{X}$ denote the probability measure over the complete path ${\SS}_{0}^{t}$ considered on the natural filtration $\mathcal{F}_{t}^{X}$ generated by $\SS$. Note that for a given initial condition $\SS_{0}$, exact realizations of ${\SS}_{0}^{t}$ can be simulated using the Stochastic Simulation Algorithm (SSA) [27].

2.2 Information transmission between two chemical species

The goal of this work is to calculate the mutual information between complete paths of two chemical species, which we denote by $\mathrm{A}$ and $\mathrm{B}$, respectively. For simplicity, we restrict ourselves to the scenario where the network $\mathcal{R}_{\SS}$ comprises only $\mathrm{A}$ and $\mathrm{B}$ and no additional species. The network evolves through an arbitrary number of reaction channels, whereas only reactions involving both $\mathrm{A}$ and $\mathrm{B}$ (as reactant or product) will lead to an exchange of information among the two. We refer to these reactions as coupling reactions. Here we consider coupling reactions that modify either $\mathrm{A}$ and $\mathrm{B}$ but not both at the same time. This involves reactions of catalytic or annihilating form such as $\mathrm{A}\rightarrow\mathrm{A}+\mathrm{B}$ or $\mathrm{A}+\mathrm{B}\rightarrow\mathrm{A}$, whereas conversions such as $\mathrm{A}\rightarrow\mathrm{B}$ are not considered. As a consequence, the total set of reactions can be split into two disjoint sets $\mathcal{R}_{\SS}=\{\mathcal{R}_{A},\mathcal{R}_{B}\}$ where $\mathcal{R}_{A}$ and $\mathcal{R}_{B}$ are the reactions that exclusively modify $\mathrm{A}$ or $\mathrm{B}$, respectively. In the following, we refer to this property as smooth coupling between $\mathrm{A}$ and $\mathrm{B}$.

Without loss of generality, we arrange the state vector such that $\SS_{t}=\pqty{A_{t},B_{t}}$. We define the partial paths ${A}_{0}^{t}$ and ${B}_{0}^{t}$ and corresponding sub-filtrations $\mathcal{F}_{t}^{A},\mathcal{F}_{t}^{B}\subset\mathcal{F}_{t}^{X}=\mathcal{F}_{t}^{AB}$ generated by $A_{t}$ and $B_{t}$. Notice that since $\mathrm{A}$ and $\mathrm{B}$ are smoothly coupled, $\mathcal{F}_{t}^{A}$ and $\mathcal{F}_{t}^{B}$ contain information exclusively about reactions in $\mathcal{R}_{A}$ and $\mathcal{R}_{B}$, respectively. Using these definitions, our goal is to characterize information transmission between paths ${A}_{0}^{t}$ and ${B}_{0}^{t}$ via their path mutual information (path-MI)

[TABLE]

In (4), the term inside the logarithm corresponds to the Radon-Nikodym derivative of the joint path measure $\mathbb{P}^{AB}$ with respect to the product of the marginal path measures $\mathbb{P}^{A}$ and $\mathbb{P}^{B}$. The latter can be thought of as the probability laws that describe the time evolution of only $\mathrm{A}$ or $\mathrm{B}$, respectively. We and others have previously shown how such marginal process models can be constructed for general reaction networks [28, 24]. While in those studies, it was used predominantly for the purpose of model reduction and stochastic simulation, it will now serve us to find explicit expressions of the Radon-Nikodym derivative in (4). In the following section, we will show how the marginal process dynamics can be obtained for the considered two-species network.

2.3 Marginal process dynamics

As a first step, we instantiate the dynamics from (3) for the considered network. This yields two coupled integral equations of the form

[TABLE]

where $\nu^{A}_{k}$ and $\nu^{B}_{k}$ are respectively the stoichiometric change coefficients corresponding to $\mathrm{A}$ and $\mathrm{B}$. Note that we have introduced $\lambda_{k}^{AB}(t)=h_{k}(A_{t},B_{t})$ to emphasize that (5) and (6) are the equations of motion for $A_{t}$ and $B_{t}$ relative to their (joint) filtration $\mathcal{F}_{t}^{AB}$. The corresponding solutions ${\SS}_{0}^{t}=\{{A}_{0}^{t},{B}_{0}^{t}\}$ admit a joint measure $\mathbb{P}^{AB}$. Our goal is to find two analogous equations, each being consistent with the marginal path measures $\mathbb{P}^{A}$ and $\mathbb{P}^{B}$, respectively. Technically, this can be achieved by restricting the two processes to depend only on their own history (i.e., $\mathcal{F}_{t}^{A}$ in case of $A_{t}$ and $\mathcal{F}_{t}^{B}$ in case of $B_{t}$) but not their joint history captured by $\mathcal{F}_{t}^{AB}$. Informally, this can be understood as ”integrating out” the dependency of one path on the other one. It can be shown [24] that relative to the filtrations $\mathcal{F}_{t}^{A}$ and $\mathcal{F}_{t}^{B}$, the counting processes (5) and (6) evolve according to

[TABLE]

where we have defined $\lambda_{k}^{A}(t)=\mathbb{E}\bqty{h_{k}(A_{t},B_{t})\mid\mathcal{F}_{t}^{A}}$ and $\lambda_{k}^{B}(t)=\mathbb{E}\bqty{h_{k}(A_{t},B_{t})\mid\mathcal{F}_{t}^{B}}$. In other words, the original propensities $h_{k}$ are replaced by their expected value conditionally on the sub-filtrations $\mathcal{F}_{t}^{A}$ and $\mathcal{F}_{t}^{B}$. In the information theory literature, these conditional expectations are commonly referred to as causal or filtering estimates, since they provide a reconstruction of the original propensities $h_{k}$ at time $t$ from a continuous path up to time $t$ (i.e., either ${A}_{0}^{t}$ or ${B}_{0}^{t}$). Notice that due to the marginalization, (7) and (8) are now decoupled from each other, such that their solutions are consistent with the product measure $\mathbb{P}^{A}\times\mathbb{P}^{B}$. We point out, however, that the dynamics are no longer Markovian due to the history-dependence of $\lambda_{k}^{A}(t)$ and $\lambda_{k}^{B}(t)$.

While (7) and (8) provide the marginal descriptions we seek for, the calculation of the marginalized propensities $\lambda_{k}^{A}(t)$ and $\lambda_{k}^{B}(t)$ is associated with certain challenges. More precisely, the expectations $\mathbb{E}\bqty{\cdot\mid\mathcal{F}_{t}^{A}}$ and $\mathbb{E}\bqty{\cdot\mid\mathcal{F}_{t}^{B}}$ are taken with respect to the conditional distributions $\pi^{A}(b,t)=P(B_{t}=b\mid\mathcal{F}_{t}^{A})$ and $\pi^{B}(a,t)=P(A_{t}=a\mid\mathcal{F}_{t}^{B})$, for which analytic expressions are not available. We have shown previously [24] that a stochastic differential equation can be found for such conditional distributions. For instance, the equation for $\pi^{B}(a,t)$ would read

[TABLE]

with $\mathcal{R}_{B|A}\subset\mathcal{R}_{B}$ denoting the set of coupling reactions affecting species $\mathrm{B}$ which are driven by $\mathrm{A}$. In (9), the terms $\mathrm{d}N_{k}(t)$ correspond to the differential version of the reaction counters $N_{k}$ at time $t$. A corresponding equation for $\pi^{A}(b,t)$ can be formulated analogously. Since the propensities obey the mass-action law from (2), we further realize that $\lambda_{k}^{A}(t)$ and $\lambda_{k}^{B}(t)$ are just functions of the moments of the associated filtering distributions. The calculation of the marginal propensities is therefore equivalent to the calculation of the moments of $\pi^{A}(b,t)$ and $\pi^{B}(a,t)$. However, depending on the stoichiometry of the considered system, the obtained moment dynamics may not be closed. In this case we can employ suitable moment-closure schemes [29] to derive approximate equations for the conditional moments. We will provide explicit expressions for the time-evolution of $\lambda_{k}^{A}(t)$ and $\lambda_{k}^{B}(t)$ for some examples in Section 3. For further information on the marginal process framework and the calculation of the conditional expectations the reader may refer to [24].

2.4 Radon-Nikodym derivatives and Jacod’s formula

Having discussed the joint and marginal dynamics of $A_{t}$ and $B_{t}$, we are now ready to study the Radon-Nikodym derivative that appears inside (4). To this end, we employ Jacod’s formula, which provides an explicit form of the Radon-Nikodym derivative for multivariate counting processes. Consider a counting process $\SS_{t}$ of the form (3) with natural filtration $\mathcal{F}_{t}^{X}$. We define path measures $\mathbb{Q}^{X}$ and $\widehat{\mathbb{Q}}^{X}$ restricted to $\mathcal{F}_{t}^{X}$ under which $\SS$ has propensitities $\eta_{k}(t)$ and $\hat{\eta}_{k}(t)$, respectively. Note that $\eta_{k}(t)$ and $\hat{\eta}_{k}(t)$ may be history-dependent in general. The Radon-Nikodym derivative of $\mathbb{Q}^{X}$ with respect to $\widehat{\mathbb{Q}}^{X}$ is given by [23]

[TABLE]

where the symbol $T_{k,j}^{-}$ is the left limit to the $j^{th}$ firing time of reaction channel $k$. Therefore, if we set $\mathbb{Q}^{X}$ and $\widehat{\mathbb{Q}}^{X}$ to $\mathbb{P}^{AB}$ and $\mathbb{P}^{A}\times\mathbb{P}^{B}$, respectively, we obtain for the Radon-Nikodym derivative inside (4)

[TABLE]

2.5 Path mutual information

We recognize that the path mutual information defined in (4) is just the expectation of the logarithm of (11), which reads

[TABLE]

Note that (12) can be simplified by realizing that the resulting (finite) sums over the reaction occurrences can be written as stochastic integrals with respect to the differential reaction counters $\mathrm{d}N_{k}(t)$. Before we state the final result, we take into account the fact that only the coupling reactions are affected by the marginalization. The propensities of all other reactions will remain unaffected when $\mathcal{F}^{AB}$ is replaced by $\mathcal{F}^{A}$ or $\mathcal{F}^{B}$, respectively. Those reactions will have two contributions in (12) with opposite sign and will thus cancel out. Taking this into account and performing some minor rearrangements yields

[TABLE]

where we recall that the symbol $\mathcal{R}_{A|B}$ denotes the coupling reactions affecting species $\mathrm{A}$ which are driven by $\mathrm{B}$, and vice versa for $\mathcal{R}_{B|A}$. The symmetric structure of (13) reflects the fact that information can be transferred both from $\mathcal{R}_{A}$ to $\mathcal{R}_{B}$ as well as $\mathcal{R}_{B}$ to $\mathcal{R}_{A}$. In order to numerically evaluate (13), we can use Gillespie’s stochastic simulation algorithm to simulate paths of the whole network $\mathcal{R}_{\SS}$ and for each of them, evaluate the integrals in (13) using the marginal propensities described in Section 2.3. Subsequently, the path-MI is obtained by averaging over all sample paths.

2.6 Connection to previous results

While our study focuses on the explicit calculation of the path-MI from (13), we want to show how our results can be related to existing theoretical work in the literature. Following [16, 17], we realize that the expectation of the Riemann integrals in (13) is equal to zero. This can be seen by exchanging the order of the expectation and time integration and realizing that $\mathbb{E}\bqty{\lambda_{k}^{AB}(s)-\lambda_{k}^{A}(s)}=0$ and $\mathbb{E}\bqty{\lambda_{k}^{AB}(s)-\lambda_{k}^{B}(s)}=0$ by definition of the marginal propensity. Additionally, we can use the Doob-Meyer decomposition theorem and expand the differential reaction counters in a predictable part and a martingale $\mathrm{d}N_{k}(s)=\lambda_{k}^{AB}(s)\mathrm{d}s+\mathrm{d}Q_{k}(s)$, where $\mathbb{E}\bqty{\mathrm{d}Q_{k}(s)\mid\mathcal{F}_{s}^{AB}}=0$. Thus, (13) can be reformulated as

[TABLE]

where $\mathcal{R}_{C}=\mathcal{R}_{A|B}\cup\mathcal{R}_{B|A}$ denotes the set of coupling reactions. Note that (14) resembles the point process mutual information given in [16] and analogous results can be found in [17, 19]. Similar considerations hold also for the transfer entropy [22], whose definition is related to (4).

3 Case Studies

In the following, we demonstrate the provided estimator of the path mutual information using three case studies. All simulations have been performed using the programming language julia [30]. The code is publicly available at https://github.com/zechnerlab/PathMI .

3.1 Protein expression network

We consider a simple two-stage model of gene expression given by the reaction network

[TABLE]

with $\mathrm{A}$ as mRNA, $\mathrm{B}$ as protein and $\gamma_{A}$, $\gamma_{B}$, $\delta_{A}$ and $\delta_{B}$ the rate constants associated to mRNA and protein synthesis and degradation. We want to quantify the path mutual information between the mRNA and protein paths. In this example the evolution of $\mathrm{A}$ is not driven by $\mathrm{B}$, thus $\mathcal{R}_{A|B}=\emptyset$. Therefore, in order to estimate $\mathbb{I}({A}_{0}^{t},{B}_{0}^{t})$, we only need to focus on reaction $\#3$ and be able to compute the marginal propensity $\lambda_{3}^{B}(t)=\gamma_{B}\hat{A}_{1}(t)$, where $\hat{A}_{1}(t)=\mathbb{E}\bqty{A_{t}|\mathcal{F}_{t}^{B}}$ denotes the first moment of the conditional distribution $\pi^{B}(a,t)$. Thus, (13) just corresponds to

[TABLE]

where $\mathrm{d}N_{3}(t)$ is the differential reaction counter of the coupling reaction. In order to evaluate $\hat{A}_{1}$, we instantiate the filtering equation from (9) for this particular example, i.e.,

[TABLE]

We can now obtain a stochastic differential equation for $\hat{A}_{1}(t)=\sum_{0}^{\infty}a\,\pi^{B}(a,t)$ by multiplying both sides of $(\ref{eq:filtering1})$ by $a$ and summing over all possible values of $a$. This yields

[TABLE]

However, since (18) depends on higher order moments, we adopt a Gamma closure at second order so that we can write the closed evolution of $\hat{A}_{2}=\mathbb{E}\bqty{A_{t}^{2}|\mathcal{F}_{t}^{B}}$ as

[TABLE]

Equations (18) and (19) can be solved alongside of the SSA simulations in order to evaluate (16). Finally, $\mathbb{I}({A}_{0}^{t},{B}_{0}^{t})$ can be estimated as a Monte Carlo over multiple SSA simulations. Fig. 1(a-b) shows simulation results for this example. The mRNA and protein levels were initialized to $\gamma_{A}/\delta_{A}$ and zero respectively. After a short initial transient, the path mutual information increases linearly with time, indicating that information is transferred from $\mathrm{A}$ to $\mathrm{B}$ at constant rate as soon as the system approaches its stationary state. Individual realizations are shown to illustrate their variability around the linear trend.

3.2 Transient induction of transcription

We again consider the reaction network introduced in (15) but in a different regime. In particular we assume that transcription is switched off at time zero such that $\gamma_{A}=0$. The initial pool of mRNA $A_{0}>0$ will then degrade as time increases and the system will ultimately converge to a zero steady state. Our goal is to quantify the information transmission between species $\mathrm{A}$ and $\mathrm{B}$ during this transient period. Note that the expressions (16), (17), (18) and (19) still apply, with the only difference that $\gamma_{A}=0$. In Fig. 1(c-d) we show the results for the initial conditions $A_{0}=100$ and $B_{0}=0$. The results show that the path-MI stops to increase after both the species are extinct. First, the mRNA gets depleted and consequently the coupling reaction $\#3$ stops firing. After this point, the $log$-jump contributions in (16) cease and the residual contribution only comes from the integral of the difference between $A_{s}=0$ and its causal estimate $\hat{A}_{1}(s)$ in the second line of (16).

3.3 A stochastic oscillatory system

In our last case study, we consider a stochastic predator-prey system inspired by the Lotka-Volterra model. The system is defined by the reaction network

[TABLE]

where the species $\mathrm{A}$ plays the role of the prey and $\mathrm{B}$ is the predator. A prey can duplicate with reproduction rate $\gamma_{A}$ and can be consumed by a predator with rate $\delta_{A}$. Instead, a predator can die with death rate $\delta_{B}$ and duplicate proportionally to the amount of prey, with rate $\gamma_{B}$.

Notice that in this example we have $\mathcal{R}_{A}=(\#1,\#2)$ and $\mathcal{R}_{B}=(\#3,\#4)$ and two coupling reactions, $\mathcal{R}_{A|B}=\#2$ and $\mathcal{R}_{B|A}=\#3$. It is worth mentioning that reaction $\#3$ satisfies the smooth coupling assumption because only $\mathrm{B}$ gets modified by it.

Note that, since in (20) the information transfer happens bi-directionally between $A$ and $B$, we require expressions for both the marginal propensities $\lambda_{2}^{A}(t)$ and $\lambda_{3}^{B}(t)$. Due to space considerations, we omit explicit expressions for $\pi^{A}$ and $\pi^{B}$ and the conditional moments, which have been obtained again under a Gamma closure. In Fig. 2(a) we show a realization of the system and the corresponding path-MI sample. The trajectories of prey and predator exhibit a characteristic oscillatory pattern. The path-MI sample grows significantly during the copy number peaks of prey and predator, because of the increased intensity of the coupling reactions. In this model it might happen that the preys go extinct in which case information transmission stops. As time increases, more and more realizations reach extinction, which causes the path-MI to saturate (Fig. 2(b)).

4 Conclusions

In this work we have presented a method to efficiently estimate the path mutual information between two chemical species whose dynamics evolve according to a stochastic reaction network. In order to derive the scheme, we employed our recently proposed marginal process framework, which allowed us to explicitly calculate the required Radon-Nikodym derivatives. We showed how the path-MI can be efficiently estimated by combining SSA with the solution of a stochastic filtering problem. We showed the efficacy of the method by applying it to three case studies with different dynamical properties. For the purposes of this work we restricted ourselves to a specific class of reaction systems, which consist of only two smoothly coupled species. In the future, we will extend the method to networks comprising more than two chemical species and arbitrary coupling reactions.

5 Acknowledgments

The authors were supported by core funding of the Max Planck Institute of Molecular Cell Biology and Genetics and the Center for Systems Biology Dresden and the Center for Advancing Electronics Dresden (cfaed).