# Path Weight Sampling: Exact Monte Carlo Computation of the Mutual
Information between Stochastic Trajectories

Manuel Reinhardt  AMOLF, Science Park 104, 1098 XG, Amsterdam, The Netherlands
Gašper Tkačik  Institute of Science and Technology Austria, 3400
Klosterneuburg, Austria  Pieter Rein ten Wolde<EMAIL_ADDRESS>AMOLF,
Science Park 104, 1098 XG, Amsterdam, The Netherlands

###### Abstract

Most natural and engineered information-processing systems transmit
information via signals that vary in time. Computing the information
transmission rate or the information encoded in the temporal characteristics
of these signals, requires the mutual information between the input and output
signals as a function of time, i.e. between the input and output trajectories.
Yet, this is notoriously difficult because of the high-dimensional nature of
the trajectory space, and all existing techniques require approximations. We
present an exact Monte Carlo technique called _Path Weight Sampling_ (PWS)
that, for the first time, makes it possible to compute the mutual information
between input and output trajectories for any stochastic system that is
described by a master equation. The principal idea is to use the master
equation to evaluate the exact conditional probability of an individual output
trajectory for a given input trajectory, and average this via Monte Carlo
sampling in trajectory space to obtain the mutual information. We present
three variants of PWS, which all generate the trajectories using the standard
stochastic simulation algorithm. While _Direct PWS_ is a brute-force method,
_Rosenbluth-Rosenbluth PWS_ exploits the analogy between signal trajectory
sampling and polymer sampling, and _Thermodynamic Integration PWS_ is based on
a reversible work calculation in trajectory space. PWS also makes it possible
to compute the mutual information between input and output trajectories for
systems with hidden internal states as well as systems with feedback from
output to input. Applying PWS to the bacterial chemotaxis system, consisting
of 182 coupled chemical reactions, demonstrates not only that the scheme is
highly efficient, but also that the number of receptor clusters is much
smaller than hitherto believed, while their size is much larger.

Complex Systems, Statistical Physics, Computational Physics

##  I Introduction

Quantifying information transmission is vital for understanding and designing
natural and engineered information-processing systems, ranging from
biochemical and neural networks, to electronic circuits and optical systems
[1, 2, 3]. Claude Shannon introduced the mutual information and the
information rate as the central measures of Information Theory more than 70
years ago [4]. These measures quantify the fidelity by which a noisy system
transmits information from its inputs to its outputs. Yet, computing these
quantities exactly remains notoriously difficult, if not impossible. This is
because the inputs and outputs are often not scalar values, but rather
temporal _trajectories_.

Most, if not all, information-processing systems transmit signal that vary in
time. The canonical measure for quantifying information transmission via time-
varying signals is the mutual information rate [4, 5, 6, 7]. It quantifies the
speed at which distinct messages are transmitted through the system, and it
depends not only on the accuracy of the input-output mapping but also on the
correlations within the input and output signals. Computing the mutual
information rate thus requires computing the mutual information between the
input and output trajectories, not between their signal values at given time
points. The rate at which this trajectory mutual information increases with
the trajectory duration in the long-time limit defines the mutual information
rate. In the absence of feedback this rate also equals the multi-step transfer
entropy [8, 9].

More generally, useful information is often contained in the temporal dynamics
of the signal. A prime example is bacterial chemotaxis, where the response
does not depend on the current ligand concentration, but rather on whether it
has changed in the recent past [10, 11]. Moreover, the information from the
input may be encoded in the temporal dynamics of the output [12, 13, 14, 15].
Quantifying information encoded in these temporal features of the signals
requires the mutual information not between two time points, i.e. the
instantaneous mutual information, but rather between input and output
trajectories [6].

Unfortunately, computing the mutual information between trajectories is
exceptionally difficult. The conventional approach requires non-parametric
distribution estimates of the input and output distributions, e.g. via
histograms of data obtained through simulations or experiments [16, 17, 18,
19, 20, 21]. These non-parametric distribution estimates are necessary because
the mutual information cannot generally be computed from summary statistics
like the mean or variance of the data alone. However, the high-dimensional
nature of trajectories makes it infeasible to obtain enough empirical data to
accurately estimate the required probability distributions. Moreover, this
approach requires the discretization of time, which becomes problematic when
the information is encoded in the precise timing of signal spikes, as, e.g.,
in neuronal systems [22]. Except for the simplest systems with a binary state
space [21], the conventional approach to estimate the mutual information via
histograms therefore cannot be transposed to trajectories.

Because there are currently no general schemes available to compute the mutual
information between trajectories exactly, approximate methods or simplified
models are typically used. While empirical distribution estimates can be
avoided by employing the K-nearest-neighbors entropy estimator [23, 24], this
method depends on a choice of metric in trajectory space and can become
unreliable for long trajectories [25]. Alternative, decoding-based information
estimates can be developed for trajectories [26], but merely provide a lower
bound of the mutual information, and it remains unclear how tight these lower
bounds are [27, 25, 28]. Analytical results are avaiable for simple systems
[29], and for linear systems that obey Gaussian statistics, the mutual
information between trajectories can be obtained from the covariance matrix
[6]. However, many information processing systems are complex and non-linear
such that the Gaussian approximation does not hold, and analytical solutions
do not exist. A more promising approach to estimate the trajectory mutual
information for chemical reaction networks has been developed by Duso and
Zechner [30] and generalized in Ref. [31]. However, the scheme relies on a
moment closure approximation and has so far only been applied to very simple
networks, seemingly being difficult to extend to complex systems.

Here, we present _Path Weight Sampling_ (PWS), an _exact_ technique to compute
the trajectory mutual information for any system described by a master
equation. Master equations are widely used to model chemical reaction networks
[32, 33, 34, 35], biological population growth [36, 37, 38], economic
processes [39, 40], and a large variety of other systems [41, 42], making our
scheme of interest to a broad class of problems.

PWS is an exact Monte Carlo scheme, in the sense that it provides an unbiased
statistical estimate of the trajectory mutual information. In PWS, the mutual
information is computed as the difference between the marginal output entropy
associated with the marginal distribution $\mathcal{P}[\bm{x}]$ of the output
trajectories $\bm{x}$, and the conditional output entropy associated with the
output distribution $\mathcal{P}[\bm{x}|\bm{s}]$ conditioned on the input
trajectory $\bm{s}$. Our scheme is inspired by the observation from Cepeda-
Humerez  _et al._ [25] that the path likelihood, i.e. the probability
$\mathcal{P}[\bm{x}|s]$, can be computed exactly from the master equation for
a _static_ input signal $s$. This makes it possible to compute the mutual
information between a discrete input and a time-varying output via a Monte
Carlo averaging procedure of the likelihoods, rather than from an empirical
estimate of the intractable high-dimensional probability distribution
functions. The scheme of Cepeda-Humerez  _et al._ [25] is however limited to
discrete input signals that do not vary in time. Here we show that the path
likelihood $\mathcal{P}[\bm{x}|\bm{s}]$ can also be computed for a dynamical
input _trajectory_ $\bm{s}$, which allows us to compute the conditional output
entropy also for time-varying inputs. While this solves the problem in part,
the marginal output entropy associated with $\mathcal{P}[\bm{x}]$ cannot be
computed with the approach of Cepeda-Humerez  _et al._ , and thus requires a
different scheme.

We show how, for time-varying input signals, the marginal probability
$\mathcal{P}[\bm{x}]$ can be obtained as a Monte Carlo average of
$\mathcal{P}[\bm{x}|\bm{s}]$ over a large number of input trajectories. How to
do this effectively is the crux of PWS. We then use the Monte Carlo estimate
for $\mathcal{P}[\bm{x}]$ to compute the marginal output entropy. We present
three variants of PWS, all of which compute the conditional entropy in the
same manner, but differ in the way this Monte Carlo averaging procedure for
computing the marginal probability $\mathcal{P}[\bm{x}]$ is carried out.

To compute $\mathcal{P}[\bm{x}]$, _Direct_ PWS (DPWS) performs a brute-force
average of the path likelihoods $\mathcal{P}[\bm{x}|\bm{s}]$ over the input
trajectories $\bm{s}$. While we show that this scheme works for simple
systems, the brute-force Monte Carlo averaging procedure becomes more
difficult for larger systems and exponentially harder for longer trajectories.

Our second and third variant of PWS are based on the realization that the
marginal probability $\mathcal{P}[\bm{x}]$ is akin to a partition function.
These schemes leverage techniques for computing free energies from statistical
physics. Specifically, the second scheme, _Rosenbluth-Rosenbluth PWS_ (RR-
PWS), exploits the observation that the computation of $\mathcal{P}[\bm{x}]$
is analogous to the calculation of the (excess) chemical potential of a
polymer, for which efficient methods have been developed [43, 44, 45]. The
third scheme, _Thermodynamic Integration PWS_ (TI-PWS), is based on the
classic free energy estimation technique of thermodynamic integration [46, 47,
48] in conjunction with a trajectory space MCMC sampler using ideas from
transition path sampling [49].

All three PWS variants make it possible to compute the mutual information
between trajectories exactly, without the need of time-discretizing the input
and output signals.

The DPWS method is presented in Section II, followed by a review of the
required concepts from the theory of Markov jump processes and master
equations. The other two PWS schemes are then presented as improvements of
DPWS in Section III.

In Section IV we show that, surprisingly, our PWS methods additionally make it
possible to compute the mutual information between input and output
trajectories of systems with hidden internal states. Hidden states correspond,
for example, to network components that merely relay, process or transform the
signal from the input to the output. Indeed, the downstream system typically
responds to the information that is encoded in this output, and not the other
internal system components. Most information processing systems contain such
hidden states, and generally we want to integrate out these latent network
components. In addition, we can generalize PWS to systems with feedback from
the output to the input as shown in Appendix C.

In Section V we apply PWS to two well-known model systems. The first is a
simple pair of coupled birth-death processes which allows us to test the
efficiency of the three PWS variants, as well as to compare the PWS results
with analytical results from the Gaussian approximation [6] and the technique
by Duso and Zechner [30]. Our second application concerns the bacterial
chemotaxis system, which is arguably the best characterized signaling system
in biology. Mattingly  _et al._ [50] recently argued that bacterial chemotaxis
in shallow gradients is information limited. Yet, to compute the information
rate from their experimental data they had to employ a Gaussian framework. PWS
makes it possible to asses the accuracy of this approximation. Our results
show that the Gaussian framework is accurate in the regime of shallow
concentration gradients as studied in Mattingly  _et al._ [50]. Yet, comparing
the PWS predictions of a model based on previous literature against their
experimental results reveals that the composition of the receptor array
differs from that hitherto believed.

##  II Monte Carlo Estimate of the Mutual Information

In this section we present the fundamental ideas of PWS. These ideas lie at
the heart of DPWS and also form the foundation of the other two more advanced
PWS variants which will be explained in subsequent sections.

###  II.1 Statement of the Problem

All information processing systems repeatedly take an input value $s$ and
produce a corresponding output $x$. Due to noise, the output produced for the
same input can be different every time, such that the system samples outputs
from the distribution $\mathrm{P}(x|s)$. In the information theoretic sense,
the device’s capabilities are fully specified by its output distributions for
all possible inputs. We consider the inputs as being distributed according to
a probability density $\mathrm{P}(s)$ such that the whole setup of signal and
device is completely described by the joint probability density
$\mathrm{P}(s,x)=\mathrm{P}(s)\,\mathrm{P}(x|s)$.

When the conditional output distributions $\mathrm{P}(x|s)$ overlap with each
other, information is lost because the input can not always be inferred
uniquely from the output (see Fig. 1). The remaining information that the
output carries about the signal on average is quantified by the mutual
information between input and output.

Figure 1: Schematic of information processing under the influence of noise.
Overlapping output distributions for different inputs lead to information
loss, because the input cannot always be uniquely inferred from the output.
The mutual information $\mathrm{I}(\mathcal{S},\mathcal{X})$ quantifies how
much information the observation of the output typically retains about the
input signal.

Mathematically, the mutual information between a random variable
$\mathcal{S}$, representing the input, and a second random variable
$\mathcal{X}$, representing the output, is defined as

$\mathrm{I}(\mathcal{S},\mathcal{X})=\iint\mathrm{d}s\,\mathrm{d}x\
\mathrm{P}(s,x)\ln\frac{\mathrm{P}(s,x)}{\mathrm{P}(s)\mathrm{P}(x)}\,,$ (1)

where the marginal output distribution is given by
$\mathrm{P}(x)=\int\mathrm{d}s\ \mathrm{P}(s,x)$. The quantity
$\mathrm{I}(\mathcal{S},\mathcal{X})$ as defined above is a non-negative real
number, representing the mutual information between $\mathcal{S}$ and
$\mathcal{X}$ in nats. The integrals in Eq. 1 run over all possible
realizations of the random variables $\mathcal{S}$ and $\mathcal{X}$. In our
case, $\mathcal{S}$ and $\mathcal{X}$ represent stochastic trajectories and so
the integrals become path integrals.

In general, the mutual information can be decomposed into two terms, a
conditional and marginal entropy. Due to the symmetry of Eq. 1 with respect to
exchange of $\mathcal{S}$ and $\mathcal{X}$, this decomposition can be written
as

$\mathrm{I}(\mathcal{S},\mathcal{X})=\mathrm{H}(\mathcal{S})-\mathrm{H}(\mathcal{S}|\mathcal{X})=\mathrm{H}(\mathcal{X})-\mathrm{H}(\mathcal{X}|\mathcal{S})\,.$
(2)

The (marginal) input entropy $\mathrm{H}(\mathcal{S})$ represents the total
uncertainty about the input, and the conditional input entropy
$\mathrm{H}(\mathcal{S}|\mathcal{X})$ describes the remaining uncertainty of
the input after having observed the output. Thus, the mutual information
$\mathrm{I}(\mathcal{S},\mathcal{X})=\mathrm{H}(\mathcal{S})-\mathrm{H}(\mathcal{S}|\mathcal{X})$
naturally quantifies the reduction in uncertainty about the input through the
observation of the output.

When analyzing data from experiments or simulations however, the mutual
information is generally estimated via
$\mathrm{I}(\mathcal{S},\mathcal{X})=\mathrm{H}(\mathcal{X})-\mathrm{H}(\mathcal{X}|\mathcal{S})$.
This is because simulation or experimental data generally provide information
about the distribution of outputs for a given input, rather than vice versa.
The accessible entropies are thus the marginal output entropy
$\mathrm{H}(\mathcal{X})$ and the conditional output entropy
$\mathrm{H}(\mathcal{X}|\mathcal{S})$, which are defined as

$\displaystyle\mathrm{H}(\mathcal{X})$ $\displaystyle=-\int\mathrm{d}x\
\mathrm{P}(x)\ln\mathrm{P}(x)$ (3)
$\displaystyle\mathrm{H}(\mathcal{X}|\mathcal{S})$
$\displaystyle=-\int\mathrm{d}s\ \mathrm{P}(s)\int\mathrm{d}x\
\mathrm{P}(x|s)\ln\mathrm{P}(x|s)\,.$ (4)

The conventional way of computing the mutual information involves generating
many samples to obtain empirical distribution estimates for $\mathrm{P}(x|s)$
and $\mathrm{P}(x)$ via histograms. However, the number of samples needs to be
substantially larger than the number of histogram bins to reduce the noise in
the bin counts. Obtaining enough samples is effectively impossible for high-
dimensional data, like signal trajectories. Moreover, any nonzero bin size
leads to a systematic bias in the entropy estimates, even in one dimension
[17]. These limitations of the conventional method make it impractical for
high-dimensional data, highlighting the need for alternative approaches to
accurately compute mutual information for trajectories.

###  II.2  Direct PWS

Figure 2: The PWS scheme to compute the mutual information between
trajectories in 4 steps. 1. Generate $N$ input-output pairs from
$\mathcal{P}[\bm{s},\bm{x}]$. 2. For each input-output pair compute the
trajectory likelihood $\mathcal{P}[\bm{x_{i}}|\bm{s_{i}}]$ using Eq. 14. 3.
Compute $\mathcal{P}[\bm{x}_{i}]$ for every output. This step differentiates
the different variants of PWS from each other. Direct PWS is presented in
Section II.2, whereas RR-PWS and TI-PWS are presented in Sections III.2 and
III.3. 4. Using the likelihoods and the marginal probabilities from the
previous steps we can estimate the mutual information using Eq. 10.

The central idea of PWS is to compute probability densities for trajectories
exactly, sidestepping the problem having to estimate them via histograms. We
exploit that for systems described by a master equation, the conditional
probability of an output trajectory for a given input trajectory can be
computed analytically. With this insight we can derive a procedure to compute
the mutual information. Specifically, we will show that

  * •

for a system described by a master equation, the trajectory likelihood
$\mathcal{P}[\bm{x}|\bm{s}]$ is a quantity that can be computed on the fly in
a stochastic simulation;

  * •

input trajectories can be generated from $\mathcal{P}[\bm{s}]$, output
trajectories for a given input $\bm{s}$ can be generated according to
$\mathcal{P}[\bm{x}|\bm{s}]$ using standard SSA (Gillespie) simulations;

  * •

by combining the two ideas above, we can derive a direct Monte Carlo estimate
for the mutual information $\mathrm{I}(\mathcal{S},\mathcal{X})$, as
illustrated in Fig. 2.

Note that we denote trajectories by bold symbols to distinguish them from
scalar quantities.

Our technique is conceptually straightforward. Using Monte Carlo simulations
we can compute averages over the configuration space of trajectories. Suppose
we have a function $f[\bm{z}]$ that takes a trajectory $\bm{z}$ and produces a
scalar value. The mean of $f[\bm{z}]$ with respect to the trajectory
distribution $\mathcal{P}[\bm{z}]$ is then

$\langle
f[\bm{z}]\rangle_{\mathcal{P}[\bm{z}]}\equiv\int\mathcal{D}[\bm{z}]\,\mathcal{P}[\bm{z}]f(\bm{z})\,.$
(5)

We write $\int\mathcal{D}[\bm{z}]$ to denote a path integral over all possible
trajectories of a given duration. We estimate $\langle
f[\bm{z}]\rangle_{\mathcal{P}[\bm{z}]}$, by generating a large number of
trajectories $\bm{z}_{1},\ldots,\bm{z}_{N}$ from $\mathcal{P}[\bm{z}]$ and
evaluating the corresponding Monte Carlo average

$\hat{f}_{N}=\frac{1}{N}\sum^{N}_{i=1}f(\bm{z}_{i})$ (6)

which converges to the true mean in the limit $N\to\infty$.

Specifically, we want to estimate the conditional and the marginal entropy to
compute the mutual information. Let us imagine that we generate $N$ input
trajectories $\bm{s}_{1},\ldots,\bm{s}_{N}$ from the distribution
$\mathcal{P}[\bm{s}]$. Next, for every input $\bm{s}_{i}$, we generate a set
of $K$ outputs $\bm{x}_{i,1},\ldots,\bm{x}_{i,K}$ from
$\mathcal{P}[\bm{x}|\bm{s}_{i}]$. Then, the Monte Carlo estimate for the
conditional entropy is

$\displaystyle\mathrm{H}(\mathcal{X}|\mathcal{S})$
$\displaystyle=-\int\mathcal{D}[\bm{s}]\
\mathcal{P}[\bm{s}]\int\mathcal{D}[\bm{x}]\
\mathcal{P}[\bm{x}|\bm{s}]\ln\mathcal{P}[\bm{x}|\bm{s}]$ (7)
$\displaystyle=-\left\langle\left\langle\ln\mathcal{P}[\bm{x}|\bm{s}]\right\rangle_{\mathcal{P}[\bm{x}|\bm{s}]}\right\rangle_{\mathcal{P}[\bm{s}]}$
$\displaystyle\approx-\frac{1}{N}\sum^{N}_{i=1}\frac{1}{K}\sum^{K}_{j=1}\ln\mathcal{P}[\bm{x}_{i,j}|\bm{s}_{i}]\,.$

Secondly, for a given output $\bm{x}$ we generate $M$ inputs
$\bm{s}^{\prime}_{1},\ldots,\bm{s}^{\prime}_{M}$ according to
$\mathcal{P}[\bm{s}]$, then we can obtain a Monte Carlo estimate for the
marginal probability of the output trajectory $\mathcal{P}[\bm{x}]$:

$\displaystyle\mathcal{P}[\bm{x}]$ $\displaystyle=\int\mathcal{D}[\bm{s}]\
\mathcal{P}[\bm{s}]\mathcal{P}[\bm{x}|\bm{s}]$ (8)
$\displaystyle=\left\langle\mathcal{P}[\bm{x}|\bm{s}]\right\rangle_{\mathcal{P}[\bm{s}]}$
$\displaystyle\approx\frac{1}{M}\sum^{M}_{j=1}\mathcal{P}[\bm{x}|\bm{s}^{\prime}_{j}]\,.$

The estimate for the marginal entropy is then given by

$\displaystyle\mathrm{H}(\mathcal{X})$
$\displaystyle=-\int\mathcal{D}[\bm{x}]\
\mathcal{P}[\bm{x}]\ln\mathcal{P}[\bm{x}]$ (9)
$\displaystyle=-\left\langle\ln\mathcal{P}[\bm{x}]\right\rangle_{\mathcal{P}[\bm{x}]}$
$\displaystyle\approx-\frac{1}{N}\sum^{N}_{i=1}\ln\mathcal{P}[\bm{x}_{i}]$
$\displaystyle\approx-\frac{1}{N}\sum^{N}_{i=1}\ln\left[\frac{1}{M}\sum^{M}_{j=1}\mathcal{P}[\bm{x}_{i}|\bm{s}^{\prime}_{i,j}]\right]\,.$

In the last step we inserted the result from Eq. 8. In this estimate, the
trajectories $\bm{x}_{1},\ldots,\bm{x}_{N}$ are sampled from
$\mathcal{P}[\bm{x}]$, i.e., by first sampling from $\mathcal{P}[\bm{s}]$ and
then from $\mathcal{P}[\bm{x}|\bm{s}]$. Finally, the mutual information is
obtained by taking the entropy difference, i.e.,
$\mathrm{I}(\mathcal{S},\mathcal{X})=\mathrm{H}(\mathcal{X})-\mathrm{H}(\mathcal{X}|\mathcal{S})$.

While this is the main idea behind PWS, it is computationally advantageous to
change the order of operations in the estimate. Specifically, computing the
difference of two averages, leads to large statistical errors. We can obtain
an improved estimate by reformulating the mutual information as a single
average of differences:

$\displaystyle\mathrm{I}(\mathcal{S},\mathcal{X})$
$\displaystyle=\int\mathcal{D}[\bm{s}]\int\mathcal{D}[\bm{x}]\
\mathcal{P}[\bm{s},\bm{x}]\ln\frac{\mathcal{P}[\bm{x}|\bm{s}]}{\mathcal{P}[\bm{x}]}$
(10)
$\displaystyle=\left\langle\ln\mathcal{P}[\bm{x}|\bm{s}]-\ln\mathcal{P}[\bm{x}]\right\rangle_{\mathcal{P}{[\bm{s},\bm{x}]}}\,.$

This equation applies to all variants of PWS. They differ, however, in the way
$\mathcal{P}[\bm{x}]$ is computed. In the brute-force version of PWS, called
_Direct PWS_ (DPWS), we use Eq. 8 to evaluate the marginal probability
$\mathcal{P}[\bm{x}]$. DPWS indeed involves two nested Monte Carlo
computations, in which $N$ pairs $(\bm{s}_{i},\bm{x}_{i})$ are generated, and
for each output $\bm{x}_{i}$, $M$ input trajectories $\\{\bm{s}\\}$ are
generated from scratch to compute $\mathcal{P}[\bm{x}]$. In Section III below,
we will present two additional variants of PWS where the brute-force estimate
of the marginal probability $\mathcal{P}[\bm{x}]$ is replaced by more
elaborate schemes. That said, DPWS is a conceptually simple, straightforward
to implement, and exact scheme to compute the mutual information.

Having explained the core ideas of our technique above, we will continue this
section with a review of the necessary concepts of master equations to
implement PWS. First, in Section II.3, we derive the formula for the
conditional probability $\mathcal{P}[\bm{x}|\bm{s}]$ which lies at the heart
of our technique. In Sections II.3 and II.4, we discuss how trajectories are
generated according to $\mathcal{P}[\bm{x}|\bm{s}]$ and $\mathcal{P}[\bm{s}]$,
which are the remaining ingredients required for using DPWS. Then, in Section
III, we will present the two other variants of PWS that improve on DPWS.

###  II.3  Driven Markov Jump Process

Throughout this article we consider systems that can be modeled by a master
equation and are being driven by a stochastic signal. The master equation
specifies the time evolution of the conditional probability distribution
$\mathrm{P}(x,t|x_{0},t_{0})$ which is the probability for the process to
reach the discrete state $x\in\Omega$ at time $t$, given that it was at state
$x_{0}\in\Omega$ at the previous time $t_{0}$. The state space $\Omega$ is
multi-dimensional if the system is made up of multiple components and
therefore $x$ and $x_{0}$ can be vectors rather than scalar values. Denoting
the transition rate at time $t$ from state $x$ to another state
$x^{\prime}\neq x$ by $w_{t}(x^{\prime},x)$, the master equation reads

$\frac{\partial\mathrm{P}(x,t)}{\partial
t}=\sum_{\begin{subarray}{c}x^{\prime}\in\Omega\\\ x^{\prime}\neq
x\end{subarray}}[w_{t}(x,x^{\prime})\mathrm{P}(x^{\prime},t)-w_{t}(x^{\prime},x)\mathrm{P}(x,t)]\,,$
(11)

where, for brevity, we suppress the dependence on the initial condition, i.e.,
$\mathrm{P}(x,t)=\mathrm{P}(x,t|x_{0},t_{0})$. By defining
$Q_{t}(x^{\prime},x)=w_{t}(x^{\prime},x)$ for $x\neq x^{\prime}$ and
$Q_{t}(x,x)=-\sum_{x^{\prime}\in\Omega\smallsetminus\\{x\\}}w_{t}(x^{\prime},x)$
the master equation simplifies to

$\frac{\partial\mathrm{P}(x,t)}{\partial
t}=\sum_{x^{\prime}\in\Omega}Q_{t}(x,x^{\prime})\mathrm{P}(x^{\prime},t)\,.$
(12)

Note that by definition the diagonal matrix element $Q_{t}(x,x)$ is the
negative exit rate from state $x$, i.e. the total rate at which probability
flows away from state $x$.

Using the master equation we can compute the probability of any trajectory. A
trajectory $\bm{x}$ is defined by a list of jump times $t_{1},\ldots,t_{n-1}$,
together with a sequence of system states $x_{0},\ldots,x_{n-1}$. The
trajectory starts at time $t_{0}$ in state $x_{0}$ and ends at time $t_{n}$ in
state $x_{n-1}$, such that its duration is $T=t_{n}-t_{0}$. At each time
$t_{i}$ (for $i=1,\ldots,n-1$) the trajectory describes an instantaneous jump
$x_{i-1}\rightarrow x_{i}$. The probability density of $\bm{x}$ is

$\displaystyle\mathcal{P}[\bm{x}]$
$\displaystyle=\mathrm{P}(x_{0})\times\left(\prod^{n-1}_{i=1}Q_{t_{i}}\left(x_{i},x_{i-1}\right)\right)$
(13)
$\displaystyle\quad\times\left(\prod^{n}_{i=1}\exp\int\limits^{t_{i}}_{t_{i-1}}\mathrm{d}t\
Q_{t}(x_{i-1},x_{i-1})\right),$

a product of the probability of the initial state $\mathrm{P}(x_{0})$, the
rates of the $n-1$ transitions $Q_{t_{i}}\left(x_{i},x_{i-1}\right)$, and the
survival probabilities for the waiting times between jumps, given by
$\exp\int^{t_{i}}_{t_{i-1}}\mathrm{d}t\ Q_{t}(x_{i-1},x_{i-1})$ for
$i=1,\ldots,n$.

####  II.3.1 Computing the Likelihood $\mathcal{P}[x|s]$

To compute likelihood or conditional probability $\mathcal{P}[\bm{x}|\bm{s}]$
of an output trajectory $\bm{x}$ for a given input trajectory $\bm{s}$, we
note that the input determines the time-dependent stochastic dynamics of the
jump process. Indeed, the transition rates at time $t$, given by
$Q_{t}(x^{\prime},x;{\bm{s}})$, depend explicitly on the input $s(t)$ at time
$t$ and may even depend on the entire history of ${\bm{s}}$ prior to $t$.

In the common case that every input trajectory $\bm{s}$ leads to a unique
transition rate matrix $Q_{t}(x^{\prime},x;\bm{s})$, i.e. the map
$\bm{s}\mapsto Q_{t}(\cdot,\cdot;\bm{s})$ is injective, the likelihood is
directly given by Eq. 13:

$\displaystyle\mathcal{P}[\bm{x}|\bm{s}]$
$\displaystyle=\mathrm{P}(x_{0}|s_{0})\times\left(\prod^{n-1}_{i=1}Q_{t_{i}}\left(x_{i},x_{i-1};\bm{s}\right)\right)$
(14)
$\displaystyle\quad\times\left(\prod^{n}_{i=1}\exp\int\limits^{t_{i}}_{t_{i-1}}\mathrm{d}t\
Q_{t}(x_{i-1},x_{i-1};\bm{s})\right)$

where $\mathrm{P}(x_{0}|s_{0})$ is the probability of the initial state
$x_{0}$ of the output given the initial state of the input $s_{0}=s(t_{0})$.

The evaluation of the trajectory likelihood is at the heart of our Monte Carlo
scheme. However, numerically computing a large product like Eq. 14 very
quickly reaches the limits of floating point arithmetic since the result is
often either too large or too close to zero to be representable as a floating
point number. Thus, to avoid numerical issues, it is vital to perform the
computations in log-space, i.e. to compute

$\ln\mathcal{P}[\bm{x}|\bm{s}]=\ln\mathrm{P}(x_{0}|s_{0})+\int^{T}_{t_{0}}\mathrm{d}t\
\mathcal{L}_{t}[\bm{s},\bm{x}]$ (15)

where

$\displaystyle\mathcal{L}_{t}[\bm{s},\bm{x}]$
$\displaystyle=Q_{t}(x(t),x(t);\bm{s})$ (16)
$\displaystyle\quad+\sum^{n-1}_{i=1}\delta(t-t_{i})\ln
Q_{t}(x_{i},x_{i-1};\bm{s})\,.$

The computation of the log-likelihood $\ln\mathcal{P}[\bm{x}|\bm{s}]$ for
given trajectories $\bm{s}$ and $\bm{x}$ according to Eqs. 15 and 16 proceeds
as follows:

  * •

At the start of the trajectory we compute the log-probability of the initial
condition $\ln\mathrm{P}(x_{0}|s_{0})$,

  * •

for every jump $x_{i-1}\rightarrow x_{i}$ in $\bm{x}$ compute the log jump
propensity $\ln Q_{t_{i}}(x_{i},x_{i-1};\bm{s})$, and

  * •

for every interval $(t_{i-1},t_{i})$ of constant output value $x(t)=x_{i-1}$
between two jumps of $\bm{x}$, we compute $\int^{t_{i}}_{t_{i-1}}\mathrm{d}t\
Q_{t}(x_{i-1},x_{i-1};\bm{s})$. This integral can be performed using standard
numerical methods such as the trapezoidal rule, which is also exact if
$Q_{t}(x(t),x(t);\bm{s})$ is a piecewise linear function of $t$ as in our
examples in Section V.

The sum of the three contributions above yields the exact log-likelihood
$\ln\mathcal{P}[\bm{x}|\bm{s}]$ as given in Eq. 15.

Thus, notably, the algorithm to compute the log-likelihood
$\ln\mathcal{P}[\bm{x}|\bm{s}]$ is both efficient and straightforward to
implement, being closely related to the standard Gillespie algorithm. The only
quantity in Eq. 15 that cannot be directly obtained from the master equation
is the log-probability of the initial state, $\ln\mathrm{P}(x_{0}|s_{0})$.

Our scheme can be applied to any system with a well-defined (non-equilibrium)
initial distribution $\mathrm{P}(s_{0},x_{0})$ as specified by, e.g. the
experimental setup. Most commonly though, one is interested in studying
information transmission for systems in steady state. Then, the initial
condition $\mathrm{P}(s_{0},x_{0})$ is the stationary distribution of the
Markov process. Depending on the complexity of the system, this distribution
can be found either analytically from the master equation [51, 52] (possibly
using simplifying approximations [53, 54]), or computationally from stochastic
simulations [55].

####  II.3.2 Sampling from $\mathcal{P}[x|s]$

Standard kinetic Monte Carlo simulations naturally produce exact samples of
the probability distribution $\mathcal{P}[\bm{x}|\bm{s}]$ as defined in Eq.
14. That is, for any signal trajectory $\bm{s}$ and initial state $x_{0}$
drawn from $\mathrm{P}(x_{0}|s_{0})$ we can use the Stochastic Simulation
Algorithm (SSA) or variants thereof to generate a corresponding trajectory
$\bm{x}$. The SSA propagates the initial condition $x_{0},t_{0}$ forward in
time according to the transition rate matrix $Q_{t}(\cdot;\bm{s})$. In the
standard Direct SSA algorithm [55] this is done by alternatingly sampling the
waiting time until the next transition, and then selecting the actual
transition.

The transition rates $Q_{t}(x^{\prime},x;\bm{s})$ of a driven master equation
are necessarily time-dependent since they include the coupling of the jump
process to the input trajectory $\bm{s}$, which itself varies in time. While
most treatments of the SSA assume that the transition rates are constant in
time, this restriction is easily lifted. Consider step $i$ of the Direct SSA
which generates the next transition time $t_{i+1}=t_{i}+\Delta t_{i}$. For
time-varying transition rates the distribution of the stochastic waiting time
$\Delta t_{i}$ is characterized by the survival function

$S_{i}(\tau)=\mathrm{P}(\Delta
t_{i}>\tau)=\exp\int^{t_{i}+\tau}_{t_{i}}\mathrm{d}t\
Q_{t}(x_{i},x_{i};\bm{s})\,.$ (17)

The waiting time can be sampled using inverse transform sampling, i.e. by
generating a uniformly distributed random number $u\in[0,1]$ and computing the
waiting time using the inverse survival function $\Delta t_{i}=S^{-1}_{i}(u)$.
Numerically, computing the inverse of the survival function requires solving
the equation

$\ln u=\int^{t_{i}+\Delta t_{i}}_{t_{i}}\mathrm{d}t\
Q_{t}(x_{i},x_{i};\bm{s})$ (18)

for the waiting time $\Delta t_{i}$. Depending on the complexity of
$Q_{t}(x_{i},x_{i}|\bm{s})$, this equation can either be solved analytically
or numerically, e.g. using Newton’s method. Hence, this method to generate
stochastic trajectories is only truly exact if we can solve Eq. 18
analytically, as in the example in Section V.1. Additionally, in some cases
more efficient variants of the SSA with time dependent rates could be used
[56, 57].

###  II.4 Input Statistics

For our mutual information estimate, we need to be able to draw samples from
the input distribution $\mathcal{P}[\bm{s}]$. Our algorithm poses no
restrictions on $\mathcal{P}[\bm{s}]$ other than the possibility to generate
sample trajectories.

For example, the input signal may be described by a continuous-time jump
process, as in Section V.1. One benefit is that it is possible to generate
exact realizations of such a process (using the SSA) and to exactly compute
the likelihood $\mathcal{P}[\bm{x}|\bm{s}]$ using Eq. 15. Specifically, the
likelihood can be exactly evaluated because the transition rates
$Q_{t}(\cdot,\cdot;\bm{s})$ for any input trajectory $\bm{s}$, while time-
dependent, are _piece-wise constant_. This implies that the integral in Eq. 15
can be evaluated analytically without approximations. Similarly, for piece-
wise constant transition rates, the inverse function of Eq. 18 can be
evaluated directly such that we can sample exact trajectories from the driven
jump process. As a result, when both input and output are described by a
master equation, PWS is a completely exact Monte Carlo scheme to compute the
mutual information.

However, the techniques described here do _not_ require the input signal
$\bm{s}$ to be described by a continuous-time jump process, or even to be
Markovian. The input signal can be any stochastic process for which
trajectories can be generated numerically. This includes continuous stochastic
processes that are found as solutions to stochastic differential equations
[58]. The application in Section V.2 provides an example.

##  III  Variants of PWS

The DPWS scheme presented in the previous section makes it possible to compute
the mutual information between trajectories of a stochastic input process and
the output of a Markov jump process. However, the number of possible
trajectories increases exponentially with trajectory length, leading to a
corresponding increase in the variance of the DPWS estimate. Hence, for long
trajectories the DPWS estimate may prove to be computationally too expensive.
To address this issue, we describe two improved variants of PWS in this
section, both based on free-energy estimators from statistical physics.

###  III.1 Marginalization Integrals in Trajectory Space

The computationally most expensive part of our scheme in Section II.2 is the
evaluation of the marginalization integral
$\mathcal{P}[\bm{x}_{i}]=\int\mathcal{D}[\bm{s}]\mathcal{P}[\bm{s},\bm{x}_{i}]$
which needs to be performed for every sample $\bm{x}_{1},\ldots,\bm{x}_{N}$.
Consequently the computational efficiency of this marginalization is essential
for the overall performance.

Marginalization is a general term to denote an operation where one or more
variables are integrated out of a joint probability distribution. For
instance, we obtain the marginal probability distribution
$\mathcal{P}[\bm{x}]$ from $\mathcal{P}[\bm{s},\bm{x}]$ by computing the
integral

$\mathcal{P}[\bm{x}]=\int\mathcal{D}[\bm{s}]\
\mathcal{P}[\bm{s},\bm{x}]=\int\mathcal{D}[\bm{s}]\
\mathcal{P}[\bm{s}]\mathcal{P}[\bm{x}|\bm{s}]\,.$ (19)

In DPWS, we use Eq. 8 to compute $\mathcal{P}[\bm{x}]$ which involves
generating independent input trajectories from $\mathcal{P}[\bm{s}]$. However,
this this is not the optimal Monte Carlo technique to perform the
marginalization. The generated input trajectories are independent from the
output trajectory $\bm{x}$. Thus, we ignore the causal connection between
$\bm{s}$ and $\bm{x}$, and we typically end up sampling trajectories
$\bm{s}^{\star}$ whose likelihoods $\mathcal{P}[\bm{x}|\bm{s}^{\star}]$ are
very small. Then, most sampled trajectories have small integral weights, and
only very few samples provide a significant contribution to the average. The
variance of the result is then very large because the effective sample size is
much smaller than the total sample size. The use of $\mathcal{P}[\bm{s}]$ as
the sampling distribution is thus only practical in cases where the dependence
of the output on the input is not too strong. It follows, perhaps
paradoxically, that this sampling scheme works best when the mutual
information is not too large 111Indeed, the mutual information
$\mathrm{I}(\mathcal{S},\mathcal{X})$ precisely quantifies how strong the
statistical dependence is between the trajectory-valued random variables
$\mathcal{S}$ and $\mathcal{X}$. From its definition
$\mathrm{I}(\mathcal{S},\mathcal{X})=\mathrm{H}(\mathcal{S})-\mathrm{H}(\mathcal{S}|\mathcal{X})$
we can understand more clearly how this affects the efficiency of the Monte
Carlo estimate. Roughly speaking, $\mathrm{H}(\mathcal{S})$ is related to the
number of distinct trajectories $\bm{s}$ that can arise from the dynamics
given by $\mathcal{P}[\bm{s}]$, while $\mathrm{H}(\mathcal{S}|\mathcal{X})$ is
related to the number of distinct trajectories $\bm{s}$ that could have lead
to a specific output $\bm{x}$, on average. Therefore, if the mutual
information is very large, the difference between these two numbers is very
large, and consequently the number of overall distinct trajectories is much
larger than the number of distinct trajectories compatible with output
$\bm{x}$. Now, if we generate trajectories according to the dynamics given by
$\mathcal{P}[\bm{s}]$, with overwhelming probability we generate a trajectory
$\bm{s}$ which is not compatible with the output trajectory $\bm{x}$, and
therefore $\mathcal{P}[\bm{x}|\bm{s}]\approx 0$. Hence, the effective number
of samples $M_{\text{eff}}$ is much smaller than the actual number of
generated trajectories $M$, i.e. $M_{\text{eff}}\ll M$. We therefore only
expect the estimate in Eq. 8 to be reliable when computing the mutual
information for systems where it is not too high. Thus, strikingly, the
difficulty of computing the mutual information is proportional to the
magnitude of the mutual information itself..

This is a well known Monte Carlo sampling problem and a large number of
techniques have been developed to solve it. The two variants of our scheme,
RR-PWS and TI-PWS, both make use of ideas from statistical physics for the
efficient computation of free energies.

$\mathcal{P}[\bm{s},\bm{x}]$ | $e^{-\mathcal{U}[\bm{s},\bm{x}]}$
---|---
$\mathcal{P}[\bm{s}]$ | $\frac{1}{\mathcal{Z}_{0}[\bm{x}]}e^{-\mathcal{U}_{0}[\bm{s}]}$
$\mathcal{P}[\bm{s}|\bm{x}]$ | $\frac{1}{\mathcal{Z}[\bm{x}]}e^{-\mathcal{U}[\bm{s},\bm{x}]}$
$1$ | $\mathcal{Z}_{0}[\bm{x}]$
$\mathcal{P}[\bm{x}]$ | $\mathcal{Z}[\bm{x}]$
$\mathcal{P}[\bm{x}|\bm{s}]$ | $e^{-\Delta\mathcal{U}[\bm{s},\bm{x}]}$
Table 1: Translation to the notation of statistical physics. The definitions
of $\mathcal{U}$ and $\mathcal{U}_{0}$ that are used here are given in Eqs. 24
and 20.

To understand how we can make use of these ideas to compute the marginal
probability $\mathcal{P}[\bm{x}]$, it is convenient to rephrase the
marginalization integral in Eq. 19 in the language of statistical physics. In
this language, $\mathcal{P}[\bm{x}]$ corresponds to the normalization
constant, or partition function, of a Boltzmann distribution for the potential

$\mathcal{U}[\bm{s},\bm{x}]=-\ln\mathcal{P}[\bm{s},\bm{x}]\,.$ (20)

In Eq. 20, we interpret $\bm{s}$ as a variable in the configuration space
whereas $\bm{x}$ is an auxiliary variable, i.e. a parameter. Note that both
$\bm{s}$ and $\bm{x}$ still represent trajectories. For this potential, the
partition function is given by

$\mathcal{Z}[\bm{x}]=\int\mathcal{D}[\bm{s}]\;e^{-\mathcal{U}[\bm{s},\bm{x}]}\,.$
(21)

The integral only runs over the configuration space, i.e. we integrate only
with respect to $\bm{s}$ but not $\bm{x}$, which remains a parameter of the
partition function. The partition function is precisely equal to the marginal
probability of the output, i.e. $\mathcal{Z}[\bm{x}]=\mathcal{P}[\bm{x}]$, as
can be verified by inserting the expression for the
$\mathcal{U}[\bm{s},\bm{x}]$. Further, the free energy is given by

$\mathcal{F}[\bm{x}]=-\ln\mathcal{Z}[\bm{x}]=-\ln\mathcal{P}[\bm{x}]$ (22)

which shows that the computation of the free energy of the trajectory ensemble
corresponding to $\mathcal{U}[\bm{s},\bm{x}]$ is equivalent to the computation
of (the logarithm of) the marginal probability $\mathcal{P}[\bm{x}]$.

Note that above we omitted any factors of $k_{\mathrm{B}}T$ since temperature
is irrelevant here. Also note that while the distribution
$\exp(-\mathcal{U}[\bm{s},\bm{x}])$ looks like the equilibrium distribution of
a canonical ensemble from statistical mechanics, this does not imply that we
can only study systems in thermal equilibrium. Indeed, PWS is used to study
information transmission in systems driven out of equilibrium by the input
signal. Thus, the notation introduced in this section is nothing else but a
mathematical reformulation of the marginalization integral to make the analogy
to statistical physics apparent and we assign no additional meaning of the
potentials and free energies introduced here.

In statistical physics it is well known that the free energy cannot be
directly measured from a simulation. Instead, one estimates the free-energy
difference

$\Delta\mathcal{F}[\bm{x}]=\mathcal{F}[\bm{x}]-\mathcal{F}_{0}[\bm{x}]=-\ln\frac{\mathcal{Z}[\bm{x}]}{\mathcal{Z}_{0}[\bm{x}]}$
(23)

between the system and a reference system with known free energy
$\mathcal{F}_{0}[\bm{x}]$. The reference system is described by the potential
$\mathcal{U}_{0}[\bm{s},\bm{x}]$ with the corresponding partition function
$\mathcal{Z}_{0}[\bm{x}]$. In our case, a natural choice of reference
potential is

$\mathcal{U}_{0}[\bm{s},\bm{x}]=-\ln\mathcal{P}[\bm{s}]$ (24)

with the corresponding partition function

$\mathcal{Z}_{0}[\bm{x}]=\int\mathcal{D}[\bm{s}]\mathcal{P}[\bm{s}]=1\,.$ (25)

This means that since $\mathcal{P}[\bm{s}]$ is a normalized probability
density function, the reference free energy is zero
($\mathcal{F}_{0}[\bm{x}]=-\ln\mathcal{Z}_{0}[\bm{x}]=0$). Hence, for the
above choice of reference system, the free-energy difference is

$\Delta\mathcal{F}[\bm{x}]=\mathcal{F}[\bm{x}]=-\ln\mathcal{P}[\bm{x}]\,.$
(26)

Note that in our case the reference potential
$\mathcal{U}_{0}[\bm{s},\bm{x}]=-\ln\mathcal{P}[\bm{s}]$ does not depend on
the output trajectory $\bm{x}$, i.e.
$\mathcal{U}_{0}[\bm{s},\bm{x}]\equiv\mathcal{U}_{0}[\bm{s}]$. It describes a
_non-interacting_ version of our input-output system where the input
trajectories evolve independently of the fixed output trajectory $\bm{x}$.

What is the interaction between the output $\bm{x}$ and the input trajectory
ensemble? We define the interaction potential
$\Delta\mathcal{U}[\bm{s},\bm{x}]$ through

$\mathcal{U}[\bm{s},\bm{x}]=\mathcal{U}_{0}[\bm{s}]+\Delta\mathcal{U}[\bm{s},\bm{x}]\,.$
(27)

The interaction potential makes it apparent that the distribution of $\bm{s}$
corresponding to the potential $\mathcal{U}[\bm{s},\bm{x}]$ is biased by
$\bm{x}$ with respect to the distribution corresponding to the reference
potential $\mathcal{U}_{0}[\bm{s}]$. By inserting the expressions for
$\mathcal{U}_{0}[\bm{s}]$ and $\mathcal{U}[\bm{s},\bm{x}]$ into Eq. 27 we see
that

$\displaystyle\Delta\mathcal{U}[\bm{s},\bm{x}]$
$\displaystyle=-\ln\mathcal{P}[\bm{x}|\bm{s}]$ (28)
$\displaystyle=-\ln\mathrm{P}(x_{0}|s_{0})-\int^{T}_{0}\mathrm{d}t\
\mathcal{L}_{t}[\bm{s},\bm{x}]$

where $\mathcal{L}_{t}[\bm{s},\bm{x}]$ is given by Eq. 15. This expression
illustrates that the interaction of the output trajectory $\bm{x}$ with the
ensemble of input trajectories is characterized by the trajectory likelihood
$\mathcal{P}[\bm{x}|\bm{s}]$. Because we can compute the trajectory likelihood
from the master equation, we can compute the interaction potential.

In this section we have introduced notation (summarized in Table 1) to show
that computing a marginalization integral is equivalent to the computation of
a free-energy difference. This picture allows us to distinguish two input
trajectory ensembles, the _non-interacting_ ensemble distributed according to
$\exp(-\mathcal{U}_{0}[\bm{s}])=\mathcal{P}[\bm{s}]$, and the _interacting_
ensemble with input distribution proportional to
$\exp(-\mathcal{U}[\bm{s},\bm{x}])\propto\mathcal{P}[\bm{s}|\bm{x}]$. For
example, the brute force estimate of $\mathcal{P}[\bm{x}]$ used in DPWS can be
written as

$\displaystyle\mathcal{P}[\bm{x}]=\frac{\mathcal{Z}[\bm{x}]}{\mathcal{Z}_{0}[\bm{x}]}$
$\displaystyle=\langle e^{-\Delta\mathcal{U}[\bm{s},\bm{x}]}\rangle_{0}$ (29)

where the notation $\langle\cdots\rangle_{0}$ refers to an average with
respect to the non-interacting ensemble. By inserting the expressions for
$\mathcal{U}_{0}$ and $\Delta\mathcal{U}$, it is easy to verify that this
estimate is equivalent to Eq. 8. As explained in Section II.2, to compute Eq.
29 using Monte Carlo, it is only necessary to sample from the non-interacting
system $\mathcal{U}_{0}[\bm{s}]$ and to compute the Boltzmann weight
$\Delta\mathcal{U}[\bm{s},\bm{x}]$ (i.e. to sample from $\mathcal{P}[\bm{s}]$
and to compute the log-likelihood $\ln\mathcal{P}[\bm{x}|\bm{s}]$). This is
indeed the DPWS scheme. However, by noting the correspondence between signal
trajectories and polymers and that Eq. 23 has the same form as the expression
for the (excess) chemical potential of a polymer, which is the free-energy
difference between the polymer of interest and the ideal chain [43, 60], more
efficient schemes can be developed, as we show next.

###  III.2 RR-PWS

Figure 3: Illustration of one step of the bootstrap particle filter in RR-PWS.
We start with a set of trajectories $\bm{s}^{k}_{[0,i-1]}$ with time span
$[\tau_{0},\tau_{i-1}]$ (left panel). In the next step we propagate these
trajectories forward in time to $\tau_{i}$, according to $\mathcal{P}[\bm{s}]$
(central panel). Then we resample the trajectories according to the Boltzmann
weights of their most recent segments, effectively eliminating or duplicating
individual segments. An example outcome of the resampling step is shown in the
right panel where the bottom trajectory was duplicated and one of the top
trajectories was eliminated. These steps are repeated for each segment, until
a set of input trajectories of the desired length is generated. The
intermediate resampling steps bias the trajectory distribution from
$\mathcal{P}[\bm{s}]$ towards $\mathcal{P}[\bm{s}|\bm{x}]$.

In Rosenbluth-Rosenbluth PWS we compute the free-energy difference
$\Delta\mathcal{F}$ between the ideal system $\mathcal{U}_{0}$ and
$\mathcal{U}$ in a _single_ simulation just like in the brute force method.
However, instead of generating $\bm{s}$ trajectories in an uncorrelated
fashion according to $\exp(-\mathcal{U}_{0}[\bm{s}])=\mathcal{P}[\bm{s}]$, we
bias our sampling distribution towards
$\exp(-\mathcal{U}[\bm{s},\bm{x}])\propto\mathcal{P}[\bm{s}|\bm{x}]$ to reduce
the sampling problems found in DPWS.

The classical scheme for biasing the sampling distribution in polymer physics
is due to Rosenbluth and Rosenbluth [61] in their study of self-avoiding
chains. A substantial improvement of the Rosenbluth algorithm was achieved by
Grassberger, by generating polymers using pruning and enrichment steps,
thereby eliminating configurations that do not significantly contribute to the
average. This scheme is known as the pruned-enriched Rosenbluth method, or
PERM [44]. While PERM is much more powerful than the standard Rosenbluth
algorithm, its main drawback is that it requires careful tuning of the pruning
and enrichment schedule to achieve optimal convergence. Therefore we have
opted to use a technique that is similar in spirit to PERM but requires less
tuning, the bootstrap particle filter [62]. We will describe how to use PWS
with a particle filter below. That said, we want to stress that the particle
filter can easily be replaced by PERM or other related methods [63]. Also
schemes inspired by variants of Forward Flux Sampling [64, 65] could be
developed.

In the methods discussed above, a polymer is grown monomer by monomer. In a
continuous-time Markov process this translates to trajectories being grown
segment by segment. To define the segments, we introduce a time discretization
$0<\tau_{1}<\tau_{2}<\cdots<\tau_{n-1}<T$. Thus, each trajectory $\bm{s}$ of
duration $T$ consists of $n$ segments where we denote the segment between
$\tau_{i}$ and $\tau_{j}$ by $\bm{s}_{[i,j]}$ (we define $\tau_{0}=0$ and
$\tau_{n}=T$). The particle filter uses the following procedure to grow an
ensemble of trajectories segment by segment:

  1. 1.

Generate $M$ starting points $s^{1}_{0},\ldots,s^{M}_{0}$ according to the
initial condition of the input signal $\mathrm{P}(s_{0})$.

  2. 2.

Iterate for $i=1,\ldots,n$:

    1. (a)

Starting with an ensemble of $M$ partial trajectories of duration $\tau_{i-1}$
(if $i=1$ an ensemble of starting points) which we label
$\bm{s}^{k}_{[0,i-1]}$ for $k=1,\ldots,M$:

$\left(\bm{s}^{1}_{[0,i-1]},\ldots,\bm{s}^{M}_{[0,i-1]}\right)\,,$ (30)

propagate each trajectory (or each starting point) forward in time from
$\tau_{i-1}$ to $\tau_{i}$. Propagation is performed according to the natural
dynamics of $\bm{s}$, i.e. generating a new segment $\bm{s}^{k}_{[i-1,i]}$
with probability

$p^{\text{gen}}_{i}(k)=\mathcal{P}\left[\bm{s}^{k}_{[i-1,i]}|\bm{s}^{k}_{[0,i-1]}\right]=e^{-\mathcal{U}_{0}\left[\bm{s}^{k}_{[i-1,i]}\right]}$
(31)

for $k=1,\ldots,M$.

    2. (b)

Compute the Boltzmann weight

$U^{k}_{i}=\Delta\mathcal{U}[\bm{s}^{k}_{[i-1,i]},\bm{x}_{[i-1,i]}]$ (32)

of each new segment. This Boltzmann weight of a segment from $\tau_{i-1}$ to
$\tau_{i}$ can be expressed as

$U^{k}_{i}=-\delta_{1i}\ln\mathrm{P}(x_{0}|s_{0})-\int^{\tau_{i}}_{\tau_{i-1}}\mathrm{d}t\
\mathcal{L}_{t}[\bm{s}^{k}_{[i-1,i]},\bm{x}_{[i-1,i]}]\,,$ (33)

see Eq. 28, and is therefore straightforward to compute from the master
equation.

    3. (c)

Sample $M$ times from the distribution

$p^{\text{select}}_{i}(k)=\frac{e^{-U^{k}_{i}}}{w_{i}}$ (34)

where the Rosenbluth weight $w_{i}$ is defined as

$w_{i}=\sum^{M}_{k=1}e^{-U^{k}_{i}}\,.$ (35)

This sampling procedure yields $M$ randomly drawn indices
$\ell^{1}_{i},\ldots,\ell^{M}_{i}$. Each $\ell^{k}_{i}$ is an index that lies
in the range from $1,\ldots,M$ and that points to one of the trajectories that
have been generated up to $\tau_{i}$. To continue the sampling procedure, we
relabel the indices such that the resampled set of trajectories is defined by
$\tilde{\bm{s}}^{k}_{[0,i]}\leftarrow\bm{s}^{\ell^{k}_{i}}_{[0,i]}$ for
$k=1,\ldots,M$. The list
$\left(\tilde{\bm{s}}^{1}_{[0,i]},\ldots,\tilde{\bm{s}}^{M}_{[0,i]}\right)$ is
subsequently used as the input for the next iteration of the algorithm.

The normalized Rosenbluth factor of the final ensemble is then given by

$\mathcal{W}=\prod^{n}_{i=1}\frac{w_{i}}{M}\,.$ (36)

As shown in Appendix A, we can derive an _unbiased_ estimate for the desired
ratio $\mathcal{Z}[\bm{x}]/\mathcal{Z}_{0}[\bm{x}]=\mathcal{P}[\bm{x}]$ based
on the Rosenbluth factor:

$\hat{\mathcal{P}}[\bm{x}]=\mathrm{P}(x_{0})\ \mathcal{W}$ (37)

with $\mathrm{P}(x_{0})$ being the probability of the initial output $x_{0}$.
The particle filter can therefore be integrated into the DPWS algorithm to
compute the marginal density $\mathcal{P}[\bm{x}]$, substituting the brute-
force estimate given in Eq. 8. We call the resulting algorithm to compute the
mutual information _RR-PWS_.

We now provide an intuitive explanation for the scheme presented above. First
note that steps 1 and 2(a) of the procedure above involve just propagating $M$
trajectories in parallel, according to
$\mathcal{P}[\bm{s}]=\exp(-\mathcal{U}_{0}[\bm{s}])$. The interesting steps
are 2(b-c) where we eliminate or duplicate some of the trajectories according
to the Boltzmann weights of the most recent segment. Note, that in general the
list of indices $(\ell^{1}_{i},\ldots,\ell^{M}_{i})$ that are sampled in step
2(c) will contain duplicates ($\ell^{k}_{i}=\ell^{k^{\prime}}_{i}$ for $k\neq
k^{\prime}$), thus cloning the corresponding trajectory. Concomitantly, the
indices $\ell^{1}_{i},\ldots,\ell^{M}_{i}$ may not include every original
index $1,\ldots,M$, therefore eliminating some trajectories. Since indices of
trajectories with high Boltzmann weight are more likely to be sampled from Eq.
34, this ensures that we are only spending computational effort on propagating
those trajectories whose Boltzmann weight is not too small. Hence, at its
heart the particle filter is an algorithm for producing samples that tend to
be distributed according to
$\exp(-\mathcal{U}_{0}[\bm{s}])\exp(-\Delta\mathcal{U}[\bm{s},\bm{x}])=\exp(-\mathcal{U}[\bm{s},\bm{x}])$,
i.e. according to the Boltzmann distribution of the _interacting_ ensemble,
see also Appendix A. For illustration of the algorithm, one iteration of the
particle filter is presented schematically in Fig. 3.

For the efficiency of the particle filter it is important to carefully choose
the number of segments $n$. When the segments are very short ($n$ large), the
accumulated weights (Eq. 33) tend to differ very little between the newly
generated segments $\bm{s}^{k}_{[i-1,i]}$. Hence, the pruning and enrichment
of the segments is dominated by noise. In contrast, when the segments are very
long, the distribution of Boltzmann weights $U^{k}_{i}$ becomes very wide.
Then only few segments contribute substantially to the corresponding
Rosenbluth weight $w_{i}$. Hence, to carefully choose $n$, we need a measure
that quantifies the variance of the trajectory weights. To this ens, we follow
Martino  _et al._ [66] and introduce an effective sample size (ESS)

$M^{\text{(eff)}}_{i}=\frac{w_{i}^{2}}{\sum^{M}_{k=1}\left(e^{-U^{k}_{i}}\right)^{2}},$
(38)

which lies in the range $1\leq M^{\text{(eff)}}_{i}\leq M$;
$M^{\text{(eff)}}_{i}=1$ if one trajectory has a much higher weight than all
the others and $M^{\text{(eff)}}_{i}=M$ if all trajectories have the same
weight. As a rule of thumb, we resample only when the $M^{\text{(eff)}}_{i}$
drops below $M/2$. Additionally, as recommended in Ref. [67], we use the
_systematic sampling_ algorithm to randomly draw the indices in step 2(c)
which helps to reduce the variance; we find, however, the improvement over
simple sampling is very minor. Using these techniques, the only parameter that
needs to be chosen by hand for the particle filter is the ensemble size $M$.

###  III.3 TI-PWS

Our third scheme, thermodynamic integration PWS (TI-PWS), is based on the
analogy of marginalization integrals with free-energy computations. As before,
we view the problem of computing the marginal probability
$\mathcal{P}[\bm{x}]$ as equivalent to that of computing the free-energy
difference between ensembles defined by the potentials
$\mathcal{U}_{0}[\bm{s},\bm{x}]$ and $\mathcal{U}[\bm{s},\bm{x}]$,
respectively. For TI-PWS, we define a potential
$\mathcal{U}_{\theta}[\bm{s},\bm{x}]$ with a continuous parameter
$\theta\in[0,1]$ that allows us to transform the ensemble from
$\mathcal{U}_{0}$ to $\mathcal{U}=\mathcal{U}_{1}$. The corresponding
partition function is

$\mathcal{Z}_{\theta}[\bm{x}]=\int\mathcal{D}[\bm{s}]\
e^{-\mathcal{U}_{\theta}[\bm{s},\bm{x}]}\,.$ (39)

For instance, for $0\leq\theta\leq 1$, we can define our potential as

$\mathcal{U}_{\theta}[\bm{s},\bm{x}]=\mathcal{U}_{0}[\bm{s},\bm{x}]+\theta\,\Delta\mathcal{U}[\bm{s},\bm{x}]\,,$
(40)

such that
$e^{-\mathcal{U}_{\theta}[\bm{s},\bm{x}]}=\mathcal{P}[\bm{s}]\mathcal{P}[\bm{x}|\bm{s}]^{\theta}$.
Note that this is the simplest choice for a continuous transformation between
$\mathcal{U}_{0}$ and $\mathcal{U}_{1}$, but by no means the only one. For
reasons of computational efficiency, it can be beneficial to choose a
different path between $\mathcal{U}_{0}$ and $\mathcal{U}_{1}$, depending on
the specific system [47]. Here we will not consider other paths however, and
derive the thermodynamic integration estimate for the potential given in Eq.
40.

To derive the thermodynamic integration estimate for the free-energy
difference, we first compute the derivative of
$\ln\mathcal{Z}_{\theta}[\bm{x}]$ with respect to $\theta$:

$\displaystyle\frac{\partial}{\partial\theta}\ln\mathcal{Z}_{\theta}[\bm{x}]$
$\displaystyle=\frac{1}{\mathcal{Z}_{\theta}[\bm{x}]}\frac{\partial}{\partial\theta}\int\mathcal{D}[\bm{s}]\
e^{-\mathcal{U}_{\theta}[\bm{s},\bm{x}]}$ (41)
$\displaystyle=-\left\langle\frac{\partial\mathcal{U}_{\theta}[\bm{s},\bm{x}]}{\partial\theta}\right\rangle_{\theta}$
$\displaystyle=-\left\langle\Delta\mathcal{U}[\bm{s},\bm{x}]\right\rangle_{\theta}\,.$

Thus, the derivative of $\ln\mathcal{Z}_{\theta}[\bm{x}]$ is an average of the
Boltzmann weight with respect to $\mathcal{P}_{\theta}[\bm{s}|\bm{x}]$ which
is the ensemble distribution of $\bm{s}$ given by

$\mathcal{P}_{\theta}[\bm{s}|\bm{x}]=\frac{1}{\mathcal{Z}_{\theta}[\bm{x}]}e^{-\mathcal{U}_{\theta}[\bm{s},\bm{x}]}\,.$
(42)

Integrating Eq. 41 with respect to $\theta$ leads to the formula for the free-
energy difference

$\Delta\mathcal{F}[\bm{x}]=-\int^{1}_{0}\mathrm{d}\theta\
\left\langle\Delta\mathcal{U}[\bm{s},\bm{x}]\right\rangle_{\theta}$ (43)

which is the fundamental identity underlying thermodynamic integration.

To compute the free-energy difference using Eq. 43, we evaluate the
$\theta$-integral numerically using Gaussian quadrature, while the inner
average $\left\langle\Delta\mathcal{U}[\bm{s},\bm{x}]\right\rangle_{\theta}$
is computed using MCMC simulations. To perform MCMC simulations in trajectory
space we use ideas from transition path sampling (TPS) [49]. As discussed in
Appendix B, the efficiency of MCMC samplers strongly depends on the proposal
moves that are employed. While better proposal moves could be conceived, we
only use the _forward shooting_ and _backward shooting_ moves of TPS [49] to
obtain the results in Section V. These moves regrow either the end, or the
beginning of a trajectory, respectively. A proposal is accepted according to
the Metropolis criterion [68].

##  IV Integrating Out Internal Components

So far the output trajectory $\bm{x}$ has been considered to correspond to the
trajectory of the system in the _full_ state space $\Omega$. Concomitantly,
the method presented is a scheme for computing the mutual information between
the input signal $\bm{s}$ and the trajectory $\bm{x}$, comprising the time
evolution of all the $n$ components in the system, $X^{1},X^{2},\ldots,X^{n}$.
Each component $X^{i}$ itself has a corresponding trajectory $\bm{x}^{i}$,
such that the full trajectory can be represented as a vector
$\bm{x}=(\bm{x}^{1},\ldots,\bm{x}^{n})$. It is indeed also the conditional
probability
$\mathcal{P}[\bm{x}|\bm{s}]=\mathcal{P}[\bm{x}^{1},\ldots,\bm{x}^{n}|\bm{s}]$
and the marginal probability
$\mathcal{P}[\bm{x}]=\mathcal{P}[\bm{x}^{1},\ldots,\bm{x}^{n}]$ of this vector
in the full state space that can be directly computed from the master
equation. In fact, it is this vector, which captures the states of all the
components in the system, that carries the most information on the input
signal $\bm{s}$, and thus has the largest mutual information. However,
typically the downstream system cannot read out the states of all the
components $X^{1},\ldots,X^{n}$. Often, the downstream system reads out only a
few components or often even just one component, the “output component”
$X^{r}$. The other components then mainly serve to transmit the information
from the input $\bm{s}$ to this readout $X^{r}$. From the perspective of the
downstream system, the other components are hidden. The natural quantity to
measure the precision of information processing is then the mutual information
$\mathrm{I}(\mathcal{S};\mathcal{X}^{r})$ between the input $\bm{s}$ and the
output component’s trajectory $\bm{x}^{r}$, not
$\mathrm{I}(\mathcal{S};\mathcal{X})$. The question then becomes how to
compute $\mathcal{P}[\bm{x}^{r}]$ and $\mathcal{P}[\bm{x}^{r}|\bm{s}]$, from
which $\mathrm{I}(\mathcal{S};\mathcal{X}^{r})$ can be obtained. Here, we
present a scheme to achieve this.

As an example, consider a chemical reaction network with species
$X^{1},\ldots,X^{n}$. Without loss of generality, we will assume that the
$n$-th component is the output component, $X^{r}=X^{n}$. The other species
$X^{1},\ldots,X^{n-1}$ are thus not part of the output, but only relay
information from the input signal $\bm{s}$ to the output signal $\bm{x}^{n}$.
To determine the mutual information $\mathrm{I}(\mathcal{S},\mathcal{X})$ we
need $\mathcal{P}[\bm{x}^{n}|\bm{s}]$, where $\bm{x}^{n}$ is the trajectory of
only the readout component $X^{n}$. However, from the master equation we can
only obtain an expression for the full conditional probability
$\mathcal{P}[\bm{x}^{1},\ldots,\bm{x}^{n}|\bm{s}]$ of all components. To
compute the value of $\mathcal{P}[\bm{x}^{n}|\bm{s}]$, we must perform the
marginalization integral

$\mathcal{P}[\bm{x}^{n}|\bm{s}]=\int\mathcal{D}[\bm{x}^{1}]\cdots\int\mathcal{D}[\bm{x}^{n-1}]\;\mathcal{P}[\bm{x}^{1},\ldots,\bm{x}^{n}|\bm{s}]\,.$
(44)

We can compute this integral using a Monte Carlo scheme as described below and
use the resulting estimate for $\mathcal{P}[\bm{x}^{n}|\bm{s}]$ to compute the
mutual information using our technique presented in Section II.2.

The marginalization of Eq. 44 entails integrating out degrees of freedom from
a known joint probability distribution. In Eq. 8 we solved the analogous
problem of obtaining the marginal probability $\mathcal{P}[\bm{x}]$ by
integrating out the input trajectories through the integral
$\mathcal{P}[\bm{x}]=\int\mathrm{d}\bm{s}\
\mathcal{P}[\bm{s},\bm{x}]=\int\mathrm{d}\bm{s}\
\mathcal{P}[\bm{s}]\mathcal{P}[\bm{x}|\bm{s}]$. As described in Section II.2,
the integral from Eq. 8 can be computed via a Monte Carlo estimate by sampling
many input trajectories from $\mathcal{P}[\bm{s}]$ and taking the average of
the corresponding conditional probabilities $\mathcal{P}[\bm{x}|\bm{s}_{i}]$.
We will show that in the case where there is no feedback from the readout
component back to the other components, a completely analogous Monte Carlo
estimate can be derived for Eq. 44. We describe this below. Additionally, in
the absence of feedback, the techniques presented in Section III above can be
employed to develop computationally more efficient schemes.

More specifically, we can evaluate Eq. 44 via a direct Monte Carlo estimate
under the condition that the stochastic dynamics of the other components
$X^{1},\ldots,X^{n-1}$ are not influenced by $X^{n}$ (i.e., no feedback from
the readout). Using the identity

$\mathcal{P}[\bm{x}^{1},\ldots,\bm{x}^{n}|\bm{s}]=\mathcal{P}[\bm{x}^{1},\ldots,\bm{x}^{n-1}|\bm{s}]\
\mathcal{P}[\bm{x}^{n}|\bm{x}^{1}_{i},\ldots,\bm{x}^{n-1}_{i},\bm{s}]$ (45)

to rewrite the integrand in Eq. 44, we are able to represent the conditional
probability $\mathcal{P}[\bm{x}^{n}|\bm{s}]$ as an average over the readout
component’s trajectory probability

$\mathcal{P}[\bm{x}^{n}|\bm{s}]=\left\langle\mathcal{P}[\bm{x}^{n}|\bm{x}^{1}_{i},\ldots,\bm{x}^{n-1}_{i},\bm{s}]\right\rangle_{\mathcal{P}[\bm{x}^{1},\ldots,\bm{x}^{n-1}|\bm{s}]}\,.$
(46)

Thus, assuming that we can evaluate the conditional probability of the readout
given all the other components,
$\mathcal{P}[\bm{x}^{n}|\bm{x}^{1}_{i},\ldots,\bm{x}^{n-1}_{i},\bm{s}]$, we
arrive at the estimate

$\mathcal{P}[\bm{x}^{n}|\bm{s}]\approx\frac{1}{M}\sum^{M}_{i=1}\mathcal{P}[\bm{x}^{n}|\bm{x}^{1}_{i},\ldots,\bm{x}^{n-1}_{i},\bm{s}]$
(47)

where the samples $\bm{x}^{1}_{i},\ldots,\bm{x}^{n-1}_{i}$ for $i=1,\ldots,M$
are drawn from $\mathcal{P}[\bm{x}^{1},\ldots,\bm{x}^{n-1}|\bm{s}]$. Notice
that the derivation of this Monte Carlo estimate is fully analogous to the
estimate in Eq. 8, but instead of integrating out the input trajectory
$\bm{s}$ we integrate out the component trajectories
$\bm{x}^{1},\ldots,\bm{x}^{n-1}$.

To obtain
$\mathcal{P}[\bm{x}^{n}|\bm{x}^{1}_{i},\ldots,\bm{x}^{n-1}_{i},\bm{s}]$ in
Eqs. 46 and 47, we note that, in absence of feedback, we can describe the
stochastic dynamics of the readout component $X^{n}$ as a jump process with
time-dependent transition rates whose time-dependence arises from the
trajectories of the other components $\bm{x}^{1},\ldots,\bm{x}^{n-1}$ and the
input input $\bm{s}$. In effect, this is a driven jump process for $X^{n}$,
driven by all upstream components $X^{1},\ldots,X^{n-1}$ and the input signal.
Specifically, denoting $\bm{u}=(\bm{x}^{1},\ldots,\bm{x}^{n-1},\bm{s})$ as the
joint trajectory representing the history of all upstream components as well
as the input signal, we can, as explained in Section II.3, write the time
dependent transition rate matrix $Q_{t}(\cdot|\bm{u})$ for the stochastic
dynamics of $X^{n}$ and use Eq. 14 to compute
$\mathcal{P}[\bm{x}^{n}|\bm{u}]=\mathcal{P}[\bm{x}^{n}|\bm{x}^{1}_{i},\ldots,\bm{x}^{n-1}_{i},\bm{s}]$.
Using Eq. 47, this then allows us to compute $\mathcal{P}[\bm{x}^{n}|\bm{s}]$.

Finally, to compute the mutual information
$\mathrm{I}(\mathcal{S};\mathcal{X}^{n})$, e.g. using the estimate in Eq. 10,
we additionally need to evaluate the marginal output probability
$\mathcal{P}[\bm{x}^{n}]$. This requires us to perform one additional
integration over the space of input trajectories $\bm{s}$:

$\displaystyle\mathcal{P}[\bm{x}^{n}]$ $\displaystyle=\int\mathcal{D}[\bm{s}]\
\mathcal{P}[\bm{s}]\mathcal{P}[\bm{x}^{n}|\bm{s}]$ (48)
$\displaystyle=\left\langle\mathcal{P}[\bm{x}^{n}|\bm{s}]\right\rangle_{\mathcal{P}[\bm{s}]}\,.$

The corresponding Monte Carlo estimate is

$\displaystyle\mathcal{P}[\bm{x}^{n}]$
$\displaystyle\approx\frac{1}{N}\sum^{N}_{i=1}\mathcal{P}[\bm{x}^{n}|\bm{s}_{i}]$
(49)
$\displaystyle\approx\frac{1}{N}\sum^{N}_{i=1}\frac{1}{M}\sum^{M}_{j=1}\mathcal{P}[\bm{x}^{n}|\bm{x}^{1}_{ij},\ldots,\bm{x}^{n-1}_{ij},\bm{s}_{i}]$

where the input trajectories $\bm{s}_{i}$ follow $\mathcal{P}[\bm{s}]$ and the
intermediate components $(\bm{x}^{1}_{ij},\ldots,\bm{x}^{n-1}_{ij})$, for
$i=1,\ldots,N$ and $j=1,\ldots,M$, follow
$\mathcal{P}[\bm{x}^{1},\ldots,\bm{x}^{n-1}|\bm{s}_{i}]$.

In summary, the scheme to obtain $\mathcal{P}[\bm{x}^{n}|\bm{u}]$ in the
presence of hidden intermediate components is analogous to that used for
computing $\mathcal{P}[\bm{x}]$ from $\mathcal{P}[\bm{x}|\bm{s}]$. In both
cases, one needs to marginalize a distribution function by integrating out
components. Indeed, the schemes presented here and in Section II.2 are bona
fide schemes to compute the mutual information between the input $\bm{s}$ and
either the trajectory of the output component $\bm{x}^{n}$ or the full output
$\bm{x}$. However, when the trajectories are sufficiently long or the
stochastic dynamics are sufficiently complex, then the free-energy schemes of
Section III may be necessary to enhance the efficiency of computing the
marginalized distribution, $\mathcal{P}[\bm{x}]$ or
$\mathcal{P}[\bm{x}^{n}|\bm{s}]$.

##  V Results

To demonstrate the power of our framework and illustrate how the techniques of
the previous sections can be used in practice, we apply PWS to two instructive
chemical reaction networks. We first consider a linearly coupled birth-death
process. This system has already been studied previously using a Gaussian
model [6], and by Duso and Zechner [30] using an approximate technique, and we
compare our results to those of these studies. This simple birth-death system
serves to illustrate the main ideas of our approach and also highlights that
linear systems can be distinctly non-Gaussian. The second example has been
chosen to demonstrate the practical applicability of our technique. We use RR-
PWS to compute the mutual information rate in the bacterial chemotaxis system,
which is a prime example of a complex information processing system consisting
of many reactions. Then we compare the computed rate against recent
experiments.

The code used to produce the PWS estimates was written in the Julia
programming language [69] and has been made freely available [70,
*pws_github]. For performing stochastic simulations we use the
DifferentialEquations.jl package [72] and biochemical reaction models are set
up with help from the ModelingToolkit.jl package [73].

###  V.1 Coupled Birth-Death Processes

As a first example system we consider a simple birth-death process
$\emptyset\rightleftharpoons\mathrm{X}$ of species $\mathrm{X}$ which is
created at rate $\rho(t)$ and decays with constant rate $\mu$ per copy of
$\mathrm{X}$. This system receives information from an input signal that
modulates the birth rate $\rho(t)$. For simplicity, we assume it is given by

$\rho(t)=\rho_{0}s(t)$ (50)

where $\rho_{0}$ is a constant and $s(t)$ is the input copy number at time
$t$. This is a simple model for gene expression, where the rate of production
of a protein $\mathrm{X}$ is controlled by a transcription factor
$\mathrm{S}$, and $\mathrm{X}$ itself has a characteristic decay rate. The
input trajectories $s(t)$ themselves are generated via a separate birth-death
process $\emptyset\rightleftharpoons\mathrm{S}$ with production rate $\kappa$
and decay rate $\lambda$.

We compute the trajectory mutual information for this system as a function of
the trajectory duration $T$ of the input and output trajectories. For
$T\rightarrow\infty$, the trajectory mutual information is expected to
increase linearly with $T$, since, on average, every additional output segment
contains the same additional amount of information on the input trajectory.
Because we are interested in the mutual information in steady state, the
initial states $(s_{0},x_{0})$ were drawn from the stationary distribution
$\mathrm{P}(s_{0},x_{0})$. This distribution was obtained using a Gaussian
approximation. This does not influence the asymptotic rate of increase of the
mutual information, but leads to a nonzero mutual information already for
$T=0$.

Figure 4: Comparing different schemes to compute the mutual information as a
function of trajectory duration for a simple coupled birth-death process with
rates $\kappa=50,\lambda=1,\rho_{0}=10,\mu=10$ and steady-state initial
condition. The top panels show example trajectories of input and output as
well as the mean (solid line) and standard deviation (shaded region). Below,
the mutual information is shown as a function of trajectory duration. The
inset shows an enlarged version of the dotted rectangle near the origin. For
short trajectories all PWS estimates agree. Yet, for longer trajectories, DPWS
and TI-PWS require a much larger number of input trajectories $M$ for
computing $\mathcal{P}[\bm{x}]$ than RR-PWS to converge. Results for the three
PWS variants are compared with the Duso and Zechner [30] estimate, and with
the linear noise approximation from Ref. [6]. We find excellent agreement
between the Duso scheme and RR-PWS. The Gaussian linear noise approximation
systematically underestimates the mutual information. All PWS estimates, as
well as the Duso approximation were computed using $N=10^{4}$ samples from
$\mathcal{P}[\bm{s},\bm{x}]$.

Figure 4 shows the mutual information as a function of the trajectory duration
$T$. We compare the three PWS variants and two approximate schemes. One is
that of Duso and Zechner [30]. To apply it, we used the code publicly provided
by the authors
222https://github.com/zechnerlab/PathMI/tree/302f03e51ad195adc6be39fa9618886c76590cc4,
and to avoid making modifications to this code, we chose a fixed initial
condition $(s_{0}=x_{0}=50)$ which causes the mutual information to be zero
for $T=0$. The figure also shows the analytical result of a Gaussian model
[6], obtained using the linear-noise approximation (see Appendix F).

We find that the efficiency of the respective PWS variants depends on the
duration of the input-output trajectories. For short trajectories all PWS
variants yield very similar estimates for the mutual information. However, for
longer trajectories the estimates of DPWS and, to a smaller degree, TI-PWS
diverge, because of poor sampling of the trajectory space in the estimate of
$\mathcal{P}[\bm{x}]$. For longer trajectories, the estimate becomes
increasingly dominated by rare trajectories, which make an exceptionally large
contribution to the average of $\mathcal{P}[\bm{x}]$. Missing these rare
trajectories with a high weight tends to increases the marginal entropy
$\mathrm{H}(\mathcal{X})$ (see Eq. 9), and thereby the mutual information;
indeed, the estimates of DPWS and TI-PWS are higher than that of RR-PWS. For
brute-force DPWS, the error decreases as we increase the number $M$ of input
trajectories per output trajectory used to estimate $\mathcal{P}[\bm{x}]$.
Similarly, for TI-PWS the error decreases as we use more MCMC samples for the
marginalization scheme. For the RR-PWS, however, already for $M=128$ the
estimate has converged; we verified that a further increase of $M$ does not
change the results.

We also find excellent agreement between the RR-PWS estimate and the
approximate result of Duso and Zechner [30]. Only very small deviations are
visible in Fig. 4. These deviations are mostly caused by the different choice
for the initial conditions. In RR-PWS, the initial conditions are drawn from
the stationary distribution, while in the Duso scheme they are fixed, such
that the mutual information computed with RR-PWS is finite while that computed
with the Duso scheme is zero. Yet, as the trajectory duration $T$ increases,
the Duso estimate slowly “catches up” with the RR-PWS result.

Fig. 4 also shows that although the Gaussian model matches the PWS result for
$T=0$, it systematically underestimates the mutual information for
trajectories of finite duration $T>0$. Interestingly, this is not a
consequence of small copy-number fluctuations: increasing the average copy
number does not significantly improve the Gaussian estimate. We leave a
detailed analysis of this observation for future work.

Figure 5: Comparing estimation bias for the different PWS variants in relation
to their CPU time requirements. Each dot represents a single mutual
information estimate with $N=10^{4}$ samples for trajectories of duration
$T=5$. Almost all the CPU time of a PWS estimate is spent on the computation
of the marginal probability $\mathcal{P}[\bm{x}]$. The bias of the marginal
probability estimate can be reduced by using a larger number $M$ of sampled
input trajectories to compute the marginalization integral, which also
increases the required CPU time. The RR-PWS estimate converges much faster
than the estimate of DPWS and TI-PWS. For DPWS and TI-PWS, the dots represents
estimates ranging from $M=2^{5}$ to $M=2^{14}$, for RR-PWS ranging from
$M=2^{3}$ to $M=2^{10}$. As the baseline of zero bias we use the converged
result from the RR-PWS estimates.

The different approaches for computing the marginal probability
$\mathcal{P}[\bm{x}]$ lead to different computational efficiencies of the
respective PWS schemes. In Fig. 5, as a benchmark, we show the magnitude of
the error of the different PWS estimates in relation to the required CPU time.
Indeed, as expected, the computation of the marginal probability poses
problems for long trajectories when using the brute force DPWS scheme. More
interestingly, while TI-PWS improves the estimate of the mutual information,
the improvement is not dramatic. Unlike the brute-force scheme, thermodynamic
integration does make it possible to generate input trajectories $\bm{s}$ that
are correlated with the output trajectories $\bm{x}$, but it still
overestimates the mutual information for long trajectories unless a very large
number of MCMC samples are used.

The RR-PWS implementation evidently outperforms the other estimates for this
system. The regular resampling steps ensure that we mostly sample input
trajectories $\bm{s}$ with non-vanishing likelihood
$\mathcal{P}[\bm{x}|\bm{s}]$, thereby avoiding the sampling problem from DPWS.
Moreover, sequential Monte Carlo techniques such as RR-PWS and FFS [64] have a
considerable advantage over MCMC techniques in trajectory sampling. With MCMC
path sampling, we frequently make small changes to an existing trajectory such
that the system moves slowly in path space, leading to poor statistics. In
contrast, in RR-PWS we generate new trajectories from scratch, segment by
segment, and these explore the trajectory space much faster.

The coupled birth-death process represents a simple yet non-trivial system
capable of information transmission. In the next section we apply PWS to a
more complex and realistic biochemical signaling network.

###  V.2 Bacterial Chemotaxis

The chemotaxis system of the bacterium Escherichia coli is a complex
information processing system. It is responsible for detecting nutrient
gradients in the cell’s environment and using that information to guide the
bacterium’s movement. Briefly, _E. coli_ navigates through its environment by
performing a biased random walk, successively alternating between so-called
runs, during which it swims with a nearly constant speed, and tumbles, during
which it randomly chooses a new direction [75]. The rates of switching between
these two states are controlled by the chemotaxis sensing system (Fig. 6a).
This system consists of receptors on the cell surface that detect the ligand,
and a downstream signaling network that processes this information by taking
the time-derivative of the signal, the ligand concentration. This derivative
is taken via two antagonistic reactions, which occur on two distinct
timescales. Attractant binding rapidly deactivates the receptor, while slow
methylation counteracts this reaction by reactivating the receptor, leading to
near perfect adaptation [76, 77, 78, 79]. Lastly, active receptors
phosphorylate the downstream messenger protein CheY, which controls the
tumbling propensity by binding the flagellar motors that propel the bacterium.

Figure 6:  The information transmission rate in the bacterial chemotaxis
system. a Cartoon of the chemotaxis network of _E. coli_. Receptors form
clusters with an associated CheA kinase. A cluster can either be active or
inactive, depending on the number of bound ligands (green dots) and methylated
sites (orange dots). Active CheA can phosphorylate CheY; phosphorylated CheY
controls the rotation direction of the flagellar motors and thereby the
movement of the bacterium. b In a shallow gradient $c(x)$, the bacterium
diffuses nearly freely in the $x$-direction. The variance of the position
increases with time, the hallmark of a random walk. The input signal is the
concentration $c(t)=c(x(t))$ as experienced by the bacterium at time $t$. c
This movement gives rise to the input statistics of the signal. The shaded
regions indicate the 75th and 95th percentiles, example trajectories are
displayed in color. This signal is non-stationary as its variance always keeps
growing. d Mutual information $I(\mathbf{C}_{T},\mathbf{Y}_{T})$ between input
trajectories $c(t)$ and output trajectories $y_{p}(t)$ as a function of
trajectory duration $T$. In each RR-PWS simulation, $N=7200$ Monte Carlo
samples were used ($M=256$ for the particle filter). e The information
transmission rate is defined as $I(\mathbf{C}_{T},\mathbf{Y}_{T})/T$ in the
limit $T\to\infty$.

In our model, the receptors are grouped in clusters (Appendix D). Each
receptor can switch between an active and an inactive conformational state,
but, in the spirit of the Monod-Wyman-Changeux model [80], the energetic cost
of having two different conformations in the same cluster is prohibitively
large. We can then speak of each cluster as either being active or inactive.
Each receptor in a cluster can bind ligand and be (de)methylated, which,
together, control the probability that the cluster is active. In the
simulations, receptor (de)methylation is modeled explicitly, because the
(de)methylation reactions are slow. In contrast, the timescale of receptor-
ligand (un)binding is much faster than the other timescales in the system,
i.e., those of the input dynamics, CheY (de)phosphorylation, and receptor
(de)methylation. The receptor-ligand binding dynamics can therefore be
integrated out without affecting information transmission, in order to avoid
wasting CPU time (Appendix D). In addition, the receptor clusters can
phosphorylate CheY, while phosphorylated CheY is dephosphorylated at a
constant rate. The dynamics of the kinase CheA and the phosphatase CheZ which
drive (de)phosphorylation are not modeled explicitly. Table 2 in Appendix D
gives the parameter values of our chemotaxis model, which are all based on
values reported in the literature. For what follows below, the key parameters
are the number of receptors per cluster, which is taken to be $N=6$ based on
Refs. [81, 82], while the number of clusters is $N_{\rm c}=N_{\rm r}/N=400$,
where $10^{3}<N_{\rm r}<10^{4}$ is an estimate for the total number of
receptors based on Ref. [83]. The translation of our model into a master
equation is explained in Appendix D.

We first asked whether this model based on the current literature can
reproduce the information transmission rate as recently measured by Mattingly
_et al._ [50]. In what follows, we call this model the “literature-based”
model.

The information transmission rate depends not only on the biochemical
chemotaxis network, but also on the dynamics of the input signal. It is
therefore important that the dynamics of this signal in our model agree with
those in the experiments of Mattingly  _et al._ [50]. In our model, the input
signal is the time-dependent ligand concentration $c(t)$ that is experienced
by the swimming bacterium. In the experiments, the cells swim in a very
shallow chemical gradient, such that their swimming behavior is, to a good
approximation, identical to that in the absence of a gradient. The cell’s
movement can thus be modeled as a (persistent) random walk. Assuming that the
gradient is oriented in $x$-direction, the autocorrelation of the
$x$-component of the velocity was experimentally found to be well described by
$V(t)\equiv\langle\delta v_{x}(t)\delta
v_{x}(0)\rangle\simeq\sigma^{2}_{v}e^{-\lambda|t|}$, where $\sigma^{2}_{v}$ is
the variance of the fluctuations in the velocity, and $\lambda^{-1}$ is the
correlation time of these fluctuations [50]. We therefore model the cell’s
velocity as $\dot{v}_{x}=-\lambda v_{x}+\xi$, which, with
$\langle\xi(t)\xi(t^{\prime})\rangle=2\sigma^{2}_{v}\lambda\delta(t-t^{\prime})$,
gives rise to the measured correlation function $V(t)$. These cells swim in a
shallow, exponential gradient $c(x)\propto e^{gx}$ with steepness $g$. The
dynamics of the ligand concentration as experienced by the bacteria are then
given by $\dot{c}=gcv_{x}$. The cell’s own swimming dynamics described by
$v_{x}$ thus give rise to the input signal $c(t)$ of the biochemical network.
It constitutes an exponential random walk, as illustrated in Fig. 6b, c (see
also Appendix E).

In our model, the output is the concentration of phosphorylated CheY, while in
the experiments of Mattingly  _et al._ [50] it is the average activity of the
receptor clusters as obtained via FRET measurements. We argue that this
difference does not significantly affect the obtained information rates, and
thus, that it is valid to compare our results to the experiments. In
particular, since the copy number of CheY is much larger than the number of
receptor clusters, the fluctuations in CheY are dominated by the extrinsic
fluctuations coming from the receptor activity noise rather than from the
intrinsic fluctuations associated with CheY (de)phosphorylation. To a good
approximation, the copy number of phosphorylated CheY, ${\rm Yp}(t)$, is thus
a deterministic function of the average receptor activity $a(t)$.
Mathematically, the mutual information $I(X;Y)$ between two stochastic
variables $X$ and $Y$ is the same as the mutual information $I(f(X);g(Y))$ for
deterministic and monotonic functions $f$ and $g$. It follows that the mutual
information between $c(t)$ and ${\rm Yp}(t)$, is nearly the same as that
between $c(t)$ and the receptor activity $a(t)$. It is therefore meaningful to
compare the information transmission rates as predicted by our PWS simulations
to those measured by Mattingly  _et al._ [50].

We use RR-PWS to exactly compute the mutual information for the literature-
based model. Specifically, we measure the mutual information
$I(\mathbf{C},\mathbf{Y_{p}};T)$ between the input trajectory of the ligand
concentration $c(t)$ and the output trajectory of phosphorylated CheY,
$y_{p}(t)$, and where each trajectory is of duration $T$. With RR-PWS it is
possible to compute $I(\mathbf{C},\mathbf{Y_{p}};\tau)$ for all $\tau\leq T$
within a single PWS simulation of duration $T$ by saving intermediate results
after each sampled segment, see Section III.2. The receptor states are hidden
internal states, and we use the technique of Section IV to integrate them out.

Figure 7: Comparison of theoretical models with experimental data for
bacterial chemotaxis system. Panels a and b show the response and noise
kernels, respectively, for the model based on literature parameters (green),
parameters fitted to experiments (blue), and experiments from Mattingly  _et
al._ [50] (orange). In panel c, the information transmission rate is shown for
each model as a function of gradient steepness, with results from the Gaussian
approximation shown alongside exact PWS calculations. The fitted model closely
matches the experiments, while the literature-based model over-estimates
information transmission rate by a factor of $\approx 4$ despite having a
lower response amplitude (panel a). This is because the literature-based model
has a large number of independents receptor clusters $N_{c}$, resulting in
much lower noise in the output (panel b). In all cases, the Gaussian
approximation matches the exact PWS results, providing support for the
accuracy of the measurements by Mattingly  _et al._ [50].

Figure 6e shows the PWS estimate of the information transmission rate for
cells swimming in gradients of different steepnesses $g$. The information
transmission rate is obtained from the PWS estimate of the trajectory mutual
information $I(\mathbf{C},\mathbf{Y_{p}};T)$, different trajectory durations
$T$. As seen in Fig. 6d, for short trajectories the mutual information
increases non-linearly with trajectory duration $T$, but in the long-duration
limit the slope becomes constant. This asymptotic rate of increase of the
mutual information with $T$ is the desired information transmission rate
$R(\mathbf{C},\mathbf{Y_{p}})$. The precise definition is given by

$R(\mathbf{C},\mathbf{Y_{p}})=\lim_{T\to\infty}\frac{I(\mathbf{C},\mathbf{Y_{p}};T)}{T}\,.$
(51)

We then compared our results for the information transmission rate of the
literature-based model to those of Mattingly  _et al._ [50]. Figure 7c shows
that the model predictions differ from the experiments by a factor of $\approx
4$. Despite this discrepancy, we believe that the agreement between experiment
and theory is, in fact, remarkable, because these predictions were made ab
initio: the model was developed based on the existing literature and we did
not fit our model to the data of Mattingly  _et al._

Yet, the question about the origin of the discrepancy remains. The difference
between their measurements and our predictions could be attributed either to
the inaccuracy of our model or to the approximation that Mattingly  _et al._
had to employ to compute the information transmission rate from experimental
data. Concerning the latter hypothesis, due to the curse of dimensionality and
experimental constraints, Mattingly  _et al._ could not directly obtain the
information transmission rate from measured time traces of the input and
output of the system. Instead, they measured three different kernels that
describe the system in the linear regime. Specifically, they obtained the
response $K(t)$ of the kinase activity to a step-change in input signal, the
autocorrelation function of the input signal $V(t)$, and the autocorrelation
$N(t)$ of the kinase activity in a constant background concentration. Then
they used a Gaussian model to compute the information transmission rate from
these measured functions $K(t)$, $V(t)$, and $N(t)$ [6, 50] (see also Appendix
F). This Gaussian model is based on a linear noise assumption and cannot
perfectly capture the true non-linear dynamics of the biochemical network.
This could be the cause for the observed discrepancies in the information
rate. We have indeed already seen in Section V.1 that there can be substantial
differences between exact computations and the Gaussian approximation for the
trajectory mutual information.

To uncover the reason for the discrepancy we first tested whether our
literature-based model reproduces the experimentally measured kernels. If the
kernels do not match, then, clearly, the discrepancy in the information rate
may be caused by the difference between our model and the experimental system,
as opposed to the inaccuracy of the Gaussian framework. Our input correlation
function, $V(t)$, is, by construction, the same as that of Mattingly  _et al._
[50]. However, we find that the response kernel $K(t)$ and the autocorrelation
function of the noise $N(t)$ of our system are different. Figure 7a, b shows
that our model reproduces the timescales of $N(t)$ and $K(t)$ as measured
experimentally. This is perhaps not surprising, because the decay of both
$N(t)$ and $K(t)$ is set by the (de)methylation rate, which has been well-
characterized experimentally. Yet, the figure also shows that our model
significantly underestimates the amplitudes of both $N(t)$ and $K(t)$.

This raises the question of whether other parameter values would allow our
model to better reproduce the measured kernels $K(t)$ and $N(t)$, and,
secondly, whether this would resolve the discrepancy in information rate
between our simulations and the experiments.

The amplitude $\sigma^{2}_{N}$ of the output noise correlation function $N(t)$
is bounded by the number of receptor clusters $N_{\rm c}$. In particular, the
variance of the receptor activity is $\sigma^{2}_{N}=\sigma^{2}_{a}/N_{\rm
c}\leq 1/4N_{\rm c}$, where $\sigma^{2}_{a}\leq 1/4$ is the variance of the
activity of a single receptor cluster. Comparing this bound to the measured
receptor noise strength $\sigma^{2}_{N}$ reveals that $N_{\rm c}$ needs to be
much smaller than our original model assumes: the number of clusters needs to
be as small as $N_{\rm c}\lesssim 10$. Indeed, Fig. 7b shows that with $N_{\rm
c}=9$, our model quantitatively fits the correlation function $N(t)$ of the
receptor activity in a constant background concentration, as measured
experimentally [50].

The amplitude of $K(t)$, i.e. the gain, depends on the ratio $K_{\rm D}^{\rm
A}/K_{\rm D}^{\rm I}$ of the dissociation constants of the receptor for ligand
binding in its active or inactive state, respectively, as well as on the
number of receptors per cluster, $N$. Both dissociation constants have been
well characterized experimentally [84, 81], but the number of receptors per
cluster has only been inferred indirectly from experiments [82, 81]. The
higher gain as measured experimentally by Mattingly  _et al._ [50] indicates
that $N$ is larger than assumed in our model: with $N=15$ our model can
quantitatively fit $K(t)$ (Fig. 7a).

We thus find that by reducing the number of clusters from $N_{\rm c}=400$ to
$N_{\rm c}=9$ while simultaneously increasing their size from $N=6$ to $N=15$,
our model is able to quantitatively fit both $N(t)$ and $K(t)$ [50] (Fig. 7).
This suggests that the number of independent receptor clusters is smaller than
hitherto believed, while their size is larger.

Finally, how accurately can our revised model reproduce the measured
information rate, and how accurate is the Gaussian framework for the
experimental system in the regime studied by Mattingly  _et al._ [50]? In the
revised model, called the “fitted model”, with $N_{\rm c}=9$ and $N=15$, all
key quantities for computing the information transmission rate within the
Gaussian framework, $V(t)$, $N(t)$ and $K(t)$, are nearly identical to the
experiments of Mattingly  _et al._ [50], see Fig. 7. Within the Gaussian
framework (see Appendix F), the information transmission rate in our model is
thus expected to be very similar to the experimentally measured one, and Fig.
7c shows that this is indeed the case. To quantify the accuracy of the
Gaussian framework, we then recomputed the information transmission rate for
the revised model, using exact PWS (see Appendix H). We found that the result
matches the Gaussian prediction very well. For these shallow and static
chemical gradients, the Gaussian model is thus highly accurate. Our analysis
validates _a posteriori_ the Gaussian framework adopted by Mattingly  _et al._
[50].

##  VI Discussion

In this manuscript, we have developed a general, practical, and flexible
method that makes it possible to compute the mutual information between
trajectories exactly. PWS is a Monte Carlo scheme based on the exact
computation of trajectory probabilities. We showed how to compute exact
trajectory probabilities from the master equation and thus how to use PWS for
any system described by a master equation. Since the master equation is
employed in many fields and in particular provides an exact stochastic model
for well-mixed chemical reaction dynamics, PWS is very broadly applicable.

The application of PWS to the bacterial chemotaxis system shows how crucial it
is to have a simulation technique that is exact. Without the latter it would
be impossible to determine whether the difference between our predictions and
the Mattingly data [50] is due to the inaccuracy of the model, the inaccuracy
of the numerical technique to simulate the model, or the approximations used
by Mattingly and coworkers in analyzing the data. In contrast, because PWS is
exact, we knew the difference between theory and experiment is either due to
the inaccuracy of the model or the approximations used to analyze the data. By
then employing the same Gaussian framework to analyze the behavior of the
model and the experimental system, we were able to establish that the
difference is due to the inaccuracy of our original model.

Our analysis indicates that the size of the receptor clusters in the E. coli
chemotaxis system, $N\approx 15$, is larger than that based on previous
estimates, $N\sim 6$ [85, 81, 82]. The early estimates of the cluster size
were based on bulk dose-response measurements with a relatively slow ligand
exchange, yielding $N\approx 6$ [85, 81]. More recent dose-response
measurements, at the single cell level and with faster ligand exchange, yield
an average that is higher, $\langle N\rangle\approx 8$, and with a broad
distribution around it, arising from cell-to-cell variability [82]. Our
estimate, $N\approx 15$, based on fitting the response kernel $K(t)$ to that
measured by Mattingly  _et al._ [50], therefore appears reasonable. At the
same time, the number of clusters, obtained by fitting the noise correlation
function $N(t)$ to the data of Mattingly  _et al._ [50] is surprisingly low,
$N_{c}\sim 10$, given the total number of receptors, $N_{r}\sim
10^{3}\text{--}10^{4}$ [83]. Interestingly, recent experiments indicate that
the receptor array is poised near the critical point [86], where receptor
switching becomes correlated over large distances. This effectively partitions
receptors into a few large domains, which may explain our fitted values for
$N$ and $N_{c}$.

It has been suggested that information processing systems are positioned close
to a critical point to maximize information transmission [87, 21], although it
has been argued that the sensing error of the E. coli chemotaxis system is
minimized for independent receptors [88]. Mattingly  _et al._ have
demonstrated that the chemotactic drift speed in shallow exponential gradients
is limited by the information transmission rate [50], but whether the system
has been optimized for information transmission, and how the latter affects
chemotactic performance in other spatio-temporal concentration profiles,
remain interesting questions for future work.

While we have focused on the computation of the mutual information between
trajectories of systems governed by a master equation, this concept can be
extended to other types of stochastic processes. The crux of PWS is the exact
evaluation of the path likelihood $\mathcal{P}[\bm{x}|\bm{s}]$ from the master
equation. Therefore, in order to use PWS with a different stochastic process,
we require the ability to compute the path likelihood
$\mathcal{P}[\bm{x}|\bm{s}]$ for its sample paths. Remarkably, for _stochastic
diffusion processes_ , which are described by a Langevin equation in a
continuous state space, the trajectory probability can be computed by
replacing the kernel $\mathcal{L}_{t}(\bm{s},\bm{x})$ in Eq. 15 with the
Onsager-Machlup function [89]. In particular, while the path probability is
not well-defined for continuous sample paths, Adib [90] shows that the
Onsager-Machlup function yields a correct expression for the path probability
of time-discretized diffusion trajectories. Therefore, PWS can be extended to
handle systems described by a Langevin equation. In such a PWS scheme, the
Gillespie simulations are replaced by standard numerical integration
techniques for Langevin equations (see e.g. Ref. [58]), and the trajectory
likelihood is evaluated on-the-fly using the Onsager-Machlup function.

PWS cannot be used to directly obtain the mutual information between
trajectories from experimental data, in contrast to model-free (yet
approximate) methods such as K-nearest-neighbors estimators [23, 24],
decoding-based information estimates [26], or schemes that compute the mutual
information from the data within the Gaussian framework [50]. PWS requires a
(generative) model based on a master equation or Langevin equation. However,
an increasingly popular approach is to estimate the mutual information from
experimental data using hidden Markov models (HMM) [91, 92, 93]. While the
construction of these HMM models is beyond the scope of this manuscript, PWS
makes it possible to compute the mutual information between time-varying
inputs and outputs exactly for these models.

We have applied PWS to compute the mutual information (rate) in steady state,
but PWS can be used equally well to study systems out of steady state. For
such systems a (non-)equilibrium initial condition $\mathrm{P}(s_{0},x_{0})$
must be specified in addition to a well-defined non-stationary probability
distribution of input trajectories $\mathcal{P}[\bm{s}]$. These distributions
are defined by the (experimental) setup and lead to a well-defined output
distribution $\mathcal{P}[\bm{x}]$ when the system is coupled to the input.
Thus, a steady state is no prerequisite for the application of PWS to study
the trajectory mutual information.

Throughout the manuscript, we have considered systems in which the output does
not feed back onto the input. In systems with feedback, the current output
influences future input, which means that we cannot straightforwardly generate
input trajectories according to $\mathcal{P}[\bm{s}]$. Moreover,
$\mathcal{P}[\bm{x}|\bm{s}]$, the central quantity of all three PWS methods,
cannot be obtained straightforwardly. Nonetheless, PWS can be extended to
systems with feedback as shown in Appendix C. As in free-energy calculations,
the trick is to define a reference system for which the marginal probability
distribution—and thus the “free energy”—is known and then compute the “free-
energy difference” between that reference system and the system of interest.

Aside from DPWS, we developed two additional variants of PWS, capitalizing on
the connection between information theory and statistical physics.
Specifically, the computation of the mutual information requires the
evaluation of the marginal probability of individual output trajectories
$\mathcal{P}[\bm{x}]$. This corresponds to the computation of a partition
function in statistical physics. RR-PWS and TI-PWS are based on techniques
from polymer and rare-event sampling to make the computation of the marginal
trajectory probability more efficient.

The different PWS variants share some characteristics yet also differ in
others. DPWS and RR-PWS are static Monte Carlo schemes in which the
trajectories are generated independently from the previous ones. These methods
are similar to static polymer sampling schemes like PERM [44] and rare-event
methods like DFFS or BG-FFS [64]. In contrast, TI-PWS is a dynamic Monte Carlo
scheme, where a new trajectory is generated from the previous trajectory. In
this regard, this method is similar to the CBMC scheme for polymer simulations
[94] and the TPS [49], TIS [95], and RB-FFS [64] schemes to harvest transition
paths. The benefit of static schemes is that the newly generated trajectories
are uncorrelated from the previous ones, which means that they are less likely
to get stuck in certain regions of path space. Concomitantly, they tend to
diffuse faster through the configuration space. Indeed, TI-PWS suffers from a
problem that is also often encountered in TPS or TIS, which is that the middle
sections of the trajectories move only slowly in their perpendicular
direction. Tricks that have been applied to TPS and TIS to solve this problem,
such as parallel tempering, could also be of use here [96].

Another distinction is that RR-PWS generates all the trajectories in the
ensemble simultaneously yet segment by segment, like DFFS, while DPWS and TI-
PWS generate only one full trajectory at the time, similar to RB-FFS, BG-FFS,
and also TPS and TIS. Consequently, RR-PWS, like DFFS, faces the risk of
_genetic drift_ , which means that, after sufficiently many resampling steps,
most paths of the ensemble will originate from the same initial seed. Thus,
when continuing to sample new segments, the old segments that are far in the
past become essentially fixed, which makes it possible to miss important paths
in the RR-PWS sampling procedure. As in DFFS, the risk of genetic drift in RR-
PWS can be mitigated by increasing the initial number of path segments.
Although we did not employ this trick here, we found that RR-PWS was by far
the most powerful scheme of the three variants studied.

Nonetheless, we expect that DPWS and TI-PWS become more efficient in systems
that respond to the input signal with a significant delay $\tau$. In these
cases, the weight of a particular output trajectory depends on the degree to
which the dynamics of the output trajectory correlates with the dynamics of
the intput trajectory a time $\tau$ earlier. Because in RR-PWS a new segment
of an output trajectory is generated based on the corresponding segment of the
input trajectory that spans the same time-interval, it may therefore miss
these correlations between the dynamics of the output and that of the input a
time $\tau$ earlier. In contrast, DPWS and TI-PWS generate full trajectories
one at the time, and are therefore more likely to capture these correlations.

Overall, PWS is a general framework for computing the mutual information
between trajectories. We presented three variants of PWS for systems described
by a master equation. Apart from these, we expect that other variants could be
developed to improve efficiency for particular applications. Because of its
flexibility and simplicity, we envision that PWS will become an important and
reliable tool for studying information transmission in dynamic stochastic
systems.

###### Acknowledgements.

We thank Bela Mulder, Tom Shimizu, Fotios Avgidis, Peter Bolhuis, and Daan
Frenkel for useful discussions and a careful reading of the manuscript, and we
thank Age Tjalma for support with obtaining the Gaussian approximation of the
chemotaxis system. This work is part of the Dutch Research Council (NWO) and
was performed at the research institute AMOLF. This project has received
funding from the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation program (grant agreement No. 885065), and
was financially supported by NWO through the “Building a Synthetic Cell
(BaSyC)” Gravitation grant (024.003.019).

##  Appendix A Justification for the Particle Filter from Section III.2

Here we justify the marginal probability estimate shown in Eq. 37, i.e. we
show that the bootstrap particle filter used in Section III.2 provides a
consistent estimate for the ratio of partition functions
$\mathcal{P}[\bm{x}]=\mathcal{Z}[\bm{x}]/\mathcal{Z}_{0}[\bm{x}]$. The result
that this estimate is also _unbiased_ is more difficult to establish; a proof
is given by Del Moral [97].

To make the derivations below easy to follow, we structure our justification
of the particle filter into three steps. We first give a brief description of
how a resampling procedure can generally be used to generate samples according
to a target distribution when only samples from a different distribution are
available. Secondly, we use these insights to explain how the resampling
procedure used in the particle filter generates trajectories that are
distributed approximately according to $\mathcal{P}[\bm{s}|\bm{x}]$, even
though we only generate trajectories according to $\mathcal{P}[\bm{s}]$.
Finally, we use this result to show that the particle filter provides a
consistent estimate for $\mathcal{P}[\bm{x}]$.

###  A.1 Sampling and Resampling

Sampling and then resampling is a strategy to use samples
$\bm{s}^{1},\ldots,\bm{s}^{M}$ from a given prior distribution $f[\bm{s}]$ to
generate _approximate_ samples from a different distribution of interest, with
density proportional to the product $h[\bm{s}]=f[\bm{s}]g[\bm{s}]$. In
general, $h[\bm{s}]$ is not normalized, and we denote the corresponding
normalized probability density by
$\hat{h}[\bm{s}]=h[\bm{s}]/\int\mathcal{D}[\bm{s}]h[\bm{s}]$. To generate
samples from $\hat{h}[\bm{s}]$, we assign each of the existing samples from
$f[\bm{s}]$ a normalized weight

$W^{k}=\frac{g[\bm{s}^{k}]}{\sum^{M}_{j=1}g[\bm{s}^{j}]}\,.$ (52)

Then, by sampling from the discrete set $\\{\bm{s}^{1},\ldots,\bm{s}^{M}\\}$
according to the weights $W^{1},\ldots,W^{M}$, we pick samples that are
approximately distributed according to $\hat{h}[\bm{s}]$. Indeed, for
$M\rightarrow\infty$ the distribution of the resulting samples approaches the
density $\hat{h}[\bm{s}]$ [98]. We use resampling at each iteration of the
algorithm of Section III.2 to regularly prune those trajectories with low
overall contribution to the marginalization integral.

###  A.2 Distribution of Trajectories in the Particle Filter

In the bootstrap particle filter, at each iteration, we start with a set of
trajectories $\bm{s}^{1}_{[0,i-1]},\ldots,\bm{s}^{M}_{[0,i-1]}$ which we
assume are approximately distributed according to
$\mathcal{P}[\bm{s}_{[0,i-1]}|\bm{x}_{[0,i-1]}]$. In each iteration of the
particle filter, the goal is to produce a set of trajectories approximately
distributed according to $\mathcal{P}[\bm{s}_{[0,i]}|\bm{x}_{[0,i]}]$.
Clearly, by iterating such a procedure, we can generate a set of trajectories
distributed approximately according to
$\mathcal{P}[\bm{s}_{[0,n]}|\bm{x}_{[0,n]}]$ for any $n>1$. Note that we
always carefully use the phrase _approximately distributed_ because, as
explained above, for finite $M$, a resampling procedure cannot generate
_exact_ samples from a probability distribution (yet the estimate for
$\mathcal{P}[\bm{x}]$ remains unbiased regardless of how good these
approximations are). We now take a closer look at one iteration of the
particle filter.

We start with the set of trajectories with a time span
$[\tau_{0},\tau_{i-1}]$, denoted by
$\left\\{\bm{s}^{1}_{[0,i-1]},\ldots,\bm{s}^{M}_{[0,i-1]}\right\\}$. These
trajectories are then propagated forward to time $\tau_{i}$, by adding a new
segment $\bm{s}^{k}_{[i-1,i]}$ to the trajectory $\bm{s}^{k}_{[0,i-1]}$ for
$k=1,\ldots,M$. Each new segment is generated from the distribution
$\mathcal{P}[\bm{s}^{k}_{[i-1,i]}|\bm{s}^{k}_{[0,i-1]}]$ such that the
propagation step results in a set of trajectories
$\\{\bm{s}^{1}_{[0,i]},\ldots,\bm{s}^{M}_{[0,i]}\\}$, distributed according to
$f[\bm{s}_{[0,i]}]=\mathcal{P}[\bm{s}_{[0,i]}|\bm{x}_{[0,i-1]}]$.

Next, we resample from the set of trajectories, with the goal of producing a
set of trajectories distributed according to the target density
$\hat{h}[\bm{s}]=\mathcal{P}[\bm{s}_{[0,i]}|\bm{x}_{[0,i]}]$. Thus, we have to
find the appropriate weighting function $g[\bm{s}_{[0,i]}]$ in order to
approximately produce samples according to the target distribution. By
choosing
$g[\bm{s}_{[0,i]}]=\exp\left\\{-\Delta\mathcal{U}[\bm{s}_{[i-1,i]},\bm{x}_{[i-1,i]}]\right\\}=\mathcal{P}[\bm{x}_{[i-1,i]}|\bm{x}_{[0,i-1]},\bm{s}_{[0,i]}]$,
we generate normalized weights

$W^{k}_{i}=\frac{\mathcal{P}[\bm{x}_{[i-1,i]}|\bm{x}_{[0,i-1]},\bm{s}^{k}_{[0,i]}]}{\sum^{M}_{j=1}\mathcal{P}[\bm{x}_{[i-1,i]}|\bm{x}_{[0,i-1]},\bm{s}^{j}_{[0,i]}]}\,,$
(53)

cf. Eq. 52. Note that this is the same choice of weighting function as in
Section III.2, Eq. 34. By comparison with the notation used there, we see that
the Boltzmann factors $U^{k}_{i}$ and Rosenbluth weights $w_{i}$ were defined
such that we can express the normalized weight equivalently as

$W^{k}_{i}=\frac{e^{-U^{k}_{i}}}{w_{i}}\,.$ (54)

Why is this choice of weighting function the correct one? First, observe that
resampling with the normalized weights of Eq. 53 produces samples
approximately distributed according to

$\displaystyle h[\bm{s}_{[0,i]}]$
$\displaystyle=f[\bm{s}_{[0,i]}]g[\bm{s}_{[0,i]}]$ (55)
$\displaystyle=\mathcal{P}[\bm{s}_{[0,i]}|\bm{x}_{[0,i-1]}]\
\mathcal{P}[\bm{x}_{[i-1,i]}|\bm{x}_{[0,i-1]},\bm{s}_{[0,i]}]\,.$

What remains to be shown is that this density $h[\bm{s}_{[0,i]}]$, when
normalized, becomes the desired target distribution
$\mathcal{P}[\bm{s}_{[0,i]}|\bm{x}_{[0,i]}]$.

To do so, we need to rewrite the expression for
$g[\bm{s}_{[0,i]}]=\mathcal{P}[\bm{x}_{[i-1,i]}|\bm{x}_{[0,i-1]},\bm{s}_{[0,i]}]$
using Bayes’ theorem

$g[\bm{s}_{[0,i]}]=\frac{\mathcal{P}[\bm{s}_{[0,i]}|\bm{x}_{[0,i-1]},\bm{x}_{[i-1,i]}]\
\mathcal{P}[\bm{x}_{[i-1,i]}|\bm{x}_{[0,i-1]}]}{\mathcal{P}[\bm{s}_{[0,i]}|\bm{x}_{[0,i-1]}]}\,.$
(56)

Notice that the first term of the numerator can be written as
$\mathcal{P}[\bm{s}_{[0,i]}|\bm{x}_{[0,i]}]$. After inserting this result into
Eq. 55, we obtain

$h[\bm{s}_{[0,i]}]=\mathcal{P}[\bm{s}_{[0,i]}|\bm{x}_{[0,i]}]\
\mathcal{P}[\bm{x}_{[i-1,i]}|\bm{x}_{[0,i-1]}]\,.$ (57)

The second term in this product is a constant, since $\bm{x}$ is fixed. The
first term is a normalized probability density for $\bm{s}_{[0,i]}$. Therefore
we find that the normalized density corresponding to $h[\bm{s}_{[0,i]}]$ is

$\hat{h}[\bm{s}_{[0,i]}]=\mathcal{P}[\bm{s}_{[0,i]}|\bm{x}_{[0,i]}]\,.$ (58)

Consequently, this is the distribution that is approximated by the set of
trajectories at the end of the $i$-th iteration of the particle filter, which
is what we wanted to show. At its heart, the particle filter is therefore an
algorithm to produce samples that are approximately distributed according to
$\mathcal{P}[\bm{s}|\bm{x}]$.

###  A.3 Marginal Probability Estimate

We now use these insights to derive an estimate of the marginal density
$\mathcal{P}[\bm{x}]$. We start by noting that the marginal density of the
$i$-th output segment, $\mathcal{P}[\bm{x}_{[i-1,i]}|\bm{x}_{[0,i-1]}]$, is
given by

$\displaystyle\mathcal{P}[\bm{x}_{[i-1,i]}|\bm{x}_{[0,i-1]}]$ (59)
$\displaystyle=\int\mathcal{D}[\bm{s}_{[0,i]}]\
\mathcal{P}[\bm{x}_{[i-1,i]},\bm{s}_{[0,i]}|\bm{x}_{[0,i-1]}]$
$\displaystyle=\int\mathcal{D}[\bm{s}_{[0,i]}]\
\mathcal{P}[\bm{s}_{[0,i]}|\bm{x}_{[0,i-1]}]\ g[\bm{s}_{[0,i]}]\,.$

The third line follows from the definition of
$g[\bm{s}_{[0,i]}]=\mathcal{P}[\bm{x}_{[i-1,i]}|\bm{x}_{[0,i-1]},\bm{s}_{[0,i]}]$.
Hence, we find that the probability
$\mathcal{P}[\bm{x}_{[i-1,i]}|\bm{x}_{[0,i-1]}]$ can be expressed as the
average

$\mathcal{P}[\bm{x}_{[i-1,i]}|\bm{x}_{[0,i-1]}]=\left\langle
g[\bm{s}_{[0,i]}]\right\rangle_{\mathcal{P}[\bm{s}_{[0,i]}|\bm{x}_{[0,i-1]}]}\,.$
(60)

In principle, this average can be computed using a Monte Carlo scheme, using
trajectories generated from $\mathcal{P}[\bm{s}_{[0,i]}|\bm{x}_{[0,i-1]}]$.
Notice that at each iteration of the particle filter, we _do_ dispose of a set
of trajectories $\bm{s}^{1}_{[0,i]},\ldots,\bm{s}^{M}_{[0,i]}$ which are
approximately distributed according to
$\mathcal{P}[\bm{s}_{[0,i]}|\bm{x}_{[0,i-1]}]$ above. Therefore, we can
compute the average Eq. 60 directly from the trajectories that are present for
each iteration of the particle filter. With the notation from Section III.2,
using $g[\bm{s}^{k}_{[0,i]}]=\exp(-U^{k}_{i})$, we thus obtain the estimate

$\mathcal{P}[\bm{x}_{[i-1,i]}|\bm{x}_{[0,i-1]}]\approx\frac{1}{M}\sum^{M}_{k=1}e^{-U^{k}_{i}}=\frac{w_{i}}{M}\,.$
(61)

The probability of the entire output trajectory $\mathcal{P}[\bm{x}]$ is given
by the product

$\mathcal{P}[\bm{x}]=\mathrm{P}(x_{0})\mathcal{P}[\bm{x}_{[0,1]}|x_{0}]\cdots\mathcal{P}[\bm{x}_{[n-1,n]}|\bm{x}_{[0,n-1]}]$
(62)

where $\mathrm{P}(x_{0})$ is the probability of the initial output state
$x_{0}$ which is assumed to be known. In conclusion, we arrive at the
following estimate for the marginal output probability

$\hat{\mathcal{P}}[\bm{x}]=\mathrm{P}(x_{0})\prod^{n}_{i=1}\frac{w_{i}}{M}$
(63)

which is precisely Eq. 37.

##  Appendix B MCMC Sampling in Trajectory Space

Thermodynamic Integration PWS in Section III.3 relies on the computation of
averages with respect to
$\mathcal{P}_{\theta}[\bm{s}|\bm{x}]\propto\exp(-\mathcal{U}_{\theta}[\bm{s},\bm{x}])$.
Sampling from these distributions using the SSA (Gillespie) algorithm is not
possible. Instead, in this section, we show different ways of how to implement
a Markov Chain Monte Carlo (MCMC) sampler in trajectory space to generate
correctly distributed trajectories.

We can build an MCMC sampler in trajectory space using the Metropolis-Hastings
algorithm. To create a Markov Chain in trajectory space, we need to find a
suitable proposal kernel, that generates a new trajectory $\bm{s}^{\prime}$
from a given trajectory $\bm{s}$ with probability
$T(\bm{s}\rightarrow\bm{s}^{\prime})$. We accept the proposal using the
Metropolis criterion with probability

$A(\bm{s}^{\prime},\bm{s})=\min\left(1,e^{\mathcal{U}_{\theta}[\bm{s},\bm{x}]-\mathcal{U}_{\theta}[\bm{s}^{\prime},\bm{x}]}\frac{T(\bm{s}^{\prime}\rightarrow\bm{s})}{T(\bm{s}\rightarrow\bm{s}^{\prime})}\right)$
(64)

to create a chain of trajectories with stationary distribution given by
$\mathcal{P}_{\theta}[\bm{s}|\bm{x}]=e^{-\mathcal{U}_{\theta}[\bm{s},\bm{x}]}/\mathcal{Z}_{\theta}[\bm{x}]$
for $0\leq\theta\leq 1$. To ensure efficient convergence of the resulting
Markov chain to its stationary distribution, the proposal kernel must balance
two conflicting requirements. To efficiently explore the state space per unit
amount of CPU time, the proposed trajectory $\bm{s}^{\prime}$ must be
sufficiently different from the original trajectory $\bm{s}$, while at the
same time it should not be so radically different that the acceptance
probability is drastically reduced. Thus, the design of the proposal kernel is
crucial for an efficient MCMC sampler, and we will discuss various strategies
to create trial trajectories. Since different types of trial moves can easily
be combined in a Metropolis-Hastings algorithm, the most efficient samplers
often incorporate multiple complementary proposal strategies to improve the
exploration speed of the trajectory space.

The simplest (and naïve) proposal kernel is to generate an entirely new
trajectory $\bm{s}^{\prime}$ independent of $\bm{s}$, by sampling directly
from $\mathcal{P}[\bm{s}]$ using the SSA. Hence, the transition kernel is
given by $T(\bm{s}\rightarrow\bm{s}^{\prime})=\mathcal{P}[\bm{s}^{\prime}]$
and a proposal $\bm{s}\rightarrow\bm{s}^{\prime}$ is accepted with probability

$\displaystyle A(\bm{s}^{\prime},\bm{s})$
$\displaystyle=\min\left(1,e^{\mathcal{U}_{\theta}[\bm{s},\bm{x}]-\mathcal{U}_{\theta}[\bm{s}^{\prime},\bm{x}]}\frac{\mathcal{P}[\bm{s}]}{\mathcal{P}[\bm{s}^{\prime}]}\right)$
(65)
$\displaystyle=\min\left(1,\frac{\mathcal{P}[\bm{x}|\bm{s}^{\prime}]^{\theta}}{\mathcal{P}[\bm{x}|\bm{s}]^{\theta}}\right)$

where the second line follows by inserting the definition of
$\mathcal{U}_{\theta}[\bm{s},\bm{x}]$ given in Eq. 40. Although this simple
scheme to completely regenerate an entire trajectory and accepting/rejecting
according to $A(\bm{s}^{\prime},\bm{s})$ creates correctly distributed
trajectories, it should not be used in simulations to compute
$\mathcal{P}[\bm{x}]$. Indeed, we get a better estimate of
$\mathcal{P}[\bm{x}]$ by just using the same number of independent sample
trajectories from $\mathcal{P}[\bm{s}]$ and using the brute-force scheme in
Eq. 8 without taking the detour of thermodynamic integration to estimate the
normalization constant.

Instead, an idea from transition path sampling is to only regenerate a part of
the old trajectory as part of the proposal kernel [99]. By not regenerating
the entire trajectory, the new trajectory $\bm{s}^{\prime}$ is going to be
correlated with the original trajectory $\bm{s}$, and correlation in general
improves the acceptance rate. The simplest way to generate trial trajectories
using a partial update is a move termed _forward shooting_ in which a time
point $\tau$ along the existing trajectory $\bm{s}$ is randomly selected, and
a new trajectory segment is regrown from this point to the end, resulting in
the proposal $\bm{s}^{\prime}$. Since the new segment is generated according
to the input statistics given by $\mathcal{P}[\bm{s}_{[T-\tau,T]}]$, the
acceptance probability for the proposed trajectory is given by Eq. 65. If the
input dynamics given by $\mathcal{P}[\bm{s}]$ are time-reversible, we can also
perform a _backward shooting_ move. Here, the beginning of $\bm{s}$ is
replaced by a new segment that is generated backwards in time. Assuming that
the initial condition is the input’s steady state distribution, the
corresponding acceptance probability of the backward shooting move is again
given by Eq. 65. Using these two moves we create an MCMC sampler where both
ends of the trajectory are flexible, and thus if the trajectory is not too
long, the chain will quickly relax to its stationary distribution. This is
indeed the MCMC sampler used to obtain the TI-PWS results for the coupled
birth-death process in Section V.1.

For long trajectories it can prove to be a problem that the middle section is
too inflexible when the proposal moves only regenerate either the beginning or
the end of a trajectory. Therefore, one could additionally try to incorporate
mid-section regrowth to make sure that also the middle parts of the trajectory
become flexible. To regrow a middle segment with duration $\tau$ of a
trajectory $\bm{s}$, we have to generate a new segment of duration $\tau$
according to the stochastic dynamics given by $\mathcal{P}[\bm{s}]$ but with
the additional condition that we have to connect _both_ endpoints of the new
segment to the existing trajectory. Although the starting point of the segment
can be freely chosen, the challenge is to ensure that the end point of the new
segment satisfies the end-point constraint. Stochastic processes that generate
trajectories under the condition of hitting a specific point after a given
duration $\tau$ are called stochastic bridging processes.

The simplest way to generate trajectories from a bridging process is by
generating a trajectory segment of length $\tau$ from the normal stochastic
process and rejecting the segment if it does not hit the correct end point
[100]. Clearly, this strategy is only feasible for very short segments and
when the state space is discrete, as otherwise almost every generated segment
will be rejected due to not hitting the correct end point. To avoid this
problem, more efficient algorithms have been developed to simulate stochastic
bridges for some types of stochastic processes. For diffusion processes,
bridges can be simulated efficiently by introducing a guiding term into the
corresponding Langevin equation [101]. For jump processes, bridges can be
simulated using particle filters [102], by a weighted stochastic simulation
algorithm (wSSA) [103], or using random time-discretization (uniformization)
[100].

Further techniques to create a trajectory space MCMC samplers have been
developed in the literature. Crooks [104] describes a scheme to create MCMC
moves for trajectories evolving in non-equilibrium dynamics, by making MCMC
moves to change the trajectories’ noise histories. In the Particle Markov
Chain Monte Carlo (PMCMC) algorithm, proposal trajectories are generated using
a particle filter and accepted with an appropriate Metropolis criterion [105].
Another class of efficient samplers for Markov jump processes can be built
using uniformization [106].

##  Appendix C Dealing with Feedback

Although so far, we have assumed the stochastic dynamics of the input to be
independent of the generated output trajectories, in principle all physical
information processing systems exhibit feedback. The physical interaction
needed to measure the input signal necessarily affects the incoming signal,
and indeed, it follows that no information can be extracted from the input
signal without any perturbation of the input dynamics. Often, it is assumed
that the amplitude of such perturbations is comparatively small and thus that
the feedback can safely be ignored. Above, the PWS scheme was derived with
this assumption. In this section, we drop the assumption and will explicitly
consider systems where the produced output perturbs the input, i.e. systems
where the output feeds back onto the input. In the following we will first
discuss the additional problems that arise when computing the mutual
information for a system with feedback, and subsequently we present a modified
version of PWS that can be used to compute the trajectory mutual information
for these systems.

###  C.1 Computing the Mutual Information with Feedback between Input and
Output

All PWS schemes presented above require the computation of the trajectory
likelihood $\mathcal{P}[\bm{x}|\bm{s}]$, a quantity that is not readily
available for systems with feedback. Indeed, as already mentioned in Section
II.3.1, for a given input trajectory $\bm{s}$, the output dynamics are no
longer described by a Markov process in a system with feedback, and therefore
we cannot find an expression for $\mathcal{P}[\bm{x}|\bm{s}]$ based on the
master equation. This implies that for systems with feedback, PWS schemes
cannot be used without modification. While it is generally not possible to
derive an expression for the conditional probability
$\mathcal{P}[\bm{x}|\bm{s}]$ in systems with feedback, we often still can
compute the joint probability density $\mathcal{P}[\bm{s},\bm{x}]$ instead.
Based on this quantity, we will present a modified PWS scheme to compute the
mutual information for systems with feedback.

Specifically, since PWS is a model-based approach to compute the mutual
information, when there is feedback from the output back to the input, we
require a complete model of the combined system. Specifically, such a model
must provide an expression for the joint probability
$\mathcal{P}[\bm{s},\bm{x}]$, describing the input dynamics and the
interaction between input and output, including the feedback.

An estimate of the mutual information that only relies on the computation of
joint probability densities $\mathcal{P}[\bm{s},\bm{x}]$ can be obtained by
expressing the mutual information as

$\mathrm{I}(\mathcal{S},\mathcal{X})=\int\mathcal{D}[\bm{s}]\int\mathcal{D}[\bm{x}]\
\mathcal{P}[\bm{s},\bm{x}]\ln\frac{\mathcal{P}[\bm{s},\bm{x}]}{\mathcal{P}[\bm{s}]\,\mathcal{P}[\bm{x}]}\,.$
(66)

Thus, the PWS scheme with feedback consists of the computation of

$\mathrm{I}(\mathcal{S},\mathcal{X})=\left\langle\ln\frac{\mathcal{P}[\bm{s},\bm{x}]}{\mathcal{P}[\bm{s}]\,\mathcal{P}[\bm{x}]}\right\rangle_{\mathcal{P}[\bm{s},\bm{x}]}$
(67)

which we want to estimate via a Monte Carlo average using samples from
$\mathcal{P}[\bm{s},\bm{x}]$. We see that while we don’t need to evaluate the
likelihood $\mathcal{P}[\bm{x}|\bm{s}]$, we now need to explicitly compute the
joint density $\mathcal{P}[\bm{s},\bm{x}]$, and two marginal densities,
$\mathcal{P}[\bm{s}]$ and $\mathcal{P}[\bm{x}]$, for each Monte Carlo sample
$(\bm{s},\bm{x})\sim\mathcal{P}[\bm{s},\bm{x}]$. While the joint density can
be evaluated directly by assumption, each of the marginalized densities can
only be computed using a nested Monte Carlo estimate.

Specifically, for PWS with feedback, we need to compute _two_ marginalization
integrals per Monte Carlo sample:

$\mathcal{P}[\bm{s}]=\int\mathcal{D}[\bm{x}]\ \mathcal{P}[\bm{s},\bm{x}]\,,$
(68)

and

$\mathcal{P}[\bm{x}]=\int\mathcal{D}[\bm{s}]\ \mathcal{P}[\bm{s},\bm{x}]\,.$
(69)

However, these marginalization integrals cannot be directly computed with the
techniques described so far. Note that while in Section III we discussed in
detail how to compute such marginalization integrals, all methods presented
there themselves require the evaluation of the likelihood
$\mathcal{P}[\bm{x}|\bm{s}]$ and cannot be used directly. Therefore, in the
following subsection, we discuss how to compute marginalization integrals for
systems with feedback.

Additionally, as discussed in Section IV, we may also need to integrate out
internal components of the master equation even when the output feeds back
onto these internal components. The technique discussed below can also be used
in this case as a way to compute the marginalization integral in Eq. 44.

###  C.2 Marginalization Integrals for Systems with Feedback

Computing marginalization integrals in systems with feedback is harder than it