Across diverse biological systems—ranging from neural networks to intracellular signaling and genetic regulatory networks—the information about changes in the environment is frequently encoded in the full temporal dynamics of the network nodes. A pressing data-analysis challenge has thus been to efficiently estimate the amount of information that these dynamics convey from experimental data. Here we develop and evaluate decoding-based estimation methods to lower bound the mutual information about a finite set of inputs, encoded in single-cell high-dimensional time series data. For biological reaction networks governed by the chemical Master equation, we derive model-based information approximations and analytical upper bounds, against which we benchmark our proposed model-free decoding estimators. In contrast to the frequently-used k-nearest-neighbor estimator, decoding-based estimators robustly extract a large fraction of the available information from high-dimensional trajectories with a realistic number of data samples. We apply these estimators to previously published data on Erk and Ca2+ signaling in mammalian cells and to yeast stress-response, and find that substantial amount of information about environmental state can be encoded by non-trivial response statistics even in stationary signals. We argue that these single-cell, decoding-based information estimates, rather than the commonly-used tests for significant differences between selected population response statistics, provide a proper and unbiased measure for the performance of biological signaling networks.
For their survival, reproduction, and differentiation, cells depend on their ability to respond and adapt to continually changing environmental conditions. Environmental information must be sensed and often transduced to the nucleus, where an appropriate response is initiated, usually by selectively up- or down-regulating the expression levels of target genes. This information flow is mediated by biochemical reaction networks, in which concentrations of various signaling molecules code for different environmental states or different response programs. This map between environmental input or response output and the internal chemical state is, however, highly stochastic, because typical networks operate with small absolute copy numbers of signaling molecules [1]; moreover, different environments can be encoded by the same signaling molecule, by differentially regulating the dynamics of its concentration [2]. This raises two fundamental questions: first, how much information the cells could, even in principle, encode in the combinatorial and possibly time-varying concentrations of multiple signaling molecules and how such information could be plausibly read out during “downstream” processing; and second, how can we quantify, in an unbiased and model-free fashion, the amount of information available to the cells from limited experimental data. Information theory provides a framework within which the theoretical study of limits to communication as well as the empirical study of actual information flows can be addressed [3]. Applications of information theory to questions in biology and, in particular, neuroscience started already in the 1950s and continue to this day, with the main focus to understand how—and with what accuracy—neural activity encodes information about the environment [4–6]. Applications of analogous techniques to biochemical signaling only started recently and represent an active area of research at the interface of physics, biology, statistics, and engineering [7–10]. Recent theoretical work analyzed the reliability of information transmission through specific reaction systems in the presence of molecular noise, e. g. , during ligand binding [11], in chemotaxis [12], gene regulation [13–19], biochemical signaling networks [20], etc. , and asked how such transmission can be maximized by tuning the reaction rates. Generally, these studies focused on steady state, by considering the information encoded in a single temporal snapshot of the reaction network at equilibrium given the input signals. Rigorous extensions to dynamical signals have been either rare and only possible for simple cases, like the BIND channel [11], or required specific operating regimes that permitted linearization and Gaussianity assumptions [12,21,22]. At its core, the analysis of signal transduction through nonlinear noisy chemical systems requires one to have control over the distribution of concentration trajectories given the (possibly) time-varying inputs; even if it were possible to calculate this distribution in principle, the curse of dimensionality puts strong limits to the manipulations required to compute the information transmission. Consequently, problems of this kind are currently considered intractable in their full generality. Empirical estimates of information transmission in biochemical networks similarly focused on the steady state [23,24], or considered only specific, hand-picked dynamical features, such as the amplitude or the frequency of the response, as information carriers [25]. Recent developments of fluorescent reporters and microfluidics have enabled the characterization of dynamical responses at a single cell resolution using large (> 104) numbers of sampled response trajectories, thereby permitting direct information estimates using generic estimators like the k-nearest-neighbors (knn) [26]. Existing approaches, however, suffer from severe limitations: they still require a prohibitive number of samples, especially when the response is distributed over multiple chemical species; or they necessitate uncontrolled assumptions about trajectory features that are supposed to be “relevant”. We recently proposed and applied decoding-based information estimators [27] as an alternative that draws on the past experiences in neuroscience [28–30] to dissect the yeast stress-response network. In this paper we provide a detailed account of the new methodology, show that it alleviates the most pressing problems of existing approaches, and benchmark it against synthetic and real data.
At their core, cellular processes consist of networks of chemical reactions. A chemical reaction network consists of a set of m molecular species {X˜1, X˜2, …, X˜m} that interact through K coupled reactions of the form: ν1k′X˜1+…+νmk′X˜m→θkν1k″X˜1+…+νmk″X˜m, k=1, …, K (1) where ν1k′, …, νmk′ and ν1k″, …, νmk″ are coefficients that determine how many molecules of each species are consumed and produced in the k-th reaction. θ1…θk∈R+ determine the rates at which the reactions occur and depend on binding affinities of chemical species, temperature and possibly the external conditions. If we assume that the system is well-stirred, in thermal equilibrium and the reaction volume is constant, it can be shown that the probability that a reaction of type k takes place in an infinitesimal time interval [t, t + dt] can be written as ak (x˜, θk) dt=θkgk (x˜) dt, where x˜=[x˜1, …, x˜m]T∈N0m contains the amounts of molecules of the m species that are present in the system at time t, and gk (x˜) =∏i=1m (x˜iνik′) counts all possibilities of choosing the required reaction molecules out of all available molecules [31,32]. θk is a constant that depends on the physical characteristics of the cell but also on the environmental conditions. Let us denote the probability that x˜ molecules of the m species are present in the system at time t∈R+ by p (x˜, t) and define the stoichiometric change vectors νk=[ν1k, …, νmk]T∈Zm, k=1, …, K, as the net changes in the amount of molecules in the reactions, i. e. νik=νik″−νik′, i=1, …, m, k=1, …, K. Then it can be shown [32] that the chemical master equation (CME) can be written as: p˙ (x˜, t) =−p (x˜, t) ∑k=1Kak (x˜, θk) +∑k=1Kp (x˜−νk, t) ak (x˜−νk, θk), (2) or in a more compact form [32] p˙ (t) =Mp (t), (3) where p (t) is a vector with components p (x˜, t), which is, in principle, infinite dimensional, and M contains the transition rates between all possible states, e. g. the transition rate from state x˜k′=x˜−νk to state x˜ is given by Mx˜, x˜k′=ak (x˜k′, θk) −δx˜, x˜k′∑qaq (x˜, θq), (4) where δ is the Kronecker delta, which is 1 when x˜=x˜k′ and 0 otherwise. The CME given in Eq (3) is an instance of a continuous-time discrete-state-space Markov Chain for a random process X that can be solved exactly only for a few simple cases. It is nevertheless possible to efficiently generate samples x of the random process X, which we will refer to as “trajectories” or “paths”, for a selected time interval, t ∈ [0, T], according to the correct probability distribution p, by the Stochastic Simulation Algorithm (SSA, or the Gillespie algorithm) [33]. To study information transmission through the biochemical networks described by the CME, we need to define the input and output signals. In the simplest setup considered here, the input U is a discrete random variable that can take on one of the q ≥ 2 possible values, U ∈ {u (1), u (2), …, u (q) }. Each input in general corresponds to a distinct set of reaction rate constants θ, but in models of real biological networks, changing input often modulates only one or a few rates in the system, e. g. , by representing the change in a key external ligand concentration, receptor activity, etc. Changes in the input are reflected in changes in the resulting trajectories of chemical species amounts, x. Typically, only a subset of chemical species could be considered as biologically-relevant “outputs” that encode the information about the environmental change: this would correspond to marginalizing p in Eq (3) over the unobserved (non-output) chemical species for the purposes of information transmission. While this is an interesting theoretical problem in its own right, here we work with simple toy examples where the output will be the trajectory, x, over the complete state space, i. e. , we assume that all chemical species in the reaction network can be fully and perfectly observed. As we explain below, this allows us to define and compute the mutual information between a discrete input, U, and the output random process X given by the CME in a straightforward fashion. We later show that this computation can be carried out also when the continuous-time process X is sampled at uniform discrete times, as would be the case with experimental measurements. Information theory introduces the mutual information as the measure of fidelity by which changes in one random variable, e. g. , the input U, can effect changes in another random variable, e. g. , X. In this sense, mutual information is simply a measure of statistical dependency (i. e. , any correlation, be it linear or not) between U and X, and can thus be written as a functional of the joint probability density function p (x, u): I (X; U) =∫X∫Up (x, u) log2 (p (x, u) pX (x) pU (u) ) dudx (5) where pU and pX are the marginal density functions for U and X, respectively, and we have generically written u and x as continuous variables; if they are discrete, integral signs are replaced by summations over the support for the corresponding probability distributions, as appropriate. Mutual information is a non-negative symmetric quantity that is measured in bits, and is zero only if X and U are statistically independent. When studying information transmission through a channel U → X specified by p (x|u), for which U serve as inputs drawn from an input distribution pU (u), it is common to rewrite Eq (5) as I (X; U) =H (U) −H (U|X) =H (X) −H (X|U), (6) where H (X) is the differential entropy of X (and analogously for H (U) ), defined as H (X) =−∫XpX (x) log2pX (x) dx, (7) and the conditional entropy, H (X|U), is H (X|U) =∫UH (X|u) pU (u) du=−∫U∫XpU (u) p (x|u) log2p (x|u) dxdu. (8) Eq (6) can be interpreted in two ways: information is either the difference between the total variability in the repertoire of responses X that the biochemical network can generate (measured by the response entropy, H (X) ) and the average variability at fixed input that is due to noise in the network (measured by the noise entropy, H (X|U) ); alternatively, information is also the entropy of the inputs, H (U), minus equivocation H (U|X), or the average uncertainty in what input was sent given that a particular response was observed. These interpretations make explicit the dependence of information both on the properties of the channel (the biochemical reaction network), as well as on the distribution of signals pU that the network receives. In this work, we will consider discrete inputs and will assume uniform pU. It is, however, also possible to compute the channel capacity C by maximizing the information flow at given p (x|u) over all possible input distributions, C=maxpUI (X; U); (9) Shannon’s classic work then proves that error-free transmission at rates higher than those given by capacity is impossible, while error-free transmission at rates below capacity can be achieved with the optimal use of the channel. Contrary to engineering, where the focus is on finding encoding and decoding schemes that best utilize a given channel, in biophysics and systems biology mutual information is used as a tool to quantify the limits to biological signal processing due to noise without needing to make assumptions about possible biochemical encoding and decoding mechanisms. The setup we consider here is one in which inputs U are drawn independently from a uniform distribution and change rarely, i. e. , at a rate that is much lower than the (inverse) timescale on which the reaction network in Eq (1) relaxes to its steady state. We assume that after an input change, we observe a fixed-time segment of the complete network dynamics, x, which is a sample path in m-dimensional discrete space, making direct calculation of information, I (X; U), by integrating / summing over all possible trajectories as implied by Eq (5) intractable. We will nevertheless show that estimates of exact information are possible if the reaction network is known, by explicitly using the transition matrix M of the Markov Chain from Eq (3) and generating exact sample paths, that is, realizations of X, using SSA. We call this model-based approach exact Monte Carlo approximation and contrast it to uncontrolled model-free estimations such as those obtained by using Gaussian approximations or k-nearest-neighbors methodology. We then introduce various decoding estimators and establish a hierarchy through which these estimates upper and lower-bound the true information, as shown in Fig 1. In the absence of a full stochastic model for the biochemical reaction network, mutual information estimation is tractable only if we make assumptions about the distribution of response trajectories given the input. We briefly summarize two approaches below: in the first, k-nearest-neighbor procedure, the space in which the response trajectories are embedded is assumed to have a particular metric; in the second, Gaussian approximation, we assume a particularly tractable functional form for the channel, p (x|u). Here and in the next section we introduce a class of decoding-based calculations that lower-bound the exact information, I (X; U), and can tractably be used as information estimators over realistically-sized data sets. Let D consist of a set of N labeled paths, typically represented in discretely sampled time, D={ (u1, x1), (u2, x2), …, (uN, xN) }, where ui and xi, for i = 1, …, N, are realizations of the random variables U ∈ {u (1), …, u (q) } and X∈Rm×d, respectively. Here, D can represent either real data (typically containing N ∼ 102 − 103 trajectories) in case of model-free information estimates, or trajectories generated by exact simulation algorithms (in which case the sample size, N, is not limiting) from the full specification of the biochemical reaction network in case of model-based approximations. The procedure of estimating the input u^ from x, such that the estimated u^ is “as close as possible” to true u for a given trajectory x, is known as decoding in information theory and neuroscience, and can equivalently be viewed as a classification task in machine learning or as an inference task in statistics. This procedure is implemented by a decoding function, u^=Fω (x); (22) F is typically parametrized by parameters ω that need to be learned from data for model-free approaches, or derived from biochemical reaction network specification in case of model-based approaches. F assigns to every xi in the dataset a corresponding “decode” u^i from the same space over which the random variable U is defined; formally, these decodes are instances of a new random variable U^. The key idea of using decoding for information estimation starts with the observation that random variables U→X→TdX→FωU^, (23) where Td represents time discretization, form a Markov chain. In other words, the distribution of U^ is conditionally independent of U and only depends on X, p (u^|x, u) =p (u^|x), and so p (u^, x, u) =pU (u) p (x|u) p (u^|x). (24) The data processing inequality [43] can be used to further extend the bounds in Eq (18): I (U; U^) ≤Iexact (U; X) ≤Iexact* (U; X), (25) where equality between the first two terms holds only if I (U; X|U^) =0. Consequently, I (U; U^) is a lower bound to the information between trajectories X and the input U [44]. Note that analogous reasoning holds for decoding directly from continuous-time trajectories X. Better decoders which increase the correspondence between the true inputs and the corresponding decoded inputs will typically provide a tighter lower bound on the information. To compute the information lower bound, we apply the decoding function to each trajectory in D in model-based approximations or to each trajectory in the testing dataset for model-free estimators that need to be learned over training data first. We subsequently construct a q × q confusion matrix, also known as an error matrix, where each element ϵij counts the fraction of realizations of x generated by an input u = u (i) that decode into u^=u (j). This matrix provides an empirical estimate of the probability distribution p (u^, u), which can thus be used to compute the information estimate: I (U^; U) =∑u, u^p (u^, u) log2p (u^, u) pU (u) pU^ (u^) ≈∑i=1q∑j=1qϵijlog2ϵij (∑kϵkj) (∑lϵil), (26) Crucially, in this estimation O (N) data points are used to empirically estimate the elements of a q × q matrix ϵ, and information estimation involves a tractable summation over these matrix elements; in contrast, direct estimates of I (U; X) would involve an intractable summation over (vastly undersampled) space for X. For typical applications where q is small, decoding thus provides an essential dimensionality reduction prior to information estimation: in a simple but biologically relevant case of two distinct stimuli (q = 2), information estimation only requires us to empirically construct a 2 × 2 confusion matrix. If required, one can apply well-known debiasing techniques for larger q [5].
We start by considering three simple chemical reaction networks for which we can obtain exact information values using the model-based approach outlined in Section Exact information calculations for fully observed reaction networks. This will allow us to precisely assess the performance of decoding-based model-free estimates, and systematically study the effects of time discretization, the number of sample trajectories, and the number of distinct discrete inputs, q. The three examples are all instances of a simple molecular birth-death process, where molecules of X˜ are created and destroyed with rates α and β, respectively: →α (U) X˜→β (U) ⌀. (36) The reaction rates, α and β, will depend in various ways on the input, U, and possibly time, as specified below. Given an initial condition, x (t = 0), the production and degradation reactions generate continuous-time stochastic trajectories, x (t), recording the number of molecules of X˜ at every time t ∈ [0, T], according to the Chemical Master Eq (3). These trajectories, or their discretized representations, are considered as the “outputs” of the example reaction networks, defining the mutual information I (X; U) that we wish to compute. In all three examples we start with the simplest case, where the random variable U can only take on two possible values, u (1) and u (2), with equal probability, pU (u (1) ) = pU (u (2) ) = 0. 5. Example 1. In this case, x (t = 0) = 0, β = 0. 01, independent of the input U, and the production rate depends on the input as α (u (1) ) = 0. 1, α (u (2) ) = 0. 07. Here, the steady state is given by Poisson distribution with mean number of molecules 〈x (t → ∞) 〉 = α/β. Steady-state is approached exponentially with the timescale that is the inverse of the degradation rate, β−1. These dynamics stylize a class of frequently observed biochemical responses where the steady-state mean expression level encodes the relevant input value. Even if the stochastic trajectories for the two possible inputs are noisy as shown in Fig 2A, we expect that the mutual information will climb quickly with the duration of the trajectory, T, since (especially in steady state) more samples provide direct evidence about the relevant input already at the level of the mean trajectories. Example 2. In this case, x (t = 0) = 0, β = 0. 01, independent of the input U, and the production rate depends on the input as α (u (1), t) = 0. 1, α (u (2), t) = 0. 05 for all t < 1000, while for t ≥ 1000 the production rate is very small and independent of input, α (u, t) = 5 ⋅ 10−4. In the early period, this network approaches input-dependent steady state with means whose differences are larger than in Example 1, but the difference decays away for t > 1000 as the network settles towards vanishingly small activity for both inputs, as shown in Fig 2B. These dynamics stylize a class of transient biochemical responses that are adapted away even if the input state persists. In this case, lengthening the observation window T will not provide significant increases in information. Example 3. In this case, x (t = 0) = 10. All reaction rates depend on the input, α (u (1) ) = 0. 1, α (u (2) ) = 0. 05, β (u (1) ) = 0. 01, β (u (2) ) = 0. 005, and are chosen so that the mean 〈x (t) 〉 = 10 is constant across time and equal for both conditions, as shown in Fig 2C. In this difficult case, inputs cannot be decoded at the level of mean responses but require sensitivity to at least second-order statistics of the trajectories. Specifically, signatures of the input are present in the autocorrelation function for x: the timescale of fluctuations and mean-reversion is two-fold faster for u1 than u2. While this case is not frequently observed in biological systems, it represents a scenario where, by construction, no information about the input is present at the level of single concentration values and having access to the trajectories is essential. Because there is no difference in the mean response, we expect linear decoding methods to provide zero bits of information about the input. This case is also interesting because of the recent focus on pulsatile stationary-state dynamics in biochemical networks [53]. These pulses, reported for transcription factors such as Msn2, NF-κB, p53, etc. , occur stochastically and, when averaged over a population of desynchronized cells, can yield a flat and featureless mean response. Information about the stimulus could, nevertheless, be encoded in either the frequency, amplitude, or other shape parameters of the pulses. While a generative description of such pulsatile dynamics goes beyond a birth-death process considered here, from the viewpoint of decoding, both pulsatile signaling and our example present an analogous problem, where the mean response is not informative about the applied input. Before proceeding, we note that our examples are not intended to be realistic models of intracellular biochemical networks, but are chosen here for their simplicity and analytical tractability, in order to benchmark model-free estimators against known “gold truth” standard. In particular, while our examples include intrinsic noise due to the stochasticity of biochemical reactions at low concentration, they do not include extrinsic noise or cell-to-cell variability which, in some systems, is known to importantly or even dominantly contribute to the total variability in the response [26,54]. The presence of such additional sources of variation by no means makes the model-free estimators inapplicable, as we show in S1 Fig where we study estimator performance in the simplest Example 1 model that includes cell-to-cell variability; it solely prevents us from comparing their performance to a tractably-computable MAP decoder result. To illustrate the use of our estimators in a realistic context, we analyzed data from two previously published papers. The first paper focused on the representation of environmental stress in the nuclear localization dynamics of several transcription factors (here we focus on data for Msn2, Dot6, and Sfp1) in budding yeast [27]. The second paper studied information transmission in biochemical signaling networks in mammalian cells (here we focus on data for ERK and Ca2+) [26]. In both cases, single-cell trajectory data were collected in hundreds or thousands of single cells sampled at sufficient resolution to represent the trajectories discretized at tens to hundreds of timepoints. Similarly, both papers estimate the information transmission in trajectories about a discrete number of environmental conditions: Ref [27] uses the linear SVM approach presented here, while Ref [26] uses the knn estimator. This makes the two datasets perfectly suited for estimator comparisons. We further note that in both datasets the trajectories can be divided into two response periods: the early “transient” response period when the external condition changes, and the late “near steady-state” response period. Typically, the transient dynamics exhibit clear differences in the trajectory means between various conditions, reminiscent of our Example 1 or early Example 2; in contrast, in the late period the response may have been adapted away, or the stimulus could be encoded only in higher-order statistics of the traces, reminiscent of the late period in Example 2 or Example 3. Fig 8 shows the raw data and summarizes our estimation results for the early and late response periods for the three translocating factors in yeast that report on the change from 2% glucose rich medium to 0. 1% glucose poor stress medium. Fig 9 similarly shows the raw data and estimation results for the early and late response periods for the signaling molecules in mammalian cells responding to multilevel inputs. Consistent with the published report [27], transient response in yeast nuclear localization signal can be decoded well with the linear SVM estimator that yields about 0. 6 bits of information per gene about the external condition. Kernelized SVM outperforms the linear method slightly by extracting an extra 0. 1-0. 2 bits of information, while knn underperforms the linear method significantly for Msn2 and Dot6 (but not for Sfp1). The Gaussian decoder estimate shows a mixed performance and the neural network estimate is the worst performer, most likely because the number of samples here is only N = 100 per input condition and neural network training is significantly impacted. It is interesting to look at the stationary responses in yeast which have not previously been analyzed in detail. First, low estimates provided by linear SVM for Msn2 and Dot6 imply that information in the stationary regime, if present, cannot be extracted by the linear classifier. Second, the Gaussian decoder also performs poorly in the stationary regime, potentially indicating that the relevant features are encoded in higher-than-pairwise order statistics of the response (e. g. , pulses could be “sparse” features as in sparse coding [58]); it is, however, hard to exclude small number of training samples as the explanation for the poor performance of the Gaussian decoder. Third, K-nearest-neighbor estimator also yields low estimates, either due to small sample number or low signal-to-noise ratio, the regime for which knn method has been observed to show reduced performance [40]. A particularly worrying feature of the knn estimates is their non-robust dependence on the length of the trajectory T. As S6 Fig shows, the performance of knn peaks at T ≈ 50 min and then drops, even well into unrealistic negative estimates for T ≈ 400 min (corresponding to the highest dimensionality d = 170 of discrete trajectories). While it is possible to make an ad hoc choice to always select trajectory duration at which the estimate peaks, the performance of kernelized SVM is, in comparison, extremely well behaved and increases monotonically with T, as theoretically expected. Finally, nonlinear SVM estimator extracts up to 0. 4 bits of information about condition per gene, more than half of the information in the early transient period. This is even though on average the response trajectories for the two conditions, 2% glucose and 0. 1% glucose, for Msn2 and Dot6 are nearly identical. For Sfp1 there is a notable difference in the mean response, which the linear estimator can use to provide a ∼ 0. 15 bits of information, yet still significantly below ∼ 0. 4 bits extracted by the nonlinear SVM. For both transient and stationary responses in yeast, our results are qualitatively in line with the expectations from the synthetic example cases—given the small number of trajectories, tightest and most robust estimates are provided by the decoding information estimator based on nonlinear (kernelized) SVM. Regardless of the decoding methodology and even without small sample corrections at N = 100 trajectories per input, our estimates are not significantly impacted by the well-known information estimation biases thanks to the dimensionality reduction that decoding provides by mapping high dimensional trajectories X back into the space for inputs U which is low dimensional; this is verified in S7 Fig by estimating the (zero) information in trajectories whose input labels have been randomly assigned. Random pulses that encode stationary environmental signals have been observed for at least 10 transcription factors in yeast [53] and for tens of transcription factors in mammalian cells [59]. Recent studies investigated the role of the pulsatile dynamics in cellular decision-making [57,60]. Nevertheless, methods for quantifying the information encoded in stochastic pulses are still in their infancy. Our nonlinear SVM decoding estimates convincingly show that there is information to be learned at the single cell level from the stationary stochastic pulsing. An interesting direction for future work is to ask whether hand-crafted features of the response trajectories (pulse frequency, amplitude, shape, etc) can extract as much information from the trajectories as the generic SVM classifier: for that, one would construct for each response trajectory a “feature vector” by hand, compute the linear SVM decoding bound information estimate from the feature vectors, and compare that to the kernelized SVM estimate over the original trajectories. This approach is a generic and operationally-defined path for finding “sufficient statistics” of the response trajectories—or a compression of the original signal to the relevant set of features—in the information-theoretic sense. A different picture emerges from the mammalian signaling network data shown in Fig 9. The key difference here is the order of magnitude larger number of sample trajectories per condition compared to yeast data. Most of the information seems linearly separable in both the early and late response periods, as evidenced by the success of the linear SVM based estimator whose performance is not improved upon by the kernelized SVM (indeed, for early ERK response period linear SVM gives a slightly higher estimate than the nonlinear version). The big winner on this dataset is the neural-network-based estimator that yields the best performance in all conditions among the decoding-based estimators, likely owning to sufficient training data. As before, the Gaussian decoder shows mixed performance which can get competitive with the best estimators under some conditions. Lastly, knn appears to do well except on the late Ca2+ data (perhaps due to low signal-to-noise ratio). It also shows counter-intuitive non-monotonic behavior with trajectory duration T in S8 Fig (cf. with Fig 2C of Ref [27], where the analysis of information conveyed in dynamical signals as a function of trajectory duration was also very revealing about signaling in yeast). Once again it is worth keeping in mind that knn is estimating the full mutual information which could be higher than the information decodable from single responses.
Increasing availability of single-cell time-resolved data should allow us to address open questions regarding the amount of information about the external world that is available in the time-varying concentrations, activation or localization patterns, and modification state of various biochemical molecules. Do full response trajectories provide more information than single temporal snapshots, as early studies suggest? Is this information gain purely due to noise averaging enabled by observing multiple snapshots, or—more interestingly—due to the ability of these intrinsically high-dimensional signals to provide a richer representation of the cellular environment? Can we isolate biologically relevant features of the response trajectories, e. g. , amplitude, frequency, pulse shape, relative phase or timing, without a priori assuming what these features are? How can cells read out the environmental state from these response trajectories and how close to the information-theoretic bounds is this readout process? More broadly, a framework for analyzing complete response trajectories in signaling or genetic regulatory networks at the single cell level could lead to architectural and functional constraints on the biological network [27], and allow us to further pursue the ideas of optimal information representation in biological systems [9]. Here, we made methodological steps towards answering these questions by focusing on two related problems: first, if we are given a full stochastic description of a biochemical reaction network, under what conditions can we theoretically compute information transmission through this network and various related bounds; second, if we are given real data with no description of the network, what are tractable schemes to estimate the information transmission. We show that when the complete state of the reaction network is observed and the inputs are discrete sets of reaction rates, there exist tractable Monte Carlo approximation schemes for the information transmission. These exact results that we compute for three simple biological network examples then serve to benchmark a family of decoding-based model-free estimators and compare their performance to the commonly-used knn estimator. We show that decoding-based estimators can closely approach the optimal decoder performance and in many cases perform better than knn, especially with typical problem dimensions (d ∼ 1 − 100) and typical number of sample trajectories (N ∼ 102 − 103). This is especially true when we ask about the combinatorial representation of the environmental state in the time trajectories of several jointly observed chemical species, as in our previous work [27], where alternative information estimation methods usually completely fail due to the high dimensionality of the input space. It is necessary to emphasize the flexibility of the decoding approach: decoding-based information estimation is based directly on the statistical problems of classification (for discrete input variable, U) or regression (for continuous input variable, U), so any classification / regression algorithm with good performance can provide the basis for information estimation. Concretely, for problems in the low data regime (small N), linear or kernelized SVM approaches appear powerful, while at larger N neural-network-based schemes can provide a better performance and thus typically a tighter information lower bound. In contrast to information approximations for which it is often impossible to assess their precision or bias (or even its sign) when the dimension, d, of the problem is large, the decoding approach yields a conservative estimate of the true information. Statistical algorithms underlying decoding-based estimations have the extra advantage that, (i), we may be able to gain biological insight by inspecting which features of the response carry the relevant stimulus information (e. g. , by looking at the linear kernels or features that neural networks extract in their various layers); (ii), pick a decoding algorithm based on features previously reported as relevant (e. g. , the Gaussian decoder for second-order statistics as in Example 3); (iii), estimate the information as a function of trajectory duration; and (iv), gain confidence in our estimates by testing their performance on withheld data. While we tested these estimators on a very restricted set of toy examples in order to be able to compare to analytically computed results, model-free decoding-based approaches are applicable more generally, e. g. , to complex, partially-observed reaction systems, or networks with significant contribution of cell-to-cell variability or extrinsic noise. By construction, decoding-based estimators only provide a lower bound to the true information. This, however, could turn out to be a smaller problem in practice than it appears in theory, especially for biochemical reaction networks. First, our extension to the Feder-Merhav bound provides us with an estimate of how large the gap between the true information and the decoded estimation can be. The bound is not tight on our examples, and can only be applied when the optimal MAP decoder can be constructed [61,62]. Second, and perhaps more importantly, information that can be decoded after single input presentations is the quantity that is likely more biologically relevant than the true channel capacity, if the organisms are under constraint to respond to the environmental changes quickly. Typically, organisms across the complexity scale operate under speed-accuracy tradeoffs [63]: faster decisions based on noisy information lead to more errors and, conversely, with enough time to integrate sensory information errors can be reduced. When speed is at a premium or relevant inputs are sparse, decisions need to be taken after single input presentations. In this case, decoding-based estimation should not be viewed as an approximate but rather as the correct methodology for the biological problem at hand. Of course, there is still the question of whether the model-free decoders that we use on real data can achieve a performance that is close to the optimal MAP decoder that represents the absolute performance limit. While there is no general way to answer this question, it appears that simple SVM decoding schemes work well when the response trajectories differ in their conditional mean, and neural networks as general approximators can be used to check for more complicated encoding features when data is plentiful. Unlike in neuroscience, there is much less clarity about what kind of read-out or decoding operations biochemical networks can mechanistically realize to mimic the functioning of our in silico decoders, and it may be challenging to biochemically implement even arbitrary linear classification of response trajectories. Until experimentally shown otherwise, it thus appears reasonable to proceed with the assumption that environmental signals can be read out from the time-dependent internal chemical state with a simple repertoire of computations. We also mention a caveat when using decoding-based estimators that rely on classification or regression methods with large expressive power, such as neural networks. While it is possible to successfully guard against overfitting within the same dataset using cross-validation, scientific insights into biological function often require generalization beyond one particular dataset. Typically, we ask for generalization at least over independent experimental replicates, but sometimes even over similar (but not same) external conditions, strains, or experimental setups. This can present a serious issue if e. g. , neural networks overfit to such systematic variations between replicates or conditions even when such variations are not biologically relevant. Regularization alone will not necessarily guard against this, unless the networks are actually trained over a subset of all data on which they will be tested. A pertinent recommendation here is to evaluate the difference in performance of expressive decoding-based estimators when trained over a subset or over all replicates, and to compare that to the generalization of less-expressive methods for which the sufficient statistics are known (e. g. , linear or Gaussian decoders). We conclude by emphasizing a simple yet important point. The decoding-based approach that we introduced here should also motivate us to look beyond methodological problems of significance and estimation, to truly biological problems of cellular decision making. Currently, data on biological regulatory processes is often analyzed by looking for “statistically significant differences” in the network response for, say, two possible network inputs. For example, one may report that the steady-state mean expression level of a certain gene is significantly larger in the stimulated vs unstimulated condition, with the statistical significance of the mean difference established through an appropriate statistical test that takes into account the number of collected population samples. While statistical significance is a necessary condition to validly report any difference in the response, it is very different from the question of whether a single cell could discriminate the two conditions given access only to its own expression levels. In caricature, population-level statistics tell us with what confidence we, as scientists having access to N samples, can discriminate between conditions given some biological readout; decoding based information estimates, on the other hand, are relevant to the N = 1 case of individual cells. We hope that further work along the latter path can clarify and quantify better the difficult constraints and conditions under which real cells need to act based on individual noisy readouts of their stochastic biochemistry.
Cells represent changes in their own state or in the state of their environment by temporally varying the concentrations of intracellular signaling molecules, mimicking in a simple chemical context the way we humans represent our thoughts and observations through temporally varying patterns of sounds that constitute speech. These time-varying concentrations are used as signals to regulate downstream molecular processes, to mount appropriate cellular responses for the environmental challenges, or to communicate with nearby cells. But how precise and unambiguous is such chemical communication, in theory and in data? On the one hand, intuition tells us that many possible environmental changes could be represented by variation in concentration patterns of multiple signaling chemicals; on the other, we know that chemical signals are inherently noisy at the molecular scale. Here we develop data analysis methodology that allows us to pose and answer these questions rigorously. Our decoding-based information estimators, which we test on simulated and real data from yeast and mammalian cells, measure how precisely individual cells can detect and report environmental changes, without making assumptions about the structure of the chemical communication and using only the amounts of data that is typically available in today’s experiments.
information entropy
neural networks
engineering and technology
electronics
applied mathematics
signaling networks
random variables
neuroscience
covariance
simulation and modeling
algorithms
mathematics
network analysis
research and analysis methods
computer and information sciences
signal decoders
probability theory
approximation methods
information theory
biology and life sciences
physical sciences
Estimating information in time-varying signals