Mutual Information Rate -- Linear Noise Approximation and Exact Computation

TL;DR

This paper evaluates the accuracy of the Gaussian approximation for mutual information rate estimation in nonlinear and discrete systems, using Path Weight Sampling for exact computation and highlighting when more precise methods are needed.

Contribution

It provides a quantitative assessment of the Gaussian approximation's limitations and demonstrates the effectiveness of Path Weight Sampling for exact mutual information rate calculation.

Findings

Gaussian approximation often fails in nearly Gaussian systems

Significant deviations occur in nonlinear diffusive systems

Path Weight Sampling offers exact mutual information estimates

Abstract

Efficient information processing is crucial for both living organisms and engineered systems. The mutual information rate, a core concept of information theory, quantifies the amount of information shared between the trajectories of input and output signals, and enables the quantification of information flow in dynamic systems. A common approach for estimating the mutual information rate is the Gaussian approximation which assumes that the input and output trajectories follow Gaussian statistics. However, this method is limited to linear systems, and its accuracy in nonlinear or discrete systems remains unclear. In this work, we assess the accuracy of the Gaussian approximation for non-Gaussian systems by leveraging Path Weight Sampling (PWS), a recent technique for exactly computing the mutual information rate. In two case studies, we examine the limitations of the Gaussian approximation. First, we focus on discrete linear systems and demonstrate that, even when the system's statistics are nearly Gaussian, the Gaussian approximation fails to accurately estimate the mutual information rate. Second, we explore a continuous diffusive system with a nonlinear transfer function, revealing significant deviations between the Gaussian approximation and the exact mutual information rate as nonlinearity increases. Our results provide a quantitative evaluation of the Gaussian approximation's performance across different stochastic models and highlight when more computationally intensive methods, such as PWS, are necessary.

Equations

$$I(S,X)=\iint P(s,x)\ln \frac{P(s,x)}{P(s)P(x)}dsdx,$$

$$I(S,X)$$

$$=H(S)-H(S|X)$$

$$=H(X)-H(X|S).$$

$$I({\mathbf{S}}_T,{\mathbf{X}}_T)=⟨\ln \frac{P({\mathbf{s}}_T,{\mathbf{x}}_T)}{P({\mathbf{s}}_T)P({\mathbf{x}}_T)}⟩$$

$$R(\mathcal{S},\mathcal{X})=\underset{T\to \infty }{\lim }\frac{1}{T}I({\mathbf{S}}_T,{\mathbf{X}}_T),$$

$$I({\mathbf{S}}_T,{\mathbf{X}}_T)=\frac{1}{2}\ln \frac{|{\mathbf{C}}_{ss}||{\mathbf{C}}_{xx}|}{|\mathbf{Z}|},$$

$$\mathbf{Z}=\left(\begin{array}{cc}{\mathbf{C}}_{ss}& {\mathbf{C}}_{sx}\\ {\mathbf{C}}_{xs}& {\mathbf{C}}_{xx}\end{array}\right)$$

$$R(\mathcal{S},\mathcal{X})=-\frac{1}{4\pi }{\int }_{-\infty }^{\infty }d\omega \ln \left(1-\frac{|S_{sx}(\omega ){|}^2}{S_{ss}(\omega )S_{xx}(\omega )}\right).$$

$${\varphi }_{sx}(\omega )=\frac{|S_{sx}(\omega ){|}^2}{S_{ss}(\omega )S_{xx}(\omega )}$$

$$\frac{\sum^{N}_{i=1}\left[\ln\operatorname{\mathrm{P}}\left(\bm{x}^{i}_{T}\Bigm{|}\bm{s}^{i}_{T}\right)-\ln\operatorname{\mathrm{P}}\left(\bm{x}^{i}_{T}\right)\right]}{N}$$

$$P({\mathbf{x}}_T)=\int d{\mathbf{s}}_{T}P({\mathbf{s}}_T)P({\mathbf{x}}_T\mid {\mathbf{s}}_T)$$

$$\ln P(\mathbf{x}\mid \mathbf{s})=-\sum _{i=1}^{n-1}\frac{1}{2\Delta t}{\left(\frac{\Delta x_i-v_i\Delta t}{\sigma (x_i)}\right)}^2+\text{const}$$

$$\dot{s}(t)$$

$$\dot{x}(t)$$

$$R_{\text{Gaussian}}=\frac{\lambda }{2}\left(\sqrt{1+\frac{\rho }{\lambda }}-1\right).$$

$$R_{\text{discrete}}=\frac{\lambda }{2}\left(\sqrt{1+2\frac{\rho }{\lambda }}-1\right)$$

$$\dot{x}(t)=\rho a(s)-\mu x(t)+{\eta }_x(t)$$

$$a(s)=\{\begin{array}{ll}\frac{s^n}{K^n+s^n}& \text{if }s\ge 0,\\ 0& \text{if }s<0.\end{array}$$

$$g=\frac{\partial\bar{x}(s)}{\partial s}\bigg{|}_{\bar{s}}=\frac{na(\bar{s})[1-a(\bar{s})]\rho}{\mu\bar{s}}\,,$$

$$R_{\text{LNA}}=\frac{\lambda }{2}\left(\sqrt{1+g^2\frac{\bar{s}\mu }{\bar{x}\lambda }}-1\right).$$

$${\hat{\varphi }}_{sx}(\omega )=\frac{|{\hat{S}}_{sx}(\omega ){|}^2}{{\hat{S}}_{ss}(\omega ){\hat{S}}_{xx}(\omega )}$$

$$<=>[\kappa ][\lambda ]S,$$

$$\dot{s}=\kappa -\lambda s(t)+{\eta }_s(t),$$

$$⟨{\eta }_s^2⟩=\kappa +\lambda \bar{s}=2\lambda \bar{s}.$$

$$S_{ss}(\omega )$$

$$->[\rho a(s)]S+X,$$

$$->[\mu ]\varnothing .$$

$$a(s)=\frac{s(t{)}^n}{K^n+s(t{)}^n}.$$

$$\dot{x}=\rho a(s)-\mu x(t)+{\eta }_x(t),$$

$$g$$

$$=n(1-\bar{a})\bar{x}/\bar{s},$$

$$r=n\bar{a}(1-\bar{a})\rho /\bar{s},$$

Full text

I Introduction

For the functioning of both living and engineered systems it is paramount that they collect and process information effectively. Increasingly, it has become clear that beyond instantaneous properties, the dynamic features of an input signal or system output often encode valuable information [1, 2, 3, 4, 5, 6]. Prime examples in biology include bacterial chemotaxis, which responds to temporal changes in concentration [7], the transcription factor NF-$\kappa$B, which encodes information about input signals in its dynamic response [8], and neuronal information processing, where information is encoded in the sequence and timing of spikes [9]. Beyond biology, dynamic input signals are critical for various sensing systems, such as those used in automated factories or self-driving cars.

To understand and evaluate the performance, potential improvements, and limitations of these systems in processing information, we need appropriate metrics that capture their full information processing capability. Information theory, introduced by Shannon [10], provides the most general mathematical framework for such metrics. The mutual information measures how much one random variable reduces uncertainty about another, quantified in bits. It is relatively straightforward to quantify the information shared between scalar properties of the input and output, as has been done in various forms [11, 12, 13, 14, 15, 16, 17, 18, 19]. However, capturing all information in dynamical properties of the input and the output is much more challenging. To do so, one must consider the information encoded in time-varying trajectories of the variables of interest. Yet, due to the high dimensionality of the trajectory space, computing the mutual information between such trajectories is notoriously difficult.

A major advance for biological systems has been the Gaussian approximation of the mutual information rate [1, 2], based on the assumption of input and output trajectories following jointly Gaussian statistics. This assumption makes it possible to compute the mutual information rate directly from the two-point correlation functions of the input and output. It is thus straightforward to apply the Gaussian approximation to experimental data. Moreover, given a mechanistic model of the underlying dynamics, the Gaussian approximation can be used to derive analytical expressions for the information rate [1, 2, 3]. Crucially, however, the assumption of Gaussian statistics restricts the method to linear systems with additive Gaussian noise [20].

Understanding when the Gaussian approximation is accurate is critical because many real-world systems, such as biological and engineered sensory systems, exhibit nonlinear dynamics. This includes features such as bimodality, discrete jumps, or heavy tails, which deviate from Gaussian dynamics. Such non-Gaussian behavior typically results from intrinsic nonlinearities in the system. Determining the degree of a system’s deviation from linearity is difficult [21, 22, 23], and the extent to which the approximation loses accuracy in nonlinear systems is unclear. Thus, although the Gaussian approximation offers a simple method to estimate information transmission, it remains an open question under what conditions this approximation is sufficiently accurate.

Until recently, addressing this question has been hard because there was no reliable benchmark for the exact information rate. Without a method to compute the true information rate of a non-Gaussian system, it is impossible to rigorously assess the accuracy of the Gaussian approximation. This gap was filled recently by the development of two independent methods [5, 4] for computing the information rate accurately even in systems that significantly deviate from Gaussian behavior. Here we leverage one of these methods: Path Weight Sampling (PWS) [5]. This Monte Carlo technique is an exact method for calculating the mutual information rate for a wide range of stochastic models. Using PWS calculations, we can directly evaluate the accuracy of the Gaussian approximation in models that exhibit explicit non-Gaussian features, and study the approximation’s robustness in typical applications.

In this article, we probe the accuracy of the Gaussian approximation for the information rate in two case studies. The first focuses on Markov jump processes, where the statistics are non-Gaussian due to their discrete nature. We first show that the Gaussian approximation fails to accurately estimate the mutual information rate in this case, even in the regime where the statistics are nearly Gaussian [4, 5]. We then apply a recently developed “discrete approximation” [4] and show it to be much more accurate. We briefly discuss this surprising failure of the Gaussian approximation, which has been more extensively analyzed in Ref. [24]. Finally, we still find regimes where neither of the approximations are highly accurate.

The second case study examines a continuous diffusive process with a nonlinear transfer function. We demonstrate how intrinsic nonlinearity can cause significant deviations between the Gaussian approximation and the true mutual information rate. By varying the degree of nonlinearity as well as the system’s response time, we provide a comprehensive quantitative understanding of the Gaussian approximation’s limitations in nonlinear systems. Additionally, we show that for such systems, the Gaussian approximation differs significantly when derived from empirical correlation functions compared to when it is analytically obtained from the nonlinear model, highlighting that the correct application of the approximation is important.

Our work translates into concrete recommendations on when to use which method for the computation of the information rate. It therefore enables researchers to more confidently determine when a simpler approximate method is sufficient, or when a more sophisticated method like PWS [5] or the method developed by Moor and Zechner [4] should be used.

II Methods

II.1 The mutual information rate

The mutual information between two random variables $S$ and $X$ is defined as

[TABLE]

or, equivalently, using Shannon entropies

[TABLE]

In the context of a noisy communication channel, $S$ and $X$ represent the messages at the sending and receiving end, respectively. Then, $I(S,X)$ is the amount of information about $S$ that is communicated when only $X$ is received. If $S$ can be perfectly reconstructed from $X$, the mutual information reaches its maximal value $I(S,X)=H(S)$. On the contrary, if $S$ and $X$ are independent, $I(S,X)=0$. The mutual information thus is always non-negative and quantifies the degree of statistical dependence between two random variables.

For systems that continuously transmit information over time, this concept must be extended to trajectories $\bm{S}_{T}=\{S(t)\mid t\in[0,T]\}$ and $\bm{X}_{T}=\{X(t)\mid t\in[0,T]\}$. The mutual information between trajectories is defined analogously as

[TABLE]

where the expected value is taken with respect to the full joint probability of both trajectories. This quantity can be interpreted as the total information that is communicated over the time interval $[0,T]$.

Note that the total amount of information communicated over the time-interval $[0,T]$ is not directly related to the instantaneous mutual information $I(S(t),X(t))$ at any instant $0\leq t\leq T$. This is because auto-correlations within the input or output sequences reduce the amount of new information transmitted in subsequent measurements. Moreover, information can be encoded in temporal features of the trajectories, which cannot be captured by an instantaneous information measure. Therefore, as previously pointed out [25, 26], the instantaneous mutual information $I(S(t),X(t))$ for any given $t$ does not provide a meaningful measure of information transmission. To correctly quantify the amount of information transmitted per unit time, we must consider entire trajectories.

For that reason, the mutual information rate is defined via the trajectory mutual information. Let the input and output of a system be given by two continuous-time stochastic processes $\mathcal{S}=\{S(t)\mid t\in\mathbb{R}\}$ and $\mathcal{X}=\{X(t)\mid t\in\mathbb{R}\}$. Then, for stationary systems, the mutual information rate between $\mathcal{S}$ and $\mathcal{X}$ is

[TABLE]

and quantifies the amount of information that can reliably be transmitted per unit time. The mutual information rate therefore represents an excellent performance measure for information processing systems.

In summary, the mutual information rate is the crucial performance metric for stochastic information processing systems. However, its information-theoretic definition does not translate into an obvious scheme for computing it. As a result, various methods have been developed to compute or approximate the mutual information rate.

II.2 Gaussian Approximation

To significantly simplify the computation of the information rate it is often assumed that the input and output trajectories obey stationary Gaussian statistics. Under this assumption Eq. 3 simplifies to

[TABLE]

where $|\bm{C}_{ss}|$ and $|\bm{C}_{xx}|$ are the determinants of the covariance matrices of the respective trajectories $\bm{S}_{T}$ and $\bm{X}_{T}$, and

[TABLE]

is the covariance matrix of their joint distribution.

In the limit that the trajectory length $N=T/\Delta$, with the discretization $\Delta$, becomes infinitely long ($N\to\infty$) and continuous ($\Delta\to 0$), the information rate as defined in Eq. 4 can be expressed in terms of the power spectral densities, or power spectra, of the Gaussian processes $\mathcal{S}$ and $\mathcal{X}$ [1, 2]:

[TABLE]

Here, $S_{ss}(\omega)$ and $S_{xx}(\omega)$ respectively are the power spectra of trajectories generated by $\mathcal{S}$ and $\mathcal{X}$, and $S_{sx}(\omega)$ is their cross-spectrum. The fraction

[TABLE]

is known as the coherence, describing the distribution of power transfer between $\mathcal{S}$ and $\mathcal{X}$ over the frequency $\omega$.

For systems that are neither Gaussian nor linear, there are two ways to still obtain an approximate Gaussian information rate. The first is to directly measure two-point correlation functions from data or simulations, and use these to numerically obtain the power spectra in Eq. 7. The second is to use Van Kampen’s linear noise approximation (LNA) and approximate the dynamics of the system to first order around a fixed point [27] which makes it possible to analytically obtain approximate power spectra, see also Appendix A. In this work, we will analyze both of these methods.

II.3 Path Weight Sampling for diffusive systems

To evaluate the accuracy of the Gaussian information rate for non-Gaussian systems, an exact method for determining the true information rate is required. Recently, a method called Path Weight Sampling (PWS) was developed, which computes the exact mutual information rate using Monte Carlo techniques without relying on approximations [5].

In Ref. [5], PWS was introduced as a computational framework for calculating the mutual information rate in systems governed by master equations. Master equations provide an exact stochastic description of continuous-time processes with discrete state spaces, commonly used in models ranging from biochemical signaling networks to population dynamics. However, many systems are not described by discrete state spaces and instead require a stochastic description based on diffusion processes or other stochastic models. Fortunately, PWS is not restricted to systems described by master equations and can be extended to a variety of stochastic models.

In general, PWS can be applied to any system that meets the following conditions: (i) sampling from the input distribution $\operatorname{\mathrm{P}}(\bm{s}_{T})$ is straightforward, (ii) sampling from the conditional output distribution $\operatorname{\mathrm{P}}(\bm{x}_{T}\mid\bm{s}_{T})$ is straightforward, and (iii) the logarithm of the conditional probability density $\ln\operatorname{\mathrm{P}}(\bm{x}_{T}\mid\bm{s}_{T})$, referred to as the path weight, can be evaluated efficiently. For any stochastic model that satisfies these three criteria, the PWS computation proceeds similarly to systems governed by master equations.

Briefly, PWS computes the trajectory Mutual Information using a Monte Carlo estimate of Eq. 3

[TABLE]

where $\bm{s}^{1}_{T},\ldots,\bm{s}^{N}_{T}$ are independently drawn from $\operatorname{\mathrm{P}}(\bm{s}_{T})$, and each $\bm{x}^{i}_{T}$ is drawn from $\operatorname{\mathrm{P}}(\bm{x}_{T}\mid\bm{s}^{i}_{T})$. As $N\to\infty$, this expression converges to the mutual information $I(\bm{S}_{T},\bm{X}_{T})$. In Eq. 9, the term $\ln\operatorname{\mathrm{P}}(\bm{x}_{T}\mid\bm{s}_{T})$ can be evaluated directly (per criterion iii), but the marginal probability $\operatorname{\mathrm{P}}(\bm{x}_{T})$ has to be computed separately for each output trajectory $\bm{x}^{i}_{T}$. Typically, this has to be done numerically via marginalization, i.e., by computing the path integral

[TABLE]

using Monte Carlo integration. Evaluating the marginalization integral efficiently is essential for computing the mutual information using PWS and discussed in detail in Ref. [5]. In summary, PWS is a generic framework that can be used beyond systems defined by a master equation as long as a suitable generative model satisfying the three conditions above is available.

For this study, we extended PWS to compute the mutual information rate for systems with diffusive dynamics, described by Langevin equations. For such systems, the aforementioned conditions are inherently fulfilled and PWS can be applied. Specifically, in a Langevin system, both the input $S(t)$ and the output $X(t)$ are stochastic processes given by the solution to a stochastic differential equation (SDE). Using stochastic integration schemes like the Euler-Mayurama method, we can straightforwardly generate realizations $s(t)$ and $x(t)$ from the corresponding stochastic process. These realizations are naturally time-discretized with the integration time step $\Delta t$. For a time-discretized trajectory $\bm{x}=(x_{1},\ldots,x_{n})$, the path weight $\ln\operatorname{\mathrm{P}}(\bm{x}|\bm{s})$ is—up to a Gaussian normalization constant—given by the Onsager-Machlup action [28]

[TABLE]

where we used $\Delta x_{i}=x_{i+1}-x_{i}$, and $v_{i}=f(x_{i},s_{i})$ is the deterministic drift, and $\sigma(x_{i})$ represents the white noise amplitude. This expression captures the likelihood of a particular trajectory, given the stochastic dynamics of the system, and serves as the path weight in the PWS computation.

III Case Studies

To investigate the conditions under which the Gaussian approximation deviates from the exact mutual information rate, we conducted two case studies. In both studies we compare the Gaussian approximation against the exact mutual information rate, computed via PWS. In the first case study we focus on a discrete linear system which is inspired by minimal motifs of cellular signaling.

III.1 Discrete reaction system

We consider a simple linear reaction system of two species, $S$ and $X$, whose dynamics are governed by 4 reactions

[TABLE]

The reaction system is linear because each reaction has at most one reactant. The trajectories of $S$ and $X$ are correlated because the production rate of $X$ depends on the copy number of $S$, and therefore information is transferred from $S$ to $X$. This set of reactions can be interpreted as a simple motif for gene expression where $S$ is a transcription factor and $X$ represents the expressed protein. In steady state, the mean copy numbers are given by $\bar{s}=\kappa\lambda^{-1}$ and $\bar{x}=\bar{s}\rho\mu^{-1}$.

The exact stochastic dynamics of this reaction system can be expressed by the chemical master equation [27]. This equation describes the time-evolution of the probability distribution over the possible copy numbers of species $S$ and $X$, capturing the noise from the discrete chemical reaction events. From this description, we can obtain the mutual information rate from $S$ to $X$ without approximations using PWS [5].

While the chemical master equation is an exact representation of the reaction system, for large copy numbers the stochastic dynamics are well approximated by a Gaussian model. The resulting Langevin equations can be systematically derived from the master equation using the LNA which yields

[TABLE]

where $s$ and $x$ are continuous variables representing the copy numbers of $S$ and $X$, and $\eta_{s},\eta_{x}$ are independent delta-correlated white noise terms with $\langle\eta^{2}_{s}\rangle=2\lambda\bar{s}$ and $\langle\eta^{2}_{x}\rangle=2\mu\bar{x}$, see Appendix A.

The Gaussian approximation of the mutual information rate is derived from the LNA description. Using this framework, Tostevin and ten Wolde [1] found an analytical expression for the mutual information rate of the motif in units of nats $s^{-1}$:

[TABLE]

More recently, Moor and Zechner [4] have derived a different expression for the mutual information rate of this reaction system by analytically approximating the relevant filtering equation, which is derived from the master equation, thus recognizing the discreteness of the state space. This approach explicitly differentiates the contributions of individual reactions to the noise amplitude of each component, while the LNA lumps their contributions together. As we will discuss in more detail below, separately accounting for the noise from each reaction better captures the information transmitted via discrete systems, making this “discrete approximation” more accurate than the Gaussian approximation for this case study. Nevertheless, the result is still based on an approximation that is only expected to be accurate for large copy numbers. The expression for the mutual information rate in the discrete approximation is remarkably similar to the expression obtained using the Gaussian framework:

[TABLE]

Note that this equation only differs from the Gaussian expression in Eq. 14 by the additional factor 2 inside the square root.

The natural—but incorrect—expectation is that for large copy numbers both approximations converge to the true mutual information rate. However, the differences between Eqs. 14 and 15 already reveal that the two approximations do not converge. Indeed, previous work has shown that even in the limit of infinite copy numbers, the Gaussian approximation only yields a lower bound which is not tight [5, 4]. We further discuss this point below.

We compared both approximations against exact PWS simulations for different parameter choices. In Fig. 1a, we vary the mean copy number of the readout, $\bar{x}$, by varying its synthesis rate $\rho$ and compute the mutual information rate using both approximations as well as PWS, while keeping the input copy number constant at $\bar{s}=100$. We observe that the Gaussian approximation via the LNA [Eq. 14] consistently underestimates the mutual information rate. This confirms that even when $\bar{s}$ and $\bar{x}$ are large, the Gaussian approximation only yields a lower bound to the information rate of the discrete linear system. In contrast, the discrete approximation [Eq. 15] coincides with the true mutual information rate obtained from PWS simulations over all output copy numbers $\bar{x}$, even for $\bar{x}\ll 1$.

In Fig. 1b, instead of varying the copy number of the output, we vary the copy number of the input by varying the production rate $\kappa$. Note that both the Gaussian approximation and the discrete approximation are independent of $\kappa$. Yet, we observe that the true mutual information rate is not. For sufficiently large input copy numbers the discrete approximation coincides with the true information rate while the Gaussian information rate remains only a lower bound. Thus, the discrete approximation is highly accurate for $\bar{s}\geq 10$. For small $\kappa$ where $\bar{s}<10$, we find that the mutual information rate deviates from both the LNA as well as the discrete approximation. Surprisingly, we find an optimal value of $\kappa$ for which the mutual information rate is maximized and exceeds both approximations. This implies that at low input copy numbers, the system is able to extract additional information from the discrete input trajectories, which is not accounted for by either of the approximations.

In all cases, we found that the Gaussian approximation deviates significantly from the true information rate for this discrete system. Seemingly paradoxically, the Gaussian approximation based on the LNA does not converge to the true information rate at high copy numbers, even though the LNA approximates the stochastic dynamics extremely well in this regime. In contrast, the discrete approximation from Moor and Zechner [4] does not suffer from this issue. It has been shown that, generally, the Gaussian approximation is a lower bound on the discrete approximation, prompting the question of which features of the trajectories are not captured by the Gaussian approximation.

In recent work [24], we found that the root cause for the deviations of the Gaussian approximation lies in how the LNA approximates the reaction noise in the chemical master equation. While the dynamics of the chemical master equation give rise to discrete sample paths, i.e., piece-wise constant trajectories connected by instantaneous discontinuous jumps, the LNA approximation yields continuous stochastic trajectories. Our results imply that a discrete sample path of $X$ carries more information about $S$ than the corresponding continuous sample path $x(t)$ would carry about $s(t)$ in the LNA. In Ref. [24] we show that this is ultimately due to the fact that in the discrete system, each reaction event is unambiguously recorded in the $X$ trajectories and thus different reactions modifying the same species can be distinguished. In contrast, in the continuous LNA description, all reactions that modify $X$ contribute to the noise term $\eta_{x}(t)$ in Eq. 13 but their contributions are lumped together and therefore cannot be distinguished from an observed $x(t)$-trajectory. Specifically, note that for the motif studied here, only the production reaction $S\to S+X$ conveys information. The decay reaction of the output $X\to\emptyset$ does not carry information on the input fluctuations since its propensity is independent of the input. Yet, it contributes to the overall fluctuations in the output. The Gaussian approximation only considers the total fluctuations in the output, while the discrete approximation correctly distinguishes between the fluctuations induced from production events and decay events. Therefore, the Gaussian approximation consistently underestimates the true information transmission, whereas the discrete approximation does not incur this systematic error. This subtle point is reflected in the difference between Eqs. 14 and 15.

III.2 Nonlinear continuous system

Next, we study a nonlinear variant of the reaction system above. In contrast to the previous case study, we deliberately avoid using discrete dynamics, as we already observed that the Gaussian approximation is generally inaccurate in such systems. Instead, we focus solely on continuous Langevin dynamics to explore how an explicitly nonlinear input-output mapping affects the accuracy of the inherently linear Gaussian approximation. We hypothesize that the accuracy of the Gaussian approximation will deteriorate as the degree of nonlinearity increases. To test this hypothesis, we analyze a simple Langevin system with adjustable nonlinearity.

The system is defined by two coupled Langevin equations, one that describes the input, and one that describes the output. The stochastic dynamics of the input $s(t)$ are given by Eq. 12. The output dynamics of $x(t)$ are given by

[TABLE]

with the Hill function

[TABLE]

This function serves as a tuneable non-linearity with the Hill coefficient $n$. As $n\to 0$, the Hill function approaches a shallow linear mapping, while for large $n$, it becomes sigmoidal and highly non-linear. As $n\to\infty$, $a(s)$ approaches the unit step function centered at $s=K$. The so-called static input-output relation specifies the mean output $\bar{x}(s)$ for a static input signal $s$ and is given by $\bar{x}(s)=\rho a(s)/\mu$. The gain of this system is then defined as the slope of this relation at $s=\bar{s}$, i.e.,

[TABLE]

as derived in Section A.3. Importantly, for $\bar{s}=K$, the gain of the system is directly proportional to the Hill coefficient $n$, i.e., the gain is directly coupled to the degree of nonlinearity.

The information rate quantifies the fidelity of the input–output mapping of the system and in dynamical systems this mapping is inherently limited by the response time scale. Figure 2 shows for various response time scales how the average output $\bar{x}(t)$ at a given time depends on the input $s(t)$ at that same time which is known as the dynamical input-output relation [31]. The solid black curve represents the static input-output relation which is determined by the instantaneous transfer function $a(s)$. As the output responds more slowly to the input, the temporal input fluctuations are averaged out increasingly, leading to a shallower response for increasing $\tau_{x}$. Indeed, for very slow time scales, the system is well approximated by a linear responder. Thus, we expect the Gaussian approximation of the information rate to be more accurate for slowly responding systems.

With the Langevin extension of PWS we can directly compute the mutual information rate of this nonlinear model whereas the Gaussian approximation can only be applied to linear systems. Therefore, to obtain the mutual information rate in the Gaussian approximation, we have to linearize the system. There are two approaches for linearizing the stochastic dynamics of this nonlinear system which result in different information estimates.

The first approach is to linearize Eq. 16 analytically via the LNA as shown in Section A.2. Within this approach, we can obtain an analytical expression for the information rate (see Section A.3),

[TABLE]

This LNA-based approach also yields a linearized dynamic input-output relation, shown as dashed lines in Fig. 2.

We observe that the linearized input-output relation closely matches the slope of the true nonlinear dynamical input-output relation at $s=\bar{s}=100$, but overall it does not correspond to a (least-squares) linear fit of the nonlinear dynamical input-output relation. For all values of $s$, the linearized input-output relation has a slope greater than or equal to the slope of the dynamical input-output relation. Empirically, the LNA thus seems to over-estimate the dynamic gain of the system. The reason may be that the LNA approximates the static input-relation (Fig. 2 black curve), and estimates the linearized dynamical input-output relation based on this static approximation only.

The second “empirical Gaussian” approach to linearize the nonlinear system avoids these issues. In this approach, we first numerically generate trajectories from the stochastic Eqs. 12 and 16 and use digital signal processing techniques to estimate the mutual information rate from the trajectories. We numerically estimate the (cross) power spectra of input and response using Welch’s method [32, see Ch. 11]. From the estimated spectral densities $\hat{S}_{\alpha\beta}(\omega)$ we compute the coherence

[TABLE]

which we use to obtain the Gaussian approximation of the mutual information rate directly using Eq. 7.

The empirical power spectra characterize the linear response of a system, but not in the same way as the LNA. While for linear systems the power spectra obtained via the LNA match the empirical power spectra [33], for a nonlinear system, the empirical power spectra and the coherence can differ from the corresponding LNA calculations. The two linearization approaches are thus not equivalent. We tested the accuracy of the Gaussian mutual information rate estimates using both linearization approaches to elucidate the differences in these approaches.

Figure 3 displays the mutual information rate obtained via two linearized approximations as well as the exact PWS result. We vary the gain $g$ and the response time-scale $\tau_{x}$, both of which significantly affect the shape of the dynamical input-output relationship. As expected, a larger gain or a faster response time lead to an increase in the mutual information rate. At large gain, the information rate naturally saturates as $a(s)$ approaches a step function. The saturation effect is clearly seen in the PWS results, and is found to be even more pronounced in the empirical Gaussian approximation. The LNA-based Gaussian approximation, however, shows no saturation. This highlights that the LNA directly linearizes the system at the level of the input-output mapping $a(s)$ which results in an approximation that is unaffected by the sigmoidal shape of $a(s)$. In contrast, the empirical approximation is affected by nonlinear saturation effects because it is computed from simulated trajectories. We thus see that both approximations yield substantially different results at large gain.

In Fig. 4 we compare the absolute deviation between the approximations and the PWS result. For small gain we see that both approximations are accurate which is not surprising since the nonlinearity is very weak in this regime. For large gain, the LNA-based approximation always overestimates the mutual information rate while the empirical Gaussian approximation always underestimates the rate. In both cases the systematic error decreases as the response timescale becomes slower. This reflects the fact that for slow responders, the dynamic input-output relation is more linear (Fig. 2) than for fast responders.

Additionally, we computed the relative deviation, see Fig. S1 in Appendix B. We find that in terms of relative error the curves for different response timescales largely overlap. In terms of relative approximation error, the gain, rather than response timescale, is the primary factor affecting the accuracy of the Gaussian approximation.

IV Discussion

We investigated the accuracy of the Gaussian approximation for the mutual information rate in two case studies, each highlighting a scenario where the approximation may be inaccurate. We were able to reliably quantify the inaccuracy in each case by comparing against the “ground truth” mutual information rate for these scenarios which was computed using a recently developed exact Monte Carlo technique called PWS [5].

We first considered linear discrete systems, which are relevant in biology due to the discrete nature of biochemical signaling networks. In our example, the Gaussian approximation cannot capture the full information rate, but only yields a lower bound. We show that a discrete approximation, developed by Moor and Zechner [4], is able to correctly estimate the mutual information rate of the network over a wide range of parameters. Since the Gaussian approximation captures the second moments of the discrete system, this finding demonstrates that a discrete system can transmit significantly more information than what would be inferred from its second moments alone. This perhaps surprising fact has been highlighted before [34, 4, 5], and it hinges on the use of a discrete reaction-based readout. As demonstrated in Ref. [24], the increased mutual information rate found for a discrete readout stems from the ability of unambiguously distinguishing between individual reaction events in the readout’s trajectory. However, it remains an open question whether biological (or other) signaling systems can effectively harness this additional information encoded in the discrete reaction events. For systems that cannot distinguish individual reaction events in downstream processing, the Gaussian framework might still accurately quantify the “accessible information”.

A notable new observation in our first case study is the deviation between the discrete approximation of the mutual information rate derived by Moor and Zechner [4] and the exact result obtained using PWS [5] for inputs with low copy number $\bar{s}$. In the discrete approximation, the mutual information rate is independent of the input copy number $\bar{s}$, but the PWS simulations show that at low copy numbers there is an optimal $s^{\star}\approx 1$ which maximizes the mutual information rate, see Fig. 1b. This surprising finding suggests that the information rate in discrete systems can be increased by reducing the copy number of the input sufficiently, such that it only switches between a few discrete input levels. Notably, in the reverse case—low output copy number $\bar{x}$ but large $\bar{s}$—the discrete approximation always remains accurate, see Fig. 1a. We leave a precise characterization of this finding for future work.

The second case study focused on a continuous but nonlinear system, where we demonstrated that the accuracy of the Gaussian approximation depends on the linearization method. Linearizing the underlying system dynamics directly via the LNA leads to an overestimation of the information rate, while estimating the system’s correlation functions empirically from data underestimates it. Regardless of the method, the Gaussian approximation is more accurate in terms of absolute deviation of the true information rate when the gain of the system is small and its response is slow compared to the timescale of the input fluctuations.

The result of our second case study—that the empirical Gaussian mutual information rate underestimates the true rate—is consistent with theoretical expectations. As shown by Mitra and Stark [23], and highlighted in [2], an empirical Gaussian estimate of the mutual information between a Gaussian input signal $S_{G}$ and a non-Gaussian output $X$ provides a lower bound on the channel capacity $C(S,X)=\max_{\operatorname{\mathrm{P}}(s)}I(S,X)$ (subject to a power constraint on $S$). Specifically, they show that $C(S,X)\geq I(S_{G};X)\geq I(S_{G};X_{G})$, where $(S_{G},X_{G})$ is a jointly Gaussian pair with the same covariance matrix as $(S_{G},X)$. For purely Gaussian systems like $(S_{G},X_{G})$, the mutual information calculated using Eq. 7 is exact and equal to the channel capacity. However, for systems that have a Gaussian input but are otherwise non-Gaussian, the mutual information is larger or equal than the corresponding Gaussian model with matching second moments, as evidenced in Fig. 4. In general, the empirical Gaussian approximation yields a lower bound on the mutual information of the nonlinear system with a Gaussian input signal, as well as a lower bound on the channel capacity of the nonlinear system 111Note that this argument does not apply to the Linear Noise Approximation (LNA). The bound specifically requires the Gaussian model to use the covariance of the full, original system. When the system is first linearized using the LNA, the resulting linear model does not retain the same covariance as the original nonlinear system. As a result, the mutual information rate calculated with the LNA is not a lower (nor an upper) bound on the true mutual information rate..

We can distill several concrete recommendations for the computation of the information rate from our analysis. For linear discrete systems, the Gaussian approximation yields a lower bound on the true information rate which may accurately quantify the information available to systems that cannot distinguish individual discrete events. Alternatively, the reaction-based discrete approximation by Moor and Zechner [4] is highly accurate, even when the copy number of the output is extremely small. However, when the copy number of the input becomes small ($\lesssim 10$), both approximations break down and one must use an exact method. Exact methods for obtaining the information rate of any stochastic reaction-based systems are PWS [5] or brute-force numerical integration of the stochastic filtering equation, as shown in [4]. For nonlinear continuous systems with small gain one can safely use the Gaussian approximation, either based on a linearization of the underlying dynamics or on empirically estimated correlation functions. Moreover, when the slowest input time-scale is more than a magnitude faster than the response time-scale, the non-linear response of the system is averaged out by the quick input fluctuations, and the Gaussian approximation yields accurate results. In this case, using a Gaussian approximation based on empirical correlation functions yields the most accurate result, and provides a lower bound for the mutual information. Finally, if the system is both highly nonlinear and has a fast response with respect to the input, one must resort to an exact method like PWS. We hope that our results will guide future research in determining the appropriate method for computing the mutual information rate.

Overall, our results greatly increase the usefulness of the Gaussian approximation for the information rate of non-Gaussian systems. The Gaussian approximation remains a useful method that can be applied directly and straightforwardly to experimental data. Here, we have quantified the prerequisites to safely use this approach. Moreover, we elucidate how an empirical Gaussian approximation constitutes a lower bound on the true information rate for systems with a sufficiently large input copy number.

Acknowledgements.

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 Gaussian approximation

Here we derive the analytical expressions for the Gaussian information rate of the networks considered in the main text. To this end we first discuss the dynamics of the input signal $S$ and its power spectrum. Then, we perform a linear approximation of the dynamics of the readout species $X$ and derive the approximate Gaussian information rate between $S$ and $X$ for the nonlinear network. Finally, we derive the Gaussian information rate of the linear network from our expression of the Gaussian information rate of the nonlinear network.

A.1 Signal

The input signal is generated by a birth-death process,

[TABLE]

Its dynamics in Langevin form are,

[TABLE]

yielding the steady state signal concentration $\bar{s}=\kappa/\lambda$. The independent Gaussian white noise process $\eta_{s}(t)$ summarizes all reactions that contribute to fluctuations in $S$. The strength of the noise term in steady state is

[TABLE]

The power spectral density, or power spectrum, of a stationary process $\mathcal{X}$ is defined as $S_{xx}(\omega)=\lim_{T\to\infty}\frac{1}{T}|\tilde{x}_{T}(\omega)|^{2}$, where $\tilde{x}(\omega)$ denotes the Fourier transform of $x(t)$. The power spectrum of a signal obeying Eq. 22 is thus given by

[TABLE]

A.2 Linear approximation

We now consider the readout $X$, which is produced via a nonlinear activation function $a(s)$:

[TABLE]

We define the activation level $a(s)$ to be a Hill function,

[TABLE]

Such a dependency, in which $K$ sets the concentration of $S$ at which the activation is half-maximal and $n$ sets the steepness, can for example arise from cooperativity between the signal molecules in activating the synthesis of $X$.

We have for the dynamics of $X$ in Langevin form

[TABLE]

with $a(s)$ given by Eq. 26. The steady state concentration of $X$ is given by $\bar{x}=\bar{a}\rho/\mu$, where we have defined the steady state activation level $\bar{a}=a(\bar{s})$. It is useful to determine the static gain of the network, which is defined as the change in the steady state of the output upon a change in the steady state of the signal:

[TABLE]

where we have defined the approximate linear activation rate

[TABLE]

and the steady state of the activation level is given by

[TABLE]

Generally, we assume that $K=\bar{s}$, which entails that in steady state the network is tuned to $\bar{a}=1/2$.

To compute the Gaussian information rate we approximate the dynamics of $X$ to first order around $\bar{x}$ via the classical linear noise approximation [27]. Within this approximation the dynamics of the deviation $\delta x(t)=x(t)-\bar{x}$ are,

[TABLE]

with the synthesis rate $r$ given by Eq. 29.

In the linear noise approximation the noise strength is a constant given by the noise strength at steady state,

[TABLE]

A.3 Information rate

Following Tostevin & Ten Wolde [1, 2], we can express the Gaussian information rate as follows,

[TABLE]

where $|K(\omega)|^{2}$ is the frequency dependent gain and $|N(\omega)|^{2}$ is the frequency dependent noise of the output process $\mathcal{X}$. If the intrinsic noise of the network is not correlated to the process that drives the signal, the power spectrum of the network output obeys the spectral addition rule [36]. In this case the frequency dependent gain and noise can be identified directly from the power spectrum of the output, because it takes the following form:

[TABLE]

For a species $X$ obeying Eq. 31, we have

[TABLE]

The Wiener Khinchin theorem states that the power spectrum of a stochastic process and its auto-correlation function are a Fourier transform pair. We thus obtain for the variance in the readout, substituting the frequency dependent gain and noise [Eq. 35] and the power spectrum of the signal [Eq. 24] in Eq. 34 and taking the inverse Fourier transform at $t=0$,

[TABLE]

where the signal variance equals its mean $\sigma_{s}^{2}=\bar{s}$, and the mean readout concentration sets the intrinsic noise $\sigma^{2}_{x|s}=\bar{x}$. We further have the static gain $g$ given by Eq. 28, and have defined the dynamic gain

[TABLE]

which is the slope of the mapping from the time-varying signal value $s(t)$ to the time-varying readout $x(t)$; for Gaussian systems $\left\langle x(t)|s(t)\right\rangle=\tilde{g}s(t)$ [2, 31, 19].

To solve the integral in Eq. 33 we exploit that

[TABLE]

Substituting the frequency dependent gain and noise given in Eq. 35 and the signal power spectrum of Eq. 24 in Eq. 33 and using Eq. 38 we obtain the information rate,

[TABLE]

where we used the noise strengths given in Eq. 23 and Eq. 32, we have the static gain $g$ of Eq. 28, and the synthesis rate $r$ of Eq. 29.

A.4 Linear network

To disambiguate differences in the information rate caused by the linear approximation of our nonlinear reaction network on the one hand and the Gaussian approximation of the underlying jump process on the other, we consider the information rate of a linear network. Any difference between the exact information rate and the Gaussian information rate must then be a result of the Gaussian approximation. To this end we use the same input signal [Eq. 22], but we consider a linear activation of the readout, i.e.

[TABLE]

such that the Langevin dynamics of $X$ are

[TABLE]

which yields the steady state concentration $\bar{x}=\bar{s}\rho/\mu$. For this linear readout, the static gain is simply set by the ratio of steady states of the input and the output, $g=\rho/\mu=\bar{x}/\bar{s}$. We can then obtain the information rate of this linear system by substitution of its static gain in Eq. 39, which yields

[TABLE]

This result is identical to that of Tostevin and ten Wolde [1, 2] (motif III).

Appendix B Relative deviation of the Gaussian approximation for a nonlinear system

Due to the definition of the mutual information, an absolute difference in information maps to a relative difference in the reduction of uncertainty. Therefore, Fig. 4 in Section III.2 focuses on the absolute deviation between the Gaussian approximation and the true mutual information. Nevertheless, the relative deviation of the Gaussian information rate from the true rate can still provide valuable insights, and we discuss it here.

In Fig. S1 we compare the relative deviation $R_{\text{Gaussian}}/R_{\text{PWS}}$ between the Gaussian approximation and the exact mutual information computed using PWS. We find that the relative deviation increases as the system gain increases, indicating that the Gaussian approximation also becomes relatively less accurate for larger gains. As already discussed above, the empirical Gaussian method consistently underestimates the true information rate, while the LNA-based approximation overestimates it.

Interestingly, we also observe that for the LNA approximation, at fast timescales the result is slightly more accurate, whereas the empirical Gaussian estimate is more accurate at slow timescales. We initially expected that in both cases slow timescales would yield better agreement with PWS, as the input-output dynamics are more linear for slow timescales, and thus better approximated by the Gaussian model. The fact that this is not the case for the LNA approximation is intriguing, indicating the need for further investigation into the interplay between timescales, system nonlinearity, and the LNA.