1. Physics of Living Systems
Download icon

Theory for the optimal detection of time-varying signals in cellular sensing systems

  1. Giulia Malaguti
  2. Pieter Rein ten Wolde  Is a corresponding author
  1. AMOLF, Science Park, Netherlands
Research Article
  • Cited 0
  • Views 712
  • Annotations
Cite this article as: eLife 2021;10:e62574 doi: 10.7554/eLife.62574

Abstract

Living cells often need to measure chemical concentrations that vary in time, yet how accurately they can do so is poorly understood. Here, we present a theory that fully specifies, without any adjustable parameters, the optimal design of a canonical sensing system in terms of two elementary design principles: (1) there exists an optimal integration time, which is determined by the input statistics and the number of receptors; and (2) in the optimally designed system, the number of independent concentration measurements as set by the number of receptors and the optimal integration time equals the number of readout molecules that store these measurements and equals the work to store these measurements reliably; no resource is then in excess and hence wasted. Applying our theory to the Escherichia coli chemotaxis system indicates that its integration time is not only optimal for sensing shallow gradients but also necessary to enable navigation in these gradients.

Introduction

Living cells continually have to respond and adapt to changes in their environment. They often do so on a timescale that is comparable to that of the environmental variations. Examples are cells that during their development differentiate in response to time-varying morphogen gradients (Durrieu et al., 2018) or cells that navigate through their environment (Tostevin and ten Wolde, 2009; Sartori and Tu, 2011; Long et al., 2016). These cells shape, via their movement, the statistics of the input signal, such that the timescale of the input fluctuations becomes comparable to that of the response. In all these cases, it is important to understand how accurately the cell can estimate chemical concentrations that vary in time.

Cells measure chemical concentrations via receptors on their surface. These measurements are inevitably corrupted by the stochastic arrival of the ligand molecules by diffusion and by the stochastic binding of the ligand to the receptor. Wiener and Kolmogorov (Extrapolation, 1950; Kolmogorov, 1992) and Kalman, 1960 have developed theories for the optimal strategy to estimate signals in the presence of noise. Their filtering theories have been employed widely in engineering, and in recent years they have also been applied to cell signaling. They have been used to show that time integration can improve the sensing of time-varying signals by reducing receptor noise, although it cannot remove this input noise completely because of signal distortion (Andrews et al., 2006; Hinczewski and Thirumalai, 2014; Becker et al., 2015). It has been shown that circadian systems can adapt their response to the statistics of the input signal, as predicted by Kalman filtering theory (Husain et al., 2019). Moreover, Wiener–Kolmogorov filtering theory has been employed to derive the optimal topology of the cellular network depending on the statistics of the input signal (Becker et al., 2015). Negative feedback and incoherent feedforward, which are common motifs in cell signaling (Alon, 2007), make it possible to predict future signal values via signal extrapolation, which is useful when the past signal contains information about the future in addition to the current signal (Becker et al., 2015).

The precision of sensing depends not only on the topology of the cellular sensing network but also on the resources required to build and operate it. Receptors and time are needed to take the concentration measurements (Berg and Purcell, 1977), downstream molecules are necessary to store the ligand-binding states of the receptor in the past, and energy is required to store these states reliably (Govern and Ten Wolde, 2014a). Many studies have addressed the question how receptors and time limit the precision of sensing static concentrations that do not vary on the timescale of cellular response (Berg and Purcell, 1977; Bialek and Setayeshgar, 2005; Wang et al., 2007; Rappel and Levine, 2008; Endres and Wingreen, 2009; Hu et al., 2010; Mora and Wingreen, 2010; Govern and Ten Wolde, 2012; Mehta and Schwab, 2012; Govern and Ten Wolde, 2014a; Govern and Ten Wolde, 2014b; Kaizu et al., 2014; Ten Wolde et al., 2016; Mugler et al., 2016; Fancher and Mugler, 2017). In addition, progress has been made in understanding how the number of readout molecules and energy set the precision of sensing static signals (Mehta and Schwab, 2012; Govern and Ten Wolde, 2014a; Govern and Ten Wolde, 2014b). Yet, what the resource requirements for sensing time-varying signals are is a wide open question. In particular, it is not known how the number of receptor and readout molecules, time, and power required to maintain a desired sensing precision depend on the strength and the timescale of the input fluctuations.

In this article, we present a theory for the optimal design of cellular sensing systems as set by resource constraints and the dynamics of the input signal. The theory applies to one of the most common motifs in cell signaling, a receptor that drives a push–pull network, which consists of a cycle of protein activation and deactivation (Goldbeter and Koshland, 1981, see Figure 1). These systems are omnipresent in prokaryotic and eukaryotic cells (Alon, 2007). Examples are GTPase cycles, as in the Ras system, phosphorylation cycles, as in MAPK cascades, and two-component systems like the chemotaxis system of Escherichia coli. Push–pull networks constitute a simple exponential filter (Hinczewski and Thirumalai, 2014; Becker et al., 2015), in which the current output depends on the current and past input (with past input values contributing to the output with a weight that decays exponentially with time back into the past). Wiener–Kolmogorov filtering theory (Extrapolation, 1950; Kolmogorov, 1992) shows that these networks are optimal for estimating signals that are memoryless (Becker et al., 2015), meaning that the past input does not contain information that is not already present in the current input. These networks are useful because they act as low-pass filters, removing the high-frequency receptor–ligand-binding noise (Andrews et al., 2006; Hinczewski and Thirumalai, 2014; Becker et al., 2015). Push–pull networks thus enable the cell to employ the mechanism of time integration, in which the cell infers the concentration not from the instantaneous number of ligand-bound receptors, but rather from the average receptor occupancy over an integration time (Berg and Purcell, 1977). Our theory gives a unified description in terms of all the cellular resources, protein copies, time, and energy, that are necessary to implement this mechanism of time integration. It does not address the sensing strategy of maximum-likelihood estimation (Endres and Wingreen, 2009; Mora and Wingreen, 2010; Lang et al., 2014; Hartich and Seifert, 2016; Ten Wolde et al., 2016) or Bayesian filtering (Mora and Nemenman, 2019).

The cell signaling network.

(a) The time-varying ligand concentration is modeled as a memoryless (Markovian) signal with mean L¯, variance σL2, and correlation time τL=λ-1. A free ligand molecule L (light blue circle) can bind at rate k1 to a free receptor R (magenta protein) on the cell membrane (black line), forming the complex RL, and unbind at rate k2 from RL. The correlation time of the receptor state is τc. The complex RL catalyzes the phosphorylation reaction, driven by adenosine triphosphate (ATP) conversion, of a downstream readout from the unphosphorylated (inactive) state x to the phosphorylated (active) state x*, with rate kf. The phosphorylated readout then spontaneously decays to the x state with rate kr. Microscopic reverse reactions of each signaling pathway are represented by dashed arrows. The relaxation time of the push–pull network is τr. (b) Free-energy landscape of a readout molecule across the activation/deactivation reactions. Fuel turnover, provided by ATP conversion, drives the activation (phosphorylation) reaction characterized by the forward rate kf and its microscopic reverse rate k-f (green arrows). Associated with this activation reaction is a free-energy drop Δμ1=logkfx¯k-fx¯*. The deactivation pathway corresponds to the spontaneous release of the inorganic phosphate; it is characterized by the rate kr and its microscopic reverse k-r (blue arrows) and corresponds to a free-energy drop Δμ2=logkrx¯*k-rx¯.

While filtering theories are powerful tools for predicting the optimal topology and response dynamics of the cellular sensing network (Andrews et al., 2006; Hinczewski and Thirumalai, 2014; Becker et al., 2015), they do not naturally reveal the resource requirements for sensing. Our theory therefore employs the sampling framework of Govern and Ten Wolde, 2014a and extends it here to time-varying signals. This framework is based on the observation that the cell estimates the current ligand concentration not from the current number of active readout molecules directly, but rather via the receptor: the cell uses its push–pull network to estimate the receptor occupancy from which the ligand concentration is then inferred (see Figure 2). To elucidate the resource requirements for time integration, the push–pull network is viewed as a device that employs the mechanism of time integration by discretely sampling, rather than continuously integrating, the state of the receptor via collisions of the readout molecules with the receptor proteins (see Figure 2). During each collision, the ligand-binding state of the receptor protein is copied into the activation state of the readout molecule (Ouldridge et al., 2017). The readout molecules thus constitute samples of the receptor state, and the fraction of active readout molecules provides an estimate of the average receptor occupancy. The readout activation states have, however, a finite lifetime, which means that this is an estimate of the (running) average receptor occupancy over this lifetime, which indeed sets the receptor integration time τr. The cell can estimate the current ligand concentration L from this estimate of the average receptor occupancy pτr over the past integration time τr because there is a unique one-to-one mapping between pτr and L. This mapping pτr(L) is the dynamic input–output relation and differs from the conventional static input–output relations used to describe the sensing of static concentrations that do not vary on the timescale of the response (Berg and Purcell, 1977; Bialek and Setayeshgar, 2005; Kaizu et al., 2014; Ten Wolde et al., 2016) in that it depends not only on the response time of the system but also on the dynamics of the input signal.

The precision of estimating a time-varying ligand concentration L.

(a) The cell estimates the current ligand concentration L=L(t) by estimating the average receptor occupancy pτr over the past integration time τr and by inverting the dynamic input–output relation pτr(L) (black solid line). The latter describes the mapping between the current concentration L(t) of the time-varying signal and the average receptor occupancy pτr over the past τr, see also (b); it depends on the timescale τL of the input signal and hence differs from the conventional static input–output relation p(Ls), which describes the mapping between the average receptor occupancy and a static ligand concentration Ls that does not vary in time (gray solid line). The squared error in the estimate of the concentration (δL^)2=σp^τr2/g~Lpτr2 depends on the variance σp^τr2 in the estimate of the average receptor occupancy p^τr and the dynamic gain g~Lpτr, which is the slope of pτr(L). Only in the limit τc,τrτL, does pτr(L) reduce to (the linearized form of) p(Ls) and does the dynamic gain g~Lpτr become the static gain gLp, which is the slope of p(Ls) at the average ligand concentration L¯. The input distribution, shown in blue, has width σL. (b) The average receptor occupancy pτr over the past integration time τr is estimated via the downstream network, which is modeled as a device that discretely samples the ligand-binding state of the receptor via its readout molecules x (Govern and Ten Wolde, 2014a); the fraction of modified readout molecules provides an estimate of pτr. The sensing error has two contributions (Equation 6): sampling error and dynamical error. The sampling error arises from the error in the estimate of pτr that is due to the stochasticity of the sampling process; it depends on the number of samples, their independence, and their accuracy. (c) The dynamical error arises because the current ligand concentration L(t) is estimated via the average receptor occupancy pτr over the past integration time τr: the latter depends on the ligand concentration in the past τr, which will, in general, deviate from the current concentration. Two different input trajectories (L1 in blue, L2 in green) ending at time t at the same value L(t) (red dot) lead to different estimates of L(t) due to their different average receptor occupancy (pτr,1>pτr,2) in the past τr.

Our theory reveals that the sensing error can be decomposed into two terms, which each depend on collective variables that reveal the resource requirements for sensing. One term, the sampling error, describes the sensing error that arises from the finite accuracy by which the receptor occupancy is estimated. This error depends on the number of receptor samples, as set by the number of receptors, readout molecules, and the integration time; their independence, as given by the receptor-sampling interval and the timescale of the receptor–ligand-binding noise; and their reliability, as determined by how much the system is driven out of thermodynamic equilibrium via fuel turnover. The other term is the dynamical error and is determined by how much the concentration in the past integration time reflects the current concentration that the cell aims to estimate; it depends on the integration time and timescale of the input fluctuations.

Our theory gives a comprehensive view on the optimal design of a cellular sensing system. Firstly, it reveals that the resource allocation principle of Govern and Ten Wolde, 2014a can be generalized to time-varying signals. There exist three fundamental resource classes – receptors and their integration time, readout molecules, and power and integration time – which each fundamentally limit the accuracy of sensing; and, in an optimally designed system, each resource class is equally limiting so that none of them is in excess and thus wasted. However, in contrast to sensing static signals, time cannot be freely traded against the number of receptors and the power to achieve a desired sensing precision: there exists an optimal integration time that maximizes the sensing precision, which arises as a trade-off between the sampling error and dynamical error. Together with the resource allocation principle, it completely specifies, without any adjustable parameters, the optimal design of the system in terms of its resources protein copies, time, and energy.

Our theory also makes a number of specific predictions. The optimal integration time decreases as the number of receptors is increased because this allows for more instantaneous measurements. Moreover, the allocation principle reveals that when the input varies more rapidly both the number of receptors and the power must increase to maintain a desired sensing precision, while the number of readout molecules does not.

Finally, we apply our theory to the chemotaxis system of E. coli. This bacterium searches for food via a run-and-tumble strategy (Berg and Brown, 1972), yielding a fluctuating input signal. In small gradients, the timescale of these input fluctuations is set by the typical run time of the bacterium, which is on the order of a few seconds (Berg and Brown, 1972; Taute et al., 2015), while the strength of these fluctuations is determined by the steepness of the gradient. Interestingly, experiments have revealed that E. coli can sense extremely shallow gradients, with a length scale of approximately 104µm (Shimizu et al., 2010), raising the question how accurately E. coli can measure the concentration and whether this is sufficient to determine whether during a run it has changed, even in these shallow gradients. To measure the concentration, the chemotaxis system employs a push–pull network to filter out the high-frequency receptor–ligand-binding noise (Sartori and Tu, 2011). Applying our theory to this system predicts that the measured integration time, on the order of 100 ms (Sourjik and Berg, 2002), is not only sufficient to enable navigation in these shallow gradients but also necessary. This suggests that this system has evolved to optimally sense shallow concentration gradients.

Results

Theory: model

We consider a single cell that needs to sense a time-varying ligand concentration L(t) (see Figure 1a). The ligand concentration dynamics is modeled as a stationary memoryless, or Markovian, signal specified by the mean (total) ligand concentration L¯, the variance σL2, and the correlation time τL=λ-1, which determines the timescale on which input fluctuations decay. It obeys Gaussian statistics (Tostevin and ten Wolde, 2010).

The concentration is measured via RT receptor proteins on the cell surface, which independently bind the ligand (Ten Wolde et al., 2016), L+Rk2 k1 RL. The correlation time of the receptor state, which is the timescale on which fluctuations in the number of ligand-bound receptors regresses to the mean, is given by τc=1/(k1L¯+k2) (Berg and Purcell, 1977; Bialek and Setayeshgar, 2005; Kaizu et al., 2014; Ten Wolde et al., 2016). It determines the timescale on which independent concentration measurements can be made.

The ligand-binding state of the receptor is read out via a push–pull network (Goldbeter and Koshland, 1981). The most common scheme is phosphorylation fueled by the hydrolysis of adenosine triphosphate (ATP) (see Figure 1b). The receptor, or an enzyme associated with it such as CheA in E. coli, catalyzes the modification of the readout, x+RL+ATP  x+RL+ADP. The active readout proteins x* can decay spontaneously or be deactivated by an enzyme, such as CheZ in E. coli, xx+Pi. Inside the living cell the system is maintained in a non-equilibrium steady state by keeping the concentrations of ATP, adenosine diphosphate (ADP), and inorganic phosphate (Pi) constant. We absorb their concentrations and the activities of the kinase and, if applicable, phosphatase in the (de)phosphorylation rates, coarse-graining the (de)modification reactions into instantaneous second-order reactions: x+RL kfkf x+RL, x krkr x. This system has a relaxation time τr=1/[(kf+k-f)RL¯+kr+k-r] (Govern and Ten Wolde, 2014a), which describes how fast fluctuations in x* relax. It determines how long x* can carry information on the ligand-binding state of the receptor; τr thus sets the integration time of the receptor state.

Theory: inferring concentration from receptor occupancy

The central idea of our theory is illustrated in Figure 2a: the cell employs the push–pull network to estimate the average receptor occupancy pτr over the past integration time τr. It then uses this estimate p^τr to infer the current concentration L via the dynamic input–output relation pτr(L), which provides a one-to-one mapping between pτr and L.

Dynamic input–output relation

The mapping pτr(L) is the dynamic input–output relation. It gives the average receptor occupancy over the past integration time τrgiven that the current value of the input signal is L=L(t) (see Figure 2a). Here, the average is not only over the noise in receptorligand binding and readout activation (Figure 2b) but also over the subensemble of past input trajectories that each end at the same current concentration L (Figure 2c; Tostevin and ten Wolde, 2010; Hilfinger and Paulsson, 2011; Bowsher et al., 2013). In contrast to the conventional static inputoutput relation p(Ls), which gives the average receptor occupancy p for a steady-state ligand concentration Ls that does not vary in time, the dynamic inputoutput relation takes into account the dynamics of the input and the finite response time of the system. It depends on all timescales in the problem: the timescale of the input, τL, the receptorligand correlation time τc, and the integration time τr. Only when τLτc,τr does the dynamic inputoutput pτr(L) become equal to the static inputoutput relation p(Ls).

Sensing error

Linearizing the dynamic input–output relation pτr(L) around the mean ligand concentration L¯ (see Figure 2a) and using the rules of error propagation, the expected error in the concentration estimate is

(1) (δL^)2=σp^τr2g~Lpτr2.

Here, σp^τr2 is the variance in the estimate p^τr of the average receptor occupancy over the past τrgiven that the current input signal is L (see Figure 2a). The quantity g~Lpτr is the dynamic gain, which is the slope of the dynamic input–output relation pτr(L); it determines how much an error in the estimate of pτr propagates to that in L. Equation 1 generalizes the expression for the error in sensing static concentrations (Berg and Purcell, 1977; Bialek and Setayeshgar, 2005; Wang et al., 2007; Mehta and Schwab, 2012; Kaizu et al., 2014; Govern and Ten Wolde, 2014a; Ten Wolde et al., 2016) to that of time-varying concentrations.

Signal-to-noise ratio

Together with the distribution of input states, the sensing error (δL^)2 determines how many distinct signal values the cell can resolve. The latter is quantified by the signal-to-noise ratio (SNR), which is defined as

(2) SNRσL2(δL^)2.

Here, σL2 is the variance of the ligand concentration L(t); because the system is stationary and time invariant, we can omit the argument in L(t) and write L=L(t). The variance σL2 is a measure for the total number of input states, such that the SNR gives the number of distinct ligand concentrations the cell can measure. Using Equation 1, it is given by

(3) SNR=g~Lpτr2σp^τr2σL2.

The SNR also yields the mutual information I(x*;L)=1/2ln(1+SNR) between the input L and output x* (Tostevin and ten Wolde, 2010).

Readout system samples receptor state

Receptor time averaging is typically conceived as a scheme in which the receptor state is averaged via the mathematical operation of an integral: pτr=1/τr0τrp(t)𝑑t. Yet, readout proteins are discrete components that interact with the receptor in a discrete and stochastic fashion. To derive the dynamic gain g~Lpτr and error in estimating pτr, σp^τr2 (Equation 3), we therefore view the push–pull network as a device that discretely samples the receptor state (see Figure 2b; Govern and Ten Wolde, 2014a). The principle is that cells employ the activation reaction x+RLx*+RL to store the state of the receptor in stable chemical modification states of the readout molecules. Readout molecules that collide with a ligand-bound receptor are modified, while those that collide with an unbound receptor are not (Figure 2b). The readout molecules serve as samples of the receptor at the time they were created, and collectively they encode the history of the receptor: the fraction of samples that correspond to ligand-bound receptors is the cell’s estimate for pτr. Indeed, this is the discrete and stochastic implementation of the mechanism of time integration. The effective number of independent samples depends not only on the creation of samples, x+RLx*+RL, but also on their decay and accuracy. Samples decay via the deactivation reaction x*x, which means that they only provide information on the receptor occupancy over the past τr. In addition, both the activation and the deactivation reaction can happen in their microscopic reverse direction, which corrupts the coding, that is, the mapping between the ligand-binding states of the receptor proteins and the activation states of the readout molecules. Energy is needed to break time reversibility and protect the coding. Furthermore, for time-varying signals, we also need to recognize that the samples correspond to the ligand concentration over the past integration time τr, which will in general differ from the current concentration L that the cell aims to estimate (see Figure 2c). While a finite τr is necessary for time integration, it will, as we show below, also lead to a systematic error in the estimate of the concentration that the cell cannot reduce by taking more receptor samples.

This analysis reveals that the dynamic gain is (see Appendix 1)

(4) g~Lpτr=gLp(1+τcτL)-1(1+τrτL)-1.

Only when τLτr,τc is the average ligand concentration over the ensemble of trajectories ending at δL(t) equal to the current concentration δL(t) (Figure 2c) and does g~Lpτr become equal to its maximal value, the static gain gLp=p(1-p)/L¯, where p is the average receptor occupancy averaged over all values of δL(t). The analysis also reveals that the error in pτr can be written as (see Appendix 1, Equation 29)

(5) σp^τr2=σp^τr2,samp+σp^τr2,dyn,

where σp^τr2,samp is a statistical error due to the stochastic sampling of the receptor and σp^τr2,dyn is a systematic error arising from the dynamics of the input, as elucidated in Figure 2b, c.

Central result

To know how the error σp^τr2 in the estimate of pτr propagates to the error (δL^)2 in the estimate of the current ligand concentration, we divide σp^τr2 by the dynamic gain g~Lpτr given by Equation 4 (see Equation 1). For the full system, the reversible push–pull network, this yields via Equation 3 the central result of our article, the SNR in terms of the total number of receptor samples, their independence, their accuracy, and the timescale on which they are generated:

(6) SNR1=(1+τcτL)2(1+τrτL)2[(L¯/σL)2p(1p)N¯I+(L¯/σL)2(1p)2N¯eff]sampling error+(1+τcτL)(1+τrτL)(1+τcτrτL(τc+τr))1dynamical error.

This expression shows that the sensing error SNR-1 can be decomposed into two distinct contributions, which each have a clear interpretation: the sampling error, arising from the stochasticity in the sampling of the receptor state, and the dynamical error, arising from the dynamics of the input.

When the timescale of the ligand fluctuations τL is much longer than the receptor correlation time τc and the integration time τr, τLτr,τc, the dynamical error reduces to zero and only the sampling error remains. Here, N¯eff is the total number of effective samples and N¯I is the number of these that are independent (Govern and Ten Wolde, 2014a). For the full system, they are given by

(7) N¯I=1(1+2τc/Δ)fI(eβΔμ11)(eβΔμ21)eβΔμ1qn˙τrpN¯N¯eff.

The quantity n˙=kfpRTx¯-k-fpRTx¯* is the net flux of x around the cycle of activation and deactivation, with RT the total number of receptor proteins and x¯ and x¯* the average number of inactive and active readout molecules, respectively. It equals the rate at which x is modified by the ligand-bound receptor; the quantity n˙/p is thus the sampling rate of the receptor, be it ligand bound or not. Multiplied with the relaxation rate τr, it yields the total number of receptor samples N¯ obtained during τr. However, not all these samples are reliable. The effective number of samples is N¯eff=qN¯, where 0<q<1 quantifies the quality of the sample. Here, β=1/(kBT) is the inverse temperature, Δμ1 and Δμ2 are the free-energy drops over the activation and deactivation reaction, respectively, with Δμ=Δμ1+Δμ2 the total drop, determined by the fuel turnover (see Figure 1b). If the system is in thermodynamic equilibrium, Δμ1=Δμ2=Δμ=0, q0 and the system cannot sense because the ligand-binding state of the receptor is equally likely to be copied into the correct modification state of the readout as into the incorrect one. In contrast, if the system is strongly driven out of equilibrium and Δμ1,Δμ2, then, during each receptor–readout interaction, the receptor state is always copied into the correct activation state of the readout; the sample quality parameter q thus approaches unity and N¯effN¯. Yet, even when all samples are reliable, they may contain redundant information on the receptor state. The factor fI is the fraction of the N¯eff samples that are independent. It reaches unity when the receptor sampling interval Δ=2τr/(N¯eff/RT) becomes larger than the receptor correlation time τc.

When the number of samples becomes very large, N¯I,N¯eff, the sampling error reduces to zero. However, the sensing error still contains a second contribution, which, following Bowsher et al., 2013, we call the dynamical error. This contribution only depends on timescales. It arises from the fact that the samples encode the receptor history and hence the ligand concentration over the past τr, which will, in general, deviate from the quantity that the cell aims to predict – the current concentration L. This contribution yields a systematic error, which cannot be eliminated by increasing the number of receptor samples, their independence, or their accuracy. It can only be reduced to zero by making the integration time τr much smaller than the ligand timescale τL (assuming τc is typically much smaller than τr,τL). Only in this regime will the ligand concentration in the past τr be similar to the current concentration and can the latter be reliably inferred from the receptor occupancy, provided the latter has been estimated accurately by taking enough samples.

Importantly, the dynamics of the input signal not only affects the sensing precision via the dynamical error but also via the sampling error. This effect is contained in the prefactor of the sampling error, (1+τc/τL)2(1+τr/τL)2, which has its origin in the dynamic gain g~Lpτr (Equation 4). It determines how the sampling error σp^τr2,samp in the estimate of pτr propagates to the error in the estimate of L (see Equation 3). Only when τc,τrτL can the readout system closely track the input signal and does g~Lpτr reach its maximal value, the static gain gLp, thus minimizing the error propagation from pτr to L.

Fundamental resources

We can use Equation 6 to identify the fundamental resources for cell sensing (Govern and Ten Wolde, 2014a) and derive Pareto fronts that quantify the trade-offs between the maximal sensing precision and these resources. A fundamental resource is a (collective) variable Qi that, when fixed to a constant, puts a non-zero lower bound on SNR-1, no matter how the other variables are varied. It is thus mathematically defined as MINQi=const(SNR-1)=f(const)>0. To find these collective variables, we numerically or analytically minimized SNR-1, constraining (combinations of) variables yet optimizing over the other variables. This reveals that the SNR is bounded by (see Appendix 2)

(8) SNR-1(1+τrτL)24(L¯/σL)2h+τrτL,

where

(9) hMIN(RTτr/τc,XT,βw˙τr).

Equations 8 and 9 show that the fundamental resources are the number of receptors RT, the integration time τr, the number of readouts XT, and the power w˙=n˙Δμ.

Figure 3a, b illustrates that RT,τr,XT,w˙ are indeed fundamental: the sensing precision is bounded by the limiting resource and cannot be enhanced by increasing another resource. Panel (a) shows that when XT is small, the maximum mutual information Imax cannot be increased by raising RT: no matter how many receptors the system has, the sensing precision is limited by the pool of readout molecules and only increasing this pool can raise Imax. Yet, when XT is large, Imax becomes independent of XT. In this regime, the number of receptors RT limits the number of independent concentration measurements and only increasing RT can raise Imax. Similarly, panel (b) shows that when the power w˙ is limiting, Imax cannot be increased by RT but only by increasing w˙. Clearly, the resources receptors, readout molecules, and energy cannot compensate each other: the sensing precision is bounded by the limiting resource.

Receptors RT, readout molecules XT, and w˙ fundamentally limit sensing, and there exists an optimal integration time τr that depends on which of the resources is limiting.

(ab) RT, XT, and w˙ are fundamental resources, with no trade-offs between them. Plotted is the maximum mutual information Imax=1/2ln(1+SNRmax), obtained by minimizing Equation 6 over p and τr, for different combinations of (a) XT and RT in the irreversible limit q1 and (b) w˙ and RT for two different values of Δμ. The sensing precision is bounded by the limiting resource, RT (solid gray lines, Equation 8 with h=RT/τr/τc), XT (dashed gray line, Equation 8 with h=XT, panel a), or w˙ (dashed gray lines, Equation 8 with h=βw˙τr or h=w˙τr/(Δμ/4), panel b). (c) Imax as a function of τr for different values of RT in the Berg–Purcell limit (q1 and XT). There exists an optimal integration time τropt that maximizes the sensing precision; τropt decreases with RT. (d) In this limit, τropt depends non-monotonically on the receptor–ligand correlation time τc: it first increases with τc to sustain time-averaging, but then drops when τropt/τc becomes of order unity and time-averaging is no longer effective (see inset). (e) τropt as a function of XT for different values of RT. When XT<RT, time averaging is not possible and the optimal system is an instantaneous responder, τropt0; when XTRT, the system reaches the Berg–Purcell regime in which Imax is limited by RT rather than XT (see panel a). (f) τropt and XT as a function of w˙. When the power w˙XT/τr is limiting, the sampling error dominates and τropt equals τL to maximize XT, minimizing the sampling error; τropt then decreases to trade part of the decrease in the sampling error for a reduction in the dynamical error such that both decrease; when the sampling interval ΔτrRT/XT becomes comparable to τc, in the region marked by the yellow bar, the sampling error is no longer limited by XT, such that τr now limits both sources of error; the two sources can therefore no longer be decreased simultaneously by increasing w˙XT/τr; the system has entered the Berg–Purcell regime, where τropt is determined by RT rather than w˙ (see panel b). Parameter values unless specified: τc/τL=102; σL/L¯T=102.

Importantly, while for sensing static concentrations the products RTτr/τc and w˙τr are fundamental (Govern and Ten Wolde, 2014a), for time-varying signals RT, w˙, and τr separately limit sensing. Consequently, neither receptors RT nor power w˙ can be traded freely against time τr to reach a desired precision, as is possible for static signals. In line with the predictions of signal filtering theories (Extrapolation, 1950; Kolmogorov, 1992; Kalman, 1960), there exists an optimal integration time τropt that maximizes the sensing precision (Andrews et al., 2006; Hinczewski and Thirumalai, 2014; Becker et al., 2015; Monti et al., 2018b; Mora and Nemenman, 2019). Interestingly, its value depends on which of the resources RT, XT, and w˙ is limiting (Figure 3c–f). We now discuss these three regimes in turn.

Receptors

Berg and Purcell, 1977 pointed out that cells can reduce the sensing error by either increasing the number of receptors or taking more measurements per receptor via the mechanism of time integration. However, Equation 8 reveals that for sensing time-varying signals time integration can never eliminate the sensing error completely, as predicted also by filtering theories (Extrapolation, 1950; Kolmogorov, 1992; Kalman, 1960). Equation 8 shows that in the Berg–Purcell regime, where receptors and their integration time are limiting and h=RTτr/τc, the sensing precision does not depend on RTτr/τc, as for static signals (Govern and Ten Wolde, 2014a), but on RT and τr separately, such that an optimal integration time τropt emerges that maximizes the sensing precision (see Figure 3c). Increasing τr improves the mechanism of time integration by increasing the number of independent samples per receptor, τr/τc, thus reducing the sampling error (Equation 6). However, increasing τr raises the dynamical error. Moreover, it lowers the dynamical gain g~Lpτr, which increases the propagation of the error in the estimate of the receptor occupancy to that of the ligand concentration. The optimal integration time τropt arises as a trade-off between these three factors.

Figure 3c also shows that the optimal integration time τropt decreases with the number of receptors RT. The total number of independent concentration measurements is the number of independent measurements per receptor, τr/τc, times the number RT of receptors, N¯I=RTτr/τc. As RT increases, less measurements τr/τc per receptor have to be taken to remove the receptor–ligand-binding noise, explaining why τropt decreases as RT increases – time integration becomes less important.

Interestingly, τropt depends non-monotonically on the receptor–ligand correlation time τc (Figure 3d). When τc increases at fixed τr, the receptor samples become more correlated. To keep the mechanism of time integration effective, τr must increase with τc. However, to avoid too strong signal distortion the cell compromises on time integration by decreasing the ratio τr/τc (see inset). When τr becomes too large, the benefit of time integration no longer pays off the cost of signal distortion. Now not only the ratio τr/τc decreases but also τr itself. The sensing system switches to a different strategy: it no longer employs time integration but becomes an instantaneous sensor.

Readout molecules

To implement time integration, the cell needs to store the receptor states in the readout molecules. When the number of readout molecules XT is limiting, the sensing precision is given by Equation 8 with h=XT. This bound is saturated when τr0. This is in marked contrast to the non-zero optimal integration τropt in the Berg–Purcell regime (see Figure 3c).

To elucidate the non-trivial behavior of τropt, Figure 3e shows τropt as a function of XT. When XT is smaller than RT, the average number of samples per receptor is less than unity. In this regime, the system cannot time integrate the receptor, and to minimize signal distortion τropt0. Yet, when XT is increased, the likelihood that two or more readout molecules provide a sample of the same receptor molecule rises, and time averaging becomes possible. Yet to obtain receptor samples that are independent, the integration time τr must be increased to make the sampling interval ΔτrRT/XT larger than the receptor correlation time τc. As XT and hence the total number of samples N¯ are increased further, the number of samples that are independent, N¯I, only continues to rise when τr increases with XT further. However, while this reduces the sampling error, it also increases the dynamical error. When the decrease in the sampling error no longer outweighs the increase in the dynamical error, τropt and the mutual information no longer change with XT (see Figure 3a). The system has entered the Berg–Purcell regime in which τropt and the mutual information are given by the optimization of Equation 8 with h=RTτr/τc (gray dashed line). In this regime, increasing XT merely adds redundant samples: the number of independent samples remains N¯I=RTτropt/τc.

Power

Time integration relies on copying the ligand-binding state of the receptor into the chemical modification states of the readout molecules (Mehta and Schwab, 2012; Govern and Ten Wolde, 2014a). This copy process correlates the state of the receptor with that of the readout, which requires work input (Ouldridge et al., 2017).

The free-energy Δμ provided by the fuel turnover drives the readout around the cycle of modification and demodification (Figure 1). The rate at which the fuel molecules do work is the power w˙=n˙Δμ, and the total work performed during the integration time τr is ww˙τr. This work is spent on taking samples of receptor molecules that are bound to ligand because only they can modify the readout. The total number of effective samples of ligand-bound receptors during τr is pN¯eff (Equation 7), which means that the work per effective sample of a ligand-bound receptor is w/(pN¯eff)=Δμ/q (Govern and Ten Wolde, 2014a).

To understand how energy limits the sensing precision, we can distinguish between two limiting regimes (Govern and Ten Wolde, 2014a). When Δμ>4kBT, the quality parameter q1 (Equation 7) and the work per sample of a ligand-bound receptor is w/(pN¯eff)=Δμ (Govern and Ten Wolde, 2014a). In this irreversible regime, the SNR bound is given by Equation 8 with h=w˙τr/(Δμ/4). The power limits the sensing accuracy not because it limits the reliability of each sample but because it limits the rate n˙=w˙/Δμ at which the receptor is sampled.

When Δμ<4kBT, the system enters the quasi-equilibrium regime in which the quality parameter qβΔμ/4 (see Equation 7, noting that in the optimal system Δμ1=Δμ2=Δμ/2). The sensing bound is now given by Equation 8 with h=βw˙τr, which is larger than h=w˙τr/(Δμ/4) in the irreversible regime (where Δμ>4kBT). The quasi-equilibrium regime minimizes the sensing error for a given power constraint (Figure 3b) because this regime maximizes the number of effective measurements per work input pN¯eff/w=q/Δμ=β/4 (Govern and Ten Wolde, 2014a).

While the sensing precision for a given power and time constraint is higher in the quasi-reversible regime, more readout molecules are required to store the concentration measurements in this regime. Noting that the flux n˙=f(1-f)XTq/τr=w˙/Δμ, it follows that in the irreversible regime (q1) the number of readout molecules consuming energy at a rate w˙ is

(10) XTirr=w˙τrΔμf(1-f),

while in the quasi-equilibrium regime (qΔμ/4) it is

(11) XTqeq=w˙τr4kBTΔμ2f(1-f).

Since in the quasi-equilibrium regime Δμ<4kBT, XTqeq>XTirr.

Equation 8 shows that the sensing precision is fundamentally bounded not by the work w=w˙τr, as observed for static signals (Govern and Ten Wolde, 2014a), but rather by the power w˙ and the integration time τr separately such that an optimal integration time τropt emerges. Figure 3f shows how τropt depends on w˙. Since the system cannot sense without any readout molecules, in the low-power regime the system maximizes XT subject to the power constraint w˙XT/τr (see Equations 10 and 11) by making τr as large as possible, which is the signal correlation time τL – increasing τropt further would average out the signal itself. As w˙ is increased, XT rises and the sampling error decreases. When the sampling error becomes comparable to the dynamical error (Equation 6), the system starts to trade a further reduction in the sampling error for a reduction in the dynamical error by decreasing τropt. The sampling error and dynamical error are now reduced simultaneously by increasing XT and decreasing τropt. This continues until the sampling interval ΔRTτr/XT becomes comparable to the receptor correlation time τc, as marked by the yellow bar. Beyond this point, Δ<τc and the sampling error is no longer limited by XT but rather by τr since τr bounds the number of independent samples per receptor, τr/τc. The system has entered the Berg–Purcell regime, where τropt is determined by the trade-off between the dynamical error and the sampling error as set by the maximum number of independent samples, RTτr/τc (Figure 3c).

Optimal design

In sensing time-varying signals, a trade-off between time averaging and signal tracking is inevitable. Moreover, the optimal integration time depends on which resource is limiting, being zero when XT is limiting and finite when RT or w˙ is limiting (Figure 3). It is therefore not obvious whether these sensing systems still obey the optimal resource allocation principle as observed for systems sensing static concentrations (Govern and Ten Wolde, 2014a).

However, Equation 8 shows that when for a given integration time τr, RTτr/τc=XT=βw˙τr, the bounds on the sensing precision as set by, respectively, the number of receptors RT, the number of readout molecules XT, and the power w˙ are equal. Each of these resources is then equally limiting sensing and no resource is in excess. We thus recover the optimal resource allocation principle:

(12) RTτr/τcXTβw˙τr.

Irrespective of whether the concentration fluctuates in time, the number of independent concentration measurements at the receptor level is RTτr/τc, which in an optimally designed system also equals the number of readout molecules XT and the energy βw˙τr that are both necessary and sufficient to store these measurements reliably.

The design principle XTβw˙τr (Equation 12) predicts that there exists a driving force Δμopt that optimizes the trade-off between the number of samples and their accuracy. Noting that βw˙τr=βn˙Δμτr=βqf(1-f)XTΔμ reveals that the principle XTβw˙τr (Equation 12) specifies Δμ for the optimal system in which f1/2 and Δμ1=Δμ2=Δμ/2 via the equation q(Δμopt)=4kBT/Δμopt, where q(Δμ) is defined in Equation 7. A numerical inspection shows that to a good approximation the solution of this equation is precisely given by the crossover from the quasi-equilibrium regime to the irreversible one: Δμopt4kBT. This can be understood by noting that in the quasi-equilibrium regime XT can, for a given power and time constraint, be reduced by increasing Δμ (Equation 11) without compromising the sensing precision (Equation 8 with h=w˙τr); in this regime, increasing Δμ increases the reliability of each sample, and a smaller number of more reliable samples precisely compensates for a larger number of less reliable ones. Yet, when Δμ becomes larger than 4kBT, the system enters the irreversible regime. Here, XT corresponding to a given w˙ and τr constraint still decreases with Δμ (Equation 10), but the sensing error now increases (Equation 8 with h=w˙τr/(Δμ/4)) because each sample has become (essentially) perfect in this regime – hence, the samples’ accuracy cannot (sufficiently) increase further to compensate for the reduction in the sampling rate n˙XT/τr.

Equation 12 holds for any integration time τr, yet it does not specify τr. The cell membrane is highly crowded, and many systems employ time integration (Berg and Purcell, 1977; Bialek and Setayeshgar, 2005; Govern and Ten Wolde, 2014a). This suggests that these systems employ time integration and accept the signal distortion that comes with it simply because there is not enough space on the membrane to increase RT. Our theory then allows us to predict the optimal integration time τropt based on the premise that RT is limiting. As Equation 8 reveals, in this limit τropt does not only depend on RT but also on τc, τL, and σL/L¯:τropt=τropt(RT,τr,τL,σL/L¯). The optimal design of the system is then given by Equation 12 but with τr given by τropt=τropt(RT,τc,τL,σL/L¯):

(13) RTτropt/τcXToptβw˙optτropt.

This design principle maximizes for a given number of receptors RT the sensing precision and minimizes the number of readout molecules XT and power w˙ needed to reach that precision.

Comparison with experiment

To test our theory, we turn to the chemotaxis system of E. coli. This system contains a receptor that forms a complex with the kinase CheA. This complex, which is coarse-grained into R (Govern and Ten Wolde, 2014a), can bind the ligand L and activate the intracellular messenger protein CheY (x) by phosphorylating it. Deactivation of CheY is catalyzed by CheZ, the effect of which is coarse-grained into the deactivation rate. This push–pull network allows E. coli to measure the current concentration, and the relaxation time of this network sets the integration time for the receptor (Sartori and Tu, 2011). The system also exhibits adaptation on longer timescales due to receptor methylation and demethylation. The push–pull network and the adaptation system together allow the cell to measure concentration gradients via a temporal derivative, taking the difference between the current concentration and the past concentration as set by the adaptation time (Segall et al., 1986). A lower bound for the error in the estimate of this difference is given by the error in the estimate of the current concentration, the central quantity of our theory. Here, we ask how accurately E. coli can estimate the latter and whether the sensing precision is sufficient to determine whether during a run the concentration has changed.

Our theory predicts that if the number of receptors is limiting then the optimal integration time τropt(RT,τc,τL,σL/L¯) is given by minimizing Equation 8 with h=RTτr/τc. The number of receptor–CheA complexes depends on the growth rate and varies between RT103 and RT104 (Li and Hazelbauer, 2004). The receptor correlation time for the binding of aspartate to the Tar receptor can be estimated from the measured dissociation constant (Vaknin and Berg, 2007) and the association rate (Danielson et al., 1994), τc10ms (Govern and Ten Wolde, 2014a). The timescale τL of the input fluctuations is set by the typical run time, which is on the order of a few seconds, τL1s (Berg and Brown, 1972; Taute et al., 2015).

This leaves one parameter to be determined, (σL/L¯)2. This is set by the spatial ligand–concentration profile and the typical length of a run. We have a good estimate of the latter. In shallow gradients, it is on the order of l50μm (Berg and Brown, 1972; Taute et al., 2015; Jiang et al., 2010; Flores et al., 2012); specifically, Figure 4 of Taute et al., 2015 shows that the typical run times are 1–2 s while the typical run speeds are 20-60μms-1, yielding a run length on the order of indeed 50 µm. We do not know the spatial concentration profiles that E. coli has experienced during its evolution. We can however get a sense of the scale by considering an exponential ligand–concentration gradient. For a profile L¯(x)=L0ex/x0 with length scale x0, the relative change in the signal over the length of a run is σL/L¯(dL¯/dx)l/L¯=l/x0. We consider the range σL/L¯l/x0<1, where σL/L¯<0.1 corresponds to shallow gradients with x0500μm in which cells move with a constant drift velocity (Shimizu et al., 2010; Flores et al., 2012).

Figure 4a shows that as the gradient becomes steeper and σL/L¯l/x0 increases the optimal integration time τropt decreases. This can be understood by noting that the relative importance of the dynamical error compared to the sampling error scales with (σL/L¯)2 (Equation 6). Shallow ingredients thus allow for a larger integration time while steep gradients necessitate a shorter one.

The optimal integration time for the chemotaxis system of E. coli.

(a) The optimal integration time τropt, obtained by numerically optimizing Equation 8 with h=RTτr/τc, as a function of the relative strength of the input noise, σL/L¯, for two different copy numbers RT of the receptor–CheA complexes; for an exponential gradient with length scale x0, the relative noise strength σL/L¯l/x0, where l50μm is the run length of E. coli. It is seen that τropt increases as σL/L¯ decreases because the relative importance of the sampling error compared to the dynamical error increases. The figure also shows that τropt decreases as RT is increased because that allows for more instantaneous measurements (see also Figure 3). The red bar indicates the range of the estimated integration time of E. coli, 50ms<τr<500ms, based on its attractant and repellent response, respectively (Sourjik and Berg, 2002), divided by the input timescale τL1s based on its typical run time of about a second (Berg and Brown, 1972; Taute et al., 2015). The panel indicates that E. coli has been optimized to detect shallow concentration gradients. (b) The signal-to-noise ratio SNRτL=(σL/δL^)2τL/τr, with (σL/δL^)2=SNR given by Equation 6, as a function of σL/L¯l/x0. To be able to detect the gradient, the SNRτL must exceed unity. The panel shows that the shallowest gradient that E. coli can detect (marked with dashed red line) has, for RT=104, a length scale of x025000μm (corresponding to σL/L¯2×10-3), which is consistent with experiments based on ramp responses (Shimizu et al., 2010). Other parameter: receptor–ligand-binding correlation time τc=10ms (Vaknin and Berg, 2007; Danielson et al., 1994).

Experiments indicate that the relaxation rate of CheY is τr-12s-1 for the attractant response and 20s-1 for the repellent response (Sourjik and Berg, 2002), such that the integration time τr50-500ms (Sourjik and Berg, 2002; Govern and Ten Wolde, 2014a). Figure 4a shows that this integration time is optimal for detecting shallow gradients. Our theory thus predicts that the E. coli chemotaxis system has been optimized for sensing shallow gradients.

To navigate, the cells must be able to resolve the signal change over a run. During a run of duration τL, the system performs τL/τr independent concentration measurements. The effective error for these measurements is the instantaneous sensing error (δL^)2 divided by the number of independent measurements τL/τr:(δL^)2/(τL/τr). Hence, the SNR for these concentration measurements is SNRτL(σL/δL^)2τL/τr.

Figure 4b shows that our theory predicts that when RT=103, the shallowest gradient that cells can resolve, defined by SNRτL=1, is l/x0σL/L¯1×10-2, corresponding to x07500μm, while when RT=104, l/x02×10-3 and x025000μm. The shallowest gradient is thus on the order of x0104μm. Shimizu et al., 2010 show that E. coli cells are indeed able to sense such very shallow gradients: Figure 2A of Shimizu et al., 2010 shows that E. coli cells can detect exponential up ramps with rate r=0.001/s; using r=vr/x0, where vr10μm/s is the run speed (Jiang et al., 2010), this corresponds to x0104μm. Importantly, the predictions of our theory (Figure 4) concern the shallowest gradient that the system with the optimal integration time can resolve. These observations indicate that the optimal integration time is not only sufficient to make navigation in these very shallow gradients possible but also necessary.

Figure 4 also shows that τropt decreases as the number of receptor–CheA complex, RT, increases because the latter allows for more instantaneous measurements, reducing the need for time integration (Figure 3c). Interestingly, the data of Li and Hazelbauer, 2004 shows that the copy numbers of the chemotaxis proteins vary with the growth rate. Clearly, it would be of interest to directly measure the response time in different strains under different growth conditions.

Discussion

Here, we have integrated ideas from Tostevin and ten Wolde, 2010; Hilfinger and Paulsson, 2011; and Bowsher et al., 2013 on information transmission via time-varying signals with the sampling framework of Govern and Ten Wolde, 2014a to develop a unified theory of cellular sensing. The theory is founded on the concept of the dynamic input–output relation pτr(L). It allows us to develop the idea that the cell employs the readout system to estimate the average receptor occupancy pτr over the past integration time τr and then exploits the mapping pτr(L) to estimate the current ligand concentration L from pτr. The theory reveals that the error in the estimate of L depends on how accurately the cell samples the receptor state to estimate pτr, and on how much pτr, which is determined by the concentration in the past τr, reflects the current ligand concentration. These two distinct sources of error give rise to the sampling error and dynamical error in Equation 6, respectively.

While the system contains no less than 11 parameters, Equation 6 provides an intuitive expression for the sensing error in terms of collective variables that have a clear interpretation. The dynamical error depends only on the timescales in the problem, most notably τr/τL. The sampling error depends on how accurately the readout system estimates pτr, which is determined by the number of receptor samples, their independence, and their accuracy; yet it also depends on τr/τL via the dynamic gain, which determines how the error in the estimate of pτr propagates to that of L. The trade-off between the sampling error and dynamical error yields an optimal integration time.

Our study reveals that the optimal integration time τropt depends in a non-trivial manner on the design of the system. When the number of readout molecules XT is smaller than the number of receptors RT, time integration is not possible and the optimal system is an instantaneous responder with τropt0. When the power w˙XT/τr, rather than XT, is limiting, τropt is determined by the trade-off between the sampling error and dynamical error. In both scenarios, however, one resource, XT or w˙, is limiting the sensing precision. In an optimally designed system, all resources are equally limiting so that no resource is wasted. This yields the resource allocation principle (Equation 12), first identified in Govern and Ten Wolde, 2014a, for sensing static concentrations. The reason it can be generalized to time-varying signals is that the principle concerns the optimal design of the readout system for estimating the receptor occupancy over a given integration time τr, which holds for any type of input: the number of independent concentration measurements at the receptor level is RTτr/τc, irrespective of how the input varies, and in an optimally designed system this also equals the number of readout molecules XT and energy βw˙τr to store these measurements reliably. We thus expect that the design principle also holds for systems that sense signals that vary more strongly in time (Mora and Nemenman, 2019).

While the allocation principle Equation 12 holds for any τr, it does not specify the optimal integration time τropt. However, our theory predicts that if the number of receptors RT is limiting, then there exists a τropt that maximizes the sensing precision for that RT (Equation 8 with h=RTτr/τc). Via the allocation principle Equation 13, RT and τropt then together determine the minimal number of readout molecules XT and power w˙ to reach that precision. The resource allocation principle, together with the optimal integration time, thus completely specifies the optimal design of the sensing system.

Applying our theory to the E. coli chemotaxis system shows that this system not only obeys the resource allocation principle (Govern and Ten Wolde, 2014a) but also that the predicted optimal integration time to measure shallow gradients is in agreement with that measured experimentally (Figure 4a). This is remarkable because there is not a single fit parameter in our theory. Moreover, Figure 4b shows that the optimal integration time is not only sufficient to enable the sensing of these shallow gradients but also necessary. This is interesting because the sensing precision could also be increased by increasing the number of receptors, readout molecules, and energy devoted to sensing – but this would be costly. Our results thus demonstrate not only that the chemotaxis system obeys the design principles as revealed by our theory but also that there is a strong selection pressure to design sensing systems optimally, that is, to maximize the sensing precision given the resource constraints.

Our theory is based on a Gaussian model and describes the optimal sensing system that minimizes the mean square error in the estimate of the ligand concentration (see Equation 1). The latter is precisely the performance criterion of Wiener–Kolmogorov (Extrapolation, 1950; Kolmogorov, 1992) and Kalman, 1960 filtering theory, which, moreover, become exact for systems that obey Gaussian statistics. In fact, since our system (including the input signal) is stationary, they predict the same optimal filter, which is an exponential filter for signals that are memoryless. The signals studied here belong to this class, and the push–pull network forms an exponential filter (Hinczewski and Thirumalai, 2014; Becker et al., 2015). This underscores the idea that our theory gives a complete description, in terms of all the required resources, for the optimal design of cellular sensing systems that need to estimate this type of signals. Furthermore, because our model is Gaussian, the goal of minimizing the mean-square error in the estimate of the input signal is equivalent to maximizing the mutual information between the input (the ligand concentration) and the output (the readout x*) (Becker et al., 2015).

In recent years, filtering theories and information theory have been applied increasingly to neuronal and cellular systems (Laughlin, 1981; Brenner et al., 2000; Fairhall et al., 2001; Andrews et al., 2006; Ziv et al., 2007; Nemenman et al., 2008; Cheong et al., 2011; Nemenman, 2012; Hinczewski and Thirumalai, 2014; Becker et al., 2015; Husain et al., 2019; Tkacik et al., 2008; Tkačik and Walczak, 2011; Dubuis et al., 2013; Monti and Wolde, 2016; Monti et al., 2018a). A key concept in these theories is that optimal sensing systems match the response to the statistics of the input. When the noise is weak, maximizing the entropy of the output distribution becomes paramount, which entails matching the shape of the input–output relation to the shape of the input distribution to generate a flat output distribution (Laughlin, 1981; Tkacik et al., 2008; Monti et al., 2018a). Yet, when the noise is large, the optimal response is also shaped by the requirement to tame the propagation of noise in the input signal (Andrews et al., 2006; Hinczewski and Thirumalai, 2014; Becker et al., 2015; Monti et al., 2018a; Monti et al., 2018b; Mora and Nemenman, 2019) or to lift the signal above the intrinsic noise in the response system (Tostevin and ten Wolde, 2010; Bowsher et al., 2013). In Appendix 3, we show that estimating the concentration from pτr is equivalent to that via readout x*. This makes it possible to connect our sampling framework, which is based on pτr(L), to filtering and information theory, which are based on x*(L). In particular, we show in this appendix how the optimal integration and dynamic gain can be understood from these ideas on matching the response to the input. We also briefly discuss in Appendix 3 the concepts from information theory that are beyond the scope of the Gaussian model considered here.

Yet, our discrete sampling framework gives a detailed description of how the optimal design of sensing systems depends on the statistics of the input signal in terms of all the required cellular resources: protein copies, time, and energy. In an optimal system, each receptor is sampled once every receptor–ligand correlation time τc, Δτc, and the number of samples per receptor is τropt/Δτropt/τc. The optimal integration time τropt for a given RT is determined by the trade-off between the age of the samples and the number required for averaging the receptor state. When the input varies more rapidly, the samples need to be refreshed more regularly: to keep the dynamical error and the dynamic gain constant, τropt must decrease linearly with τL (see Equation 6). Yet, only decreasing τropt would inevitably increase the sampling error σp^τr2,samp in estimating the receptor occupancy because the sampling interval ΔRTτropt/XTopt would become smaller than τc, creating redundant samples. To keep the sensing precision constant, the number of receptors RT needs to be raised with τL-1, such that the sampling interval ΔRTτropt/XTopt remains of order τc and the decrease in the number of samples per receptor, τropt/τc, is precisely compensated for by the increase in RT. The total number of independent concentration measurements, RTτropt/τc, and hence the number of readout molecules XTopt to store these, does indeed not change. In contrast, the required power βw˙optRT/τc rises (Equation 12): each receptor molecule is sampled each τc at Δμopt4kBT, and the increase in RT raises the sampling rate n˙=w˙opt/ΔμoptXTopt/τropt. Our theory thus predicts that when the input varies more rapidly the number of receptors and the power must rise to maintain a required sensing precision, while the number of readout molecules does not.

The fitness benefit of a sensing system does not only depend on the sensing precision but also on the energetic cost of maintaining and running the system. In principle, the cell can reduce the sensing error arbitrarily by increasing RT and decreasing τr. Our resource allocation principle (Equation 12) shows that then not only the number of readout molecules needs to be raised but also the power. Clearly, improving the sensing precision comes at a cost: more copies of the components of the sensing system need to be synthesized every cell cycle, and more energy is needed to run the system. Our theory (i.e., Equation 6) makes it possible to derive the Pareto front that quantifies the trade-off between the maximal sensing precision and the cost of making the sensing system (see Figure 5). Importantly, the design of the optimal system at the Pareto front obeys, to a good approximation, our resource allocation principle (Equation 12). This is because this principle specifies the optimal ratios of RT, XT, w˙, and τr given the input statistics, and these ratios are fairly insensitive to the costs of the respective resources: resources that are in excess cannot improve sensing and are thus wasted, no matter how cheap they are. It probably explains why our theory, without any fit parameters, not only predicts the integration time that allows E. coli to sense shallow gradients (Figure 4) but also the number of receptor and readout molecules (Govern and Ten Wolde, 2014a).

The benefit of a sensing system depends on the sensing precision it can achieve and the cost of making it.

The Pareto front characterizes the trade-off between the maximal sensing precision, quantified by the maximal mutual information Imax(x*;L), and the cost of making the sensing system, C=RT+cXXT, where cX is the relative cost of making a readout versus a receptor protein, here taken to be cX=1. System designs below the Pareto front are suboptimal and can be improved by reducing the cost, that is, the number of proteins, and / or increasing the sensing precision. The optimal systems at the Pareto front obey, to a good approximation, the allocation principle Equation 12. The Pareto front, formed by the maximal value Imax(x*;L) of I(x*;L)=1/2ln(1+SNR) as a function of C, is obtained by minimizing Equation 6 over p,τr,RT,XT subject to the constraint C=RT+XT; the quality parameter is qopt0.76 corresponding to Δμopt4kBT; τc/τL=102; σL/L¯T=102.

In our study, we have limited ourselves to a canonical push–pull motif. Yet, the work of Govern and Ten Wolde, 2014a indicates that our results hold more generally, pertaining also to systems that employ cooperativity, negative or positive feedback, or multiple layers, as the MAPK cascade. While multiple layers and feedback change the response time, they do not make time integration more efficient in terms of readout molecules or energy (Govern and Ten Wolde, 2014a). And provided it does not increase the input correlation time (Skoge et al., 2011; Ten Wolde et al., 2016), cooperative ligand binding can reduce the sensing error per sample, but the resource requirements in terms of readout molecules and energy per sample do not change (Govern and Ten Wolde, 2014a). In all these systems, time integration requires that the history of the receptor is stored, which demands protein copies and energy.

Lastly, in this article we have studied the resource requirements for estimating the current concentration via the mechanism of time integration. However, to understand how E. coli navigates in a concentration gradient, we do not only have to understand how the system filters the high-frequency ligand-binding noise via time averaging but also how on longer timescales the system adapts to changes in the ligand concentration (Sartori and Tu, 2011). This adaptation system also exhibits a trade-off between accuracy, speed, and power (Lan et al., 2012; Sartori and Tu, 2015). Intriguingly, simulations indicate that the combination of sensing and adaptation allows E. coli not only to accurately estimate the current concentration but also the future ligand concentration (Becker et al., 2015). It will be interesting to see whether an optimal resource allocation principle can be formulated for systems that need to predict future ligand concentrations.

Materials and methods

Methods are described in Appendices 1–3. Appendix 1 derives the central result of our article (Equation 6). Appendix 2 derives the fundamental resources and the corresponding sensing limits (Equations 8 and 9). Appendix 3 describes how the optimal gain and integration time can be understood using ideas from filtering and information theory.

Appendix 1

Signal-to-noise ratio

Here, we provide the derivation of the central result of this article, Equation 6 of the main text. The derivation starts from the SNR, given in Equation 2. Here, σL2 is the width of the input distribution, while (δL^)2 is the error in the estimate of the concentration. The latter is derived from the dynamic input–output relation pτr(L), which is the mapping between the average receptor occupancy over the past integration time τr and the current ligand concentration L (see Figure 2). Concretely, the error (δL^)2 is given by Equation 1, where σp^τr2 is the error in the estimate of the average receptor occupancy over the past integration time τr and g~Lpτr is the dynamic gain, which is the slope of the dynamic input–output relation pτr(L). Below, we first derive the dynamic gain g~Lpτr and then the error in the estimate of the receptor occupancy σp^τr2.

Dynamic input–output relation

The dynamic input–output relation pτr(L) is the average receptor occupancy pτr over the past integration time τr, given that the current ligand concentration L(t)=L. The cell estimates pτr via its receptor readout system, which is a device that takes samples of the receptor: the readout molecules at time t constitute samples of the ligand-binding state of the receptor at earlier sampling times ti (see Figure 2). More specifically, the cell estimates pτr from the number of active readout molecules x*(L(t))=x*(L):

(14) p^τr(L)=x*(L)N¯,

where N¯ is the average of the number of samples N taken during the integration time τr. Hence, the dynamic input–output relation is

(15) pτr(L)En(ti)L(t),

where n(ti)=0,1 is the receptor occupancy at time ti, E denotes the expectation over the sampling times ti, and L(t) denotes an average over receptor–ligand binding noise and the subensemble of ligand trajectories that each end at L(t) (see Figure 2c); the quantity n(ti)L(t) is indeed the average receptor occupancy at time ti, given that the ligand concentration at time t is L(t)=L. Importantly, the receptor samples can also decay via the deactivation of x*. Taking this into account, the probability that a readout molecule at time t provides a sample of the receptor at an earlier time ti is p(ti|sample)=e-(t-ti)/τr/τr (Govern and Ten Wolde, 2014a). Averaging the receptor occupancy over the sampling times ti then yields

(16) pτr(L)=-t𝑑tin(ti)L(t)e-(t-ti)/τrτr.

Dynamic gain

When the current ligand concentration L(t) deviates from its mean L¯ by δL(t)L(t)-L¯, then pτr deviates on average from its mean p (the average receptor occupancy over all δL(t)) by

(17) δpτrpτr-p=Eδn(ti)δL(t)=-t𝑑tiδn(ti)δL(t)e-(t-ti)/τrτr.

Here, E denotes again the expectation over the sampling times ti, and δn(ti)δL(t)n(ti)δL(t)-p is the average deviation in the receptor occupancy n(ti) at time ti from its mean p, given that the ligand concentration at time t is δL(t) (see Figure 2c). We can compute it within the linear-noise approximation (Gardiner, 2009):

(18) δn(ti)δL(t)=ρn-ti𝑑tδL(t)δL(t)e-(ti-t)/τc,

where ρn=p(1-p)/(L¯Tτc) and δL(t)δL(t) is the average ligand concentration at time t, given that the ligand concentration at time t is δL(t). The latter is given by Bowsher et al., 2013

(19) δL(t)δL(t)=δL(t)e-|t-t|/τL.

Combining Equations 17–19 yields the following expression for the average change in the average receptor occupancy pτr, given that the ligand at time t is δL(t):

(20)δpτr=p(1p)L¯T(1+τcτL)1(1+τrτL)1δL(t),(21)g~Lpτr δL(t).

Hence, the dynamic gain is

(22)g~Lpτr=p(1p)L¯(1+τcτL)1(1+τrτL)1,(23)=gLp(1+τcτL)1(1+τrτL)1.

The dynamic gain is the slope of the dynamic input–output relation pτr(L) (see Figure 2a). It yields the average change in the receptor occupancy pτr over the past integration time τr when the change in the ligand concentration at time t is δL(t). It depends on all the timescales in the problem and only reduces to the static gain gLp=p(1-p)/L¯ when the integration time τr and the receptor correlation time τc are both much shorter than the ligand correlation time τL. The dynamic gain determines how much an error in the estimate of pτr propagates to the estimate of L(t).

Error in receptor occupancy

We can derive the variance in the estimate of the receptor occupancy over the past integration time τr, σp^τr2, directly from Equation 14 for the system in the irreversible limit (Malaguti and Ten Wolde, 2019). While this derivation is illuminating, it is also lengthy. For the fully reversible system studied here, we follow a simpler route. Since the average number of samples N¯ over the integration time τr is constant, it follows from Equation 14 that

(24) σp^τr2=σx*|L2N¯2,

where σx*|L2 is the variance in the number of phosphorylated readout molecules, conditioned on the signal at time t being L(t)=L. The conditional variance (Tostevin and ten Wolde, 2010)

(25) σx*|L2=σx*2-g~Lx*2σL2

is the full variance σx*2 of x* minus the variance g~Lx*2σL2 that is due to the signal variations, given by the dynamic gain g~Lx*2 from L to x* times the signal variance σL2.

The full variance of the readout σx*2 in Equation 25 can be obtained from the linear-noise approximation (Gardiner, 2009), see Malaguti and Ten Wolde, 2019:

(26) σx2=f(1f)XT+ρ2μ(μ+μ)[p(1p)RT+ρ2σL2(λ+μ+μ)μ(λ+μ)(λ+μ)].

In this expression, μ=τc-1=k1L¯+k2 is the inverse of the receptor correlation time τc; p=RL¯/RT=k1L¯/(k2+k1L¯)=k1L¯τc is the probability that a receptor is bound to ligand; ρ=RTk1(1p)=p(1p)RTμ/L¯; μ=τr1=(kf+kf)pRT+kr+kr is the inverse of the integration time τr; f=x¯/xT=(kfpRT+kr)τr is the fraction of phosphorylated readout; and ρ=kfXT(1-f)-k-fXTf=n˙/(pRT) is the total flux n˙ around the cycle of readout activation and deactivation divided by the total number pRT of ligand-bound receptors: it is the rate at which each receptor is sampled, be it ligand bound or not. For what follows below, we note that the quality parameter q=(eΔμ1-1)(eΔμ2-1)/(eΔμ-1)=ρpRTτr/(f(1-f)XT)=n˙τr/(f(1-f)XT).

To get σp^τr2 from Equations 24 and 25, we need not only σx2 (Equation 26) but also the average number of samples N¯ and the dynamic gain g~Lx*2. The average number of samples taken during the integration time τr is N¯=n˙τr/p=f(1-f)XTq/p=ρRT/μ, and the effective number of reliable samples is N¯eff=qN¯. Since pτr(L)=Ex*L/N¯, where Ex*L is the average number of active readout molecules for a given input L(t)=L and N¯ is a constant independent of L, it follows that

(27) g~Lx*=g~LpτrN¯=g~LpτrRTρμ,

with g~Lpτr the dynamic gain from L to pτr, given by Equation 22. Equation 27 can be verified via another route that does not rely on the sampling framework because we also know that g~Lx*=σL,x*2/σL2 (Tostevin and ten Wolde, 2010), where the co-variance σL,x*2 can be obtained from the linear-noise approximation (Malaguti and Ten Wolde, 2019; Gardiner, 2009). Combining Equations 24–27 yields

(28) σp^τr2=p(1p)N¯eff+p(1p)RT(1+τr/τc)+p2N¯eff+g~Lpτr2σL2[(1+τcτL)(1+τrτL)(1+τcτrτL(τc+τr))1].

This can be rewritten using the expression for the fraction of independent samples, which, assuming that τrτc, is fI=1/(1+2τc/Δ), with Δ=2τrRT/N¯eff the effective spacing between the samples (Govern and Ten Wolde, 2014a):

(29) σp^τr2=p(1-p)fIN¯eff+p2N¯effσp^τr2,samp+g~Lpτr2σL2[(1+τcτL)(1+τrτL)(1+τcτrτL(τc+τr))-1]σp^τr2,dyn,

Here, σp^τr2,samp is the sampling error in the estimate of pτr (Malaguti and Ten Wolde, 2019); it is a statistical error, which arises from the finite cellular resources to sample the state of the receptor, protein copies, time, and energy (see Figure 2b). The other contribution, σp^τr2,dyn, is the dynamical error in the estimate of pτr (Malaguti and Ten Wolde, 2019); it is a systematic error that arises from the input dynamics and only depends on the average receptor occupancy and the timescales of the input, receptor, and readout (see Figure 2c); it neither depends on the number of protein copies nor on the energy necessary to sample the receptor.

Final result: SNR

Combining Equations 29 and 22 with Equation 3 yields the principal result of our work (Equation 6) of the main text.

Appendix 2

Fundamental resources

To identify the fundamental resources limiting the sensing accuracy and derive the corresponding sensing limits (Equations 8 and 9), it is helpful to rewrite the SNR in terms of collective variables that illuminate the cellular resources. For that, we start from Equation 6 of the main text and split the first term on the right-hand side and exploit the expression for the effective number of independent samples N¯I=1/(1+2τc/Δ)N¯eff with Δ=2τrRT/N¯eff. We then sum up the last two terms on the right-hand side and use that N¯eff=qN¯=qn˙τr/p:

SNR1=(1+τcτL)2(1+τrτL)2[(L¯/σL)2N¯effp(1p)2+(L¯/σL)2p(1p)RT(1+τr/τc)](30)+(1+τcτL)(1+τrτL)(1+τcτrτL(τc+τr))1=(1+τcτL)2(1+τrτL)2[(L¯/σL)2(1p)2qn˙τrcoding noise+(L¯/σL)2p(1p)RT(1+τr/τc)receptor input noise](31)+(1+τcτL)(1+τrτL)(1+τcτrτL(τc+τr))1dynamical error.

The second term in between the square brackets describes the contribution to the sensing error that comes from the stochasticity in the concentration measurements at the receptor level. The first term in between the square brackets, the coding noise, describes the contribution that arises in storing these measurements into the readout molecules.

From Equation 30, the fundamental resources and the corresponding sensing limits (Equations 8 and 9) can be derived. Specifically, when the number of receptors and their integration are limiting, the coding noise in Equation 30 is zero; exploiting that typically τcτr,τL and that the contribution to the sensing error from the receptor input noise is minimized for p1/2, this yields Equation 8 with h=RTτr/τc. When the number of readout molecules XT is limiting, the receptor input noise is zero and q1; noting that n˙=f(1-f)XTq/τr and that the contribution from the coding noise is minimized when f1/2 and p0, and again exploiting that τcτr,τL, this yields Equation 8 with h=XT. When the power w˙=n˙Δμ is limiting, then the receptor input noise is (again) zero. The coding noise is minimized for a given power constraint w˙ when Δμ1=Δμ2=Δμ/2, but two regimes can be distinguished based on the total free-energy drop Δμ. When Δμ>4kBT, the system is in the irreversible regime and q1 (see Equation 7); Equation 30 shows that the error is then bounded by Equation 8 with h=w˙τr/(Δμ/4), using τcτr,τL and p0. Yet, the sensing error is minimized in the quasi-equilibrium regime, where Δμ1=Δμ2=Δμ/20 and qβΔμ/4, yielding Equation 8 with h=βw˙τr.

Appendix 3

The optimal gain and optimal integration time

The theory of the main text (Equation 6) is based on the idea that the cell uses its push–pull network to estimate the receptor occupancy pτr(L) from which the current ligand concentration L is then inferred by inverting the dynamic input–output relation pτr(L). Yet, as we show here, this framework is equivalent to the idea that the cell estimates the concentration from the output x*, using the dynamic input–output relation x*(L). Here, we use this observation to analyze our system using ideas from filtering and information theory. But first we demonstrate this correspondence.

To show that estimating the concentration from p^τr is equivalent to that from estimating it from x*, we first note that because the average number of samples N¯ is constant, σx*|L2=σp^τr2N¯2 while the gain from L to x* is g~Lx*2=g~Lpτr2N¯2. Consequently, the absolute error (δL^)2 in estimating the concentration via x*, (δL^)2=σx*|L2/g~Lx*2, is the same as that of Equation 1: because the instantaneous number of active readout molecules x* reflects the average receptor occupancy pτr over the past τr, estimating the ligand concentration from x* is no different from inferring it from the average receptor occupancy p^τr=x*/N¯.

To make the connection with information and filtering theory, we note that in our Gaussian model the conditional distribution of δx* given δL is given by Tostevin and ten Wolde, 2010

(32) p(δx|δL)=12πσx|L2e(δxg~LxδL)22σx|L2,

where g~Lx*δL=δxL is the average value of δx* given that δL(t)=δL, and σx*|L2 is the variance of this distribution (see also Equation 25).

The relative error, the inverse of the SNR (see Equation 2), is

(33) SNR1=(δL^)2σL2=σx|L2g~Lx2σL2.

As mentioned in the main text, the SNR also yields the mutual information I(x*;L)=1/2ln(1+SNR) between the input L and output x* (Tostevin and ten Wolde, 2010).

The notion of an optimal integration time or optimal dynamic gain is well known from filtering and information theory (Andrews et al., 2006; Hinczewski and Thirumalai, 2014; Becker et al., 2015; Monti et al., 2018a; Monti et al., 2018b; Mora and Nemenman, 2019). To elucidate the optimal gain and integration time in our system, we combine the above equation with Equations 25 and 26 to write the relative error as

(34) SNR1=f(1f)XTg~Lx2σL2readout switching noise+gRLx21/(1+τr/τc)p(1p)RTg~Lx2σL2receptor input noise+(1+τcτL)(1+τrτL)(1+τcτrτL(τc+τr))1dynamical error,

where gRLx*=ρ/μ is the static gain from RL to x*. Written in this form, the trade-offs in maximizing the mutual information I(x*;L) (and minimizing the relative error in estimating the concentration) become apparent: increasing the dynamic gain g~Lx* by decreasing the integration time τr raises the slope of the input–output relation x*(L), which helps to lift the transmitted signal above the intrinsic binomial switching noise of the readout, f(1-f)XT. Also, the dynamical error is minimized by minimizing τr and maximizing g~Lx*. Yet, for the second term, which describes how noise in the input signal arising from receptor switching, p(1-p)RT, is propagated to the output x*, there exists an optimal integration time that minimizes this term: while decreasing τr increases the dynamic gain, which helps to raise the signal above the noise, it also impedes time averaging of this switching noise, described by the factor 1/(1+τr/τc).

The mutual information is I(x*;L)=H(x*)-H(x*|L), with H(x*) the entropy of the marginal output distribution and H(x*|L) the entropy of the output distribution conditioned on the input. Hence, information theory shows that in the weak noise limit, information transmission is optimal when the entropy of the output distribution is maximized (Laughlin, 1981; Tkacik et al., 2008). Our system obeys this principle. Since the dynamic gain g~Lx*=ρρτL2τcτr/[(τc+τL)(τr+τL)]RTXT, the amplification of the signal rises with RT and XT. Since the standard deviation of the noise added to the transmitted signal coming from the stochastic receptor and readout activation scales with RT and XT, respectively, it is clear that the SNR increases with RT and XT. In the limit that RT,XT, the relative error SNR-1 is only set by the dynamical error, which can be reduced to zero by τr0, exploiting that typically τcτL. This is the weak-noise limit in which the mutual information I(x*;L) is maximized by maximizing the entropy of the output distribution H(x*). Indeed, τr0 corresponds to maximizing the gain, which maximizes the width of the output distribution, in this limit equal to σx2=g~Lx*2σL2 (see Equation 25), and thereby the entropy of the output distribution H(x*)=1/2ln(2πeσx2).

Finally, we note that our Gaussian model is linear such that the central control parameter, besides protein copies and energy, is the integration time or the dynamic gain, which sets the slope of the linear input–output relation. While Wiener–Kolmogorov and Kalman filtering are exact only for these Gaussian models, information theory also applies to non-linear systems with non-Gaussian statistics. It has been used to show that neuronal systems (Laughlin, 1981; Brenner et al., 2000; Fairhall et al., 2001; Nemenman et al., 2008; Tkacik et al., 2010), signaling and gene networks (Segall et al., 1986; Tkacik et al., 2008; Tkačik and Walczak, 2011; Nemenman, 2012; Dubuis et al., 2013), and circadian systems (Monti and Wolde, 2016; Monti et al., 2018a) can maximize information transmission by optimizing the shape of the input–output relation (Laughlin, 1981; Brenner et al., 2000; Fairhall et al., 2001; Tkacik et al., 2008; Monti et al., 2018a); by desensitization, that is, adapting the output to the mean input via incoherent feedforward or negative feedback (Segall et al., 1986); by gain control, that is, adapting the output to the variance of the input by capitalizing on a steep response function and temporal correlations in the input (Nemenman, 2012); by removing coding redundancy via temporal decorrelation (Nemenman et al., 2008); by optimizing the tiling of the output space via the topology of the network (Tkačik and Walczak, 2011; Dubuis et al., 2013); or by exploiting cross-correlations between the signals (Tkacik et al., 2010; Monti and Wolde, 2016).

Data availability

All data generated or analysed during this study are included in the manuscript and supporting files.

References

  1. Book
    1. Extrapolation WN
    (1950)
    Extrapolation, Interpolation, and Smoothing of Stationary Time Series: With Engineering Applications
    MIT Press.
  2. Book
    1. Gardiner CW
    (2009)
    Stochastic Methods: A Handbook for the Natural and Social Sciences
    Berlin: Springer-Verlag.
  3. Book
    1. Kolmogorov AN
    (1992) Probability theory and mathematical statistics
    In: Watanabe S, Prokhorov J. V, editors. Selected Works of A. N. Kolmogorov. Netherlands: Springer Science & Business Media. pp. 8–14.
    https://doi.org/10.1007/BFb0078455

Decision letter

  1. Raymond E Goldstein
    Reviewing Editor; University of Cambridge, United Kingdom
  2. Detlef Weigel
    Senior Editor; Max Planck Institute for Developmental Biology, Germany

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

Malaguti and ten Wolde combine analysis of sensing limits for time-varying signals with previous work on resource allocation to minimize the costs of sensing. The work makes a solid contribution towards understanding the principles of resource allocation, and the conclusion that E. coli chemotaxis is optimized for shallow gradients should stimulate further discussion and work.

Decision letter after peer review:

Thank you for submitting your article "Theory for the optimal detection of time-varying signals in cellular sensing systems" for consideration by eLife. Your article has been reviewed by two peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by Detlef Weigel as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

Malaguti and ten Wolde combine analysis of Mora-Nemenman (PRL 2019) on Berg-Purcell type sensing limits for time-varying signals with previous work of one of the authors (ten Wolde, PNAS 2014) on resource allocation to minimize the costs of sensing. The previous work discussed tradeoffs for the costs of sensing a constant signal. Here it is extended to time-varying signals (modeled as colored-noise Gaussian). Although the authors exaggerate the novelty of their analysis, they make a solid contribution towards understanding the principles of resource allocation. The conclusion that E. coli chemotaxis is optimized for shallow gradients seems reasonable and should stimulate further discussion and work. Despite flaws in establishing context, the paper deserves publication

Revisions for this paper:

While enthusiastic about the results in the manuscript, the reviewers found that there are issues involving claims of novelty, and the overall presentation of mathematical results that will need considerable revision.

1) The novelty of their analysis is exaggerated.

The basic problem of optimally estimating the state of a time-varying signal probed by noisy measurements is textbook material for engineers. Here, because everything is linearized about an operating point and because noise is approximated as Gaussian, there is an exact analytic theory going back to Kalman and Bucy in 1961. The result underlies the Mora-Nemenman analysis (although they, too, did not reference that work for the Gaussian approximation used for their concrete results) and that done here. In engineering, the general problem optimally estimating a stochastic signal via noisy measurements was already considered by Kolmogorov and Wiener in the 1940s and formulated in the time domain, for linear systems driven by white noise (the approximation used here) by Kalman and Bucy.

It is a weakness that the connections are not mentioned specifically by referencing. Appealing to known results not only underlines how different disciplines often need to tackle the same problems, it allows for the use of textbook results and can shorten a paper. As a corollary, authors who map their problem onto known results should not be penalized for doing so, when the application is new and important (as here).

2) Specific sentences have language that seems exaggerated, in the light of historical work on filtering theory:

"Our theory is based on a new concept, the dynamic input-output relation pτr(L)."

- Dynamical models (with dynamical input-output relations) are a central element of filtering theory.

".Our theory reveals that the sensing error can be decomposed into two terms [sampling (sensing) error and dynamical error]."

- The framework of filtering theory assumes noisy measurements of stochastic signals.

"Our theory illuminates how the optimal design depends on the timescale of the input τL."

- The statement is true but should be framed in the context of other work, including in systems biology / biophysics, which comes to the same conclusion. For example, the work of Laughlin in 1981, and later Nemenman and Bialek, show that the input-output "gain function" should be adapted to the statistics of the input (including time scales of variation). Bialek's 2012 Biophysics book discusses many examples.

3) In its present form, the paper will be essentially unreadable by the vast majority of the eLife readership; the writing style is more appropriate to the Physical Review than to eLife. Indeed, much of the paper is written like a technical note that builds on previous work (Berg/Purcell, etc.) without sufficient explanation for the general reader. This criticism extends to the explanations of the underlying model, the lack of definition of terminology, and the notation.

Here are but a few examples of the points above. The authors are encouraged to rethink the entire presentation de novo for maximum improvement.

a) The definition of a "push-pull" network should be given in the paper.

b) For the purpose of readability by a diverse audience some care should be taken to define technical terms explicitly (e.g. "Markovian"), and to explain various statements more clearly. Appearing so early in the paper, such condensed statements will be off-putting to the general reader.

c) The whole notation of "inverting the mapping" that is central to the theory is not well-explained.

d) Regarding notation: consider, for example, the caption of Figure 2 and Equation 1, in which there is the quantity σ2p̂tr|L . This is simply too complicated, and its meaning will be utterly opaque to the general reader.

4) In many ways it seems that if the authors presaged the development of the theory by indicating the run-and-tumble context early on it would help the reader to understand the precise motivations behind their lengthy calculations, and to give an idea of the time scales.

https://doi.org/10.7554/eLife.62574.sa1

Author response

Revisions for this paper:

While enthusiastic about the results in the manuscript, the reviewers found that there are issues involving claims of novelty, and the overall presentation of mathematical results that will need considerable revision.

1) The novelty of their analysis is exaggerated.

The basic problem of optimally estimating the state of a time-varying signal probed by noisy measurements is textbook material for engineers. Here, because everything is linearized about an operating point and because noise is approximated as Gaussian, there is an exact analytic theory going back to Kalman and Bucy in 1961. The result underlies the Mora-Nemenman analysis (although they, too, did not reference that work for the Gaussian approximation used for their concrete results) and that done here. In engineering, the general problem optimally estimating a stochastic signal via noisy measurements was already considered by Kolmogorov and Wiener in the 1940s and formulated in the time domain, for linear systems driven by white noise (the approximation used here) by Kalman and Bucy.

It is a weakness that the connections are not mentioned specifically by referencing. Appealing to known results not only underlines how different disciplines often need to tackle the same problems, it allows for the use of textbook results and can shorten a paper. As a corollary, authors who map their problem onto known results should not be penalized for doing so, when the application is new and important (as here).

We fully acknowledge the point that is made here and we agree that we should have described the connection with filtering theory much more clearly, especially given the fact that we have applied filtering theory ourselves to biochemical signalling before (Becker, Mugler, Ten Wolde, 2015).

We indeed employ a Gaussian model, and for such a Gaussian model Kalman and Wiener-Kolmogorov filtering theory become exact, as the reviewers correctly point out. In our previous work, we studied different types of input signals (Markovian and non-Markovian), and then used Wiener-Kolmogorov filtering theory to derive the optimal topology of the sensing network for the different types of input signals (Becker et al., 2015). In our current work, we consider one type of input signal—a stationary Markovian signal. The optimal filter for this type of input signal is a low-pass, exponential filter and the canonical network motif studied in our current study implements precisely such a filter (see Becker et al., 2015). This filter is useful, because it allows the system to time integrate the input signal and hence filter out the high-frequency input noise arising here from receptor-ligand binding. We agree with the reviewers that these ideas are well established and that we should describe them more clearly in our manuscript.

Yet, while filtering theory provides a powerful approach to elucidating the optimal topology and response dynamics of the sensing network, as determined by the statistics of the input signal, it does not naturally reveal the resource requirements for sensing. To this end, we have generalised the sampling framework that we have developed previously (Govern and Ten Wolde, 2014). This framework is particularly well suited for elucidating the resource requirements for time integration, because it starts from the observation that the receptor-readout system consists of discrete molecules that interact with the receptor in a stochastic fashion: our description is thus based on the idea that the readout system is a sampling device that implements the mechanism of time integration not by continuously integrating the state of the receptor, but rather by discretely and stochastically sampling it, via the collisions of the readout molecules with the receptor proteins. The novelty of this manuscript is that we have now extended this sampling framework, for the first time, to time-varying signals.

In the revised manuscript, we have thoroughly rewritten the Introduction and added two new paragraphs to the Discussion. In the Introduction, we now mention that the problem of optimally predicting time-varying signals is a classic problem, for which different theories have been developed. We then also briefly review the application of filtering theories to biochemical signalling, adding relevant references. We then emphasise that we consider one class of input signals and one canonical network motif, which is known to implement a filter that is optimal for this type of signals, and raise the central question of our manuscript, which is what the cellular resource requirements are for optimally implementing such a network. We end the Introduction with a brief discussion of the main ideas of our theory and our findings.

In the revised Introduction, we also note that our theory applies to systems that employ the mechanism of time integration, and not to systems that employ Maximum-Likelihood sensing or Bayesian filtering, which underlies the analysis of Mora and Nemenman, 2019.

In the new paragraphs in the Discussion section, we first emphasise that our model is Gaussian and that the performance criterion of our theory, minimizing the mean-square error, is identical to that of Wiener and Kalman filtering theory, which are exact for Gaussian models. We then mention that for our Gaussian model minimizing the mean-square error is equivalent to maximizing the mutual information, thus making a connection between our theory, filtering theories, and information theory. In the next paragraph, we then discuss the concepts that have emerged from these theories and we describe how our findings relate to these concepts; here we have also added a new Appendix 3, where we work out this connection in more mathematical detail.

We thank the reviewers and the editor for this important comment, because we agree that addressing it puts our work into a broader context.

2) Specific sentences have language that seems exaggerated, in the light of historical work on filtering theory:

"Our theory is based on a new concept, the dynamic input-output relation pτr(L)."

- Dynamical models (with dynamical input-output relations) are a central element of filtering theory.

We fully agree that dynamic input-output relations are a key concept in filtering theory and we acknowledge that similar input-output relations have been studied before in the biophysics literature (Bialek et al., IEEE, 2006; Tostevin and Ten Wolde, 2010; Nemenman, 2012; Hinczewski and Thirumalai, 2014; Becker et al., 2015). We merely wanted to emphasise that our dynamic input-output relation pτr(L) differs fundamentally from the conventional static input-output relations that are typically considered in the context of sensing static signals (Berg and Purcell, 1977; Bialek and Setayesghar, 2005; Kaizu et al., 2014; Mugler et al., 2016).

We have completely rewritten this paragraph. As we describe in our response to the previous point, we now first introduce filtering theory, then describe the novel aspect of the current manuscript, which is the extension of our sampling framework to time-varying signals, and then mention the dynamic input-output relation, contrasting it with the static input-output relation used to describe the sensing of static signals.

"Our theory reveals that the sensing error can be decomposed into two terms [sampling (sensing) error and dynamical error]."

- The framework of filtering theory assumes noisy measurements of stochastic signals.

We agree with this, but our decomposition is more specific, showing how the sampling error and the dynamical error depend on the cellular resources protein copies, time, and energy. Our decomposition identifies the combinations of cellular resources, the “collective variables”, that control the sensing precision. It is this decomposition that allows us to make specific predictions on the optimal design of sensing systems that maximize the sensing precision given resource constraints. This is the novel contribution of our paper.

We now emphasise this more clearly.

"Our theory illuminates how the optimal design depends on the timescale of the input τL."

- The statement is true but should be framed in the context of other work, including in systems biology / biophysics, which comes to the same conclusion. For example, the work of Laughlin in 1981, and later Nemenman and Bialek, show that the input-output "gain function" should be adapted to the statistics of the input (including time scales of variation). Bialek's 2012 Biophysics book discusses many examples.

The idea of an optimal integration time, and the related ideas of matching the gain and the input-output relation to the statistics of the input, are indeed well known. In fact, we had cited a number of papers that also show the existence of an optimal integration time (Becker et al., 2015; Monti et al., 2018; Mora and Nemenman, 2019). The aim of our paragraph in the Discussion, was to discuss how the input timescale governs the optimal design in terms of the cellular resources receptors, readout molecules, and energy—this is the novel aspect of our work. But we agree with the reviewers that it is important to connect our observations to previous ideas on the optimal gain and input-output relations that are matched to the statistics of the input. We therefore now first discuss here, in the Discussion section, these ideas and then describe how our findings relate to them. After we have discussed this connection, we emphasise that our sampling framework gives a detailed description of the optimal design in terms of the required cellular resources, and then discuss how this depends on the statistics of the input signal.

More concretely, we first point out that our model is Gaussian and that the performance criterion of our theory is minimizing the mean-squared error, which is precisely the performance criterion of Wiener and Kalman filtering. We emphasise that the push-pull network is an exponential filter and that this is the optimal filter for the memoryless (Markovian) signals as studied here. We also point out that for our Gaussian model minimizing the mean-square error is equivalent to maximizing the mutual information, thus making the connection between our work, filtering theory, and information theory. In the next paragraph, we then describe concepts that have emerged from filtering and information theory, and describe how our findings relate to these concepts. Here, we also refer to a new Appendix 3, in which we work this connection out in more (mathematical) detail. In this appendix we also discuss concepts for optimizing information transmission in non-linear systems, which are beyond the scope of the current model.

3) In its present form, the paper will be essentially unreadable by the vast majority of the eLife readership; the writing style is more appropriate to the Physical Review than to eLife. Indeed, much of the paper is written like a technical note that builds on previous work (Berg/Purcell, etc.) without sufficient explanation for the general reader. This criticism extends to the explanations of the underlying model, the lack of definition of terminology, and the notation.

Here are but a few examples of the points above. The authors are encouraged to rethink the entire presentation de novo for maximum improvement.

We acknowledge that the audience of eLife is broader than that of Physical Review and we have rewritten the manuscript drastically. We now try to avoid technical terms like “Markovian”, “mapping”, “push-pull” network as much as possible, and where this cannot be avoided we have explained them for a broader audience. Indeed, we have carefully verified whether all the necessary concepts that are needed to follow our analysis are explained for a broad audience. For example, in the description of our model and our theory we now describe in more detail the concept of the input correlation time, receptor correlation time, readout relaxation time, the idea of time integration via discrete sampling of the receptor state. We have also simplified notation. Most importantly, perhaps, we have tried to rewrite the manuscript in a more lucid style so that it is easier to follow for a broad audience (see, for example, the description of the optimal design principle, Equation 12). With these improvements, we strongly believe that the revised manuscript is accessible for the broad readership of eLife, also because the previous, well cited, paper on which the current manuscript builds, is written with a similar level of detail and published in a journal with a similarly broad audience (Govern and Ten Wolde, 2014).

a) The definition of a "push-pull" network should be given in the paper.

While we could replace “push-pull network” with “cellular network” or simply “network”, we think it does help to keep this term, because it specifically refers to the cycle of readout phosphorylation and dephosphorylation downstream of the receptor. We have therefore decided to keep this term. We do, however, now explain it in the Introduction when we introduce it for the first time.

b) For the purpose of readability by a diverse audience some care should be taken to define technical terms explicitly (e.g. "Markovian"), and to explain various statements more clearly. Appearing so early in the paper, such condensed statements will be off-putting to the general reader.

Since “Markovian” only appears twice in the manuscript, we have changed “Markovian” by “memoryless” with an explanation. The concepts of “correlation time”, “integration time”, and “relaxation time”, are important and used widely in the manuscript; we have therefore added explanations the first time we introduce these, see Section “Theory: Model”. We have also clarified other technical terms, such as time averaging, sampling device, coding.

c) The whole notation of "inverting the mapping" that is central to the theory is not well-explained.

We have clarified that.

d) Regarding notation: consider, for example, the caption of Figure 2 and Equation 1, in which there is the quantity σ2p̂tr|L This is simply too complicated, and its meaning will be utterly opaque to the general reader.

We have now simplified the notation as much as possible, while trying to keep clear to what specific quantity the symbol actually refers. For example, we have simplified σ2p̂tr|L to σ2p̂tr; we needed to retain the subscript τr because this quantity denotes the variance not of the receptor occupancy p but rather that of the receptor occupancy pτr over the integration time τr.

4) In many ways it seems that if the authors presaged the development of the theory by indicating the run-and-tumble context early on it would help the reader to understand the precise motivations behind their lengthy calculations, and to give an idea of the time scales.

We would like to emphasise that the analysis is much more generic than the specific application to the E. coli chemotaxis system. We therefore would like to maintain the distinction between the general theory and the application of it to E. coli chemotaxis. However, we also agree that introducing this system earlier helps to understand the motivation and main ideas behind our theory. We have therefore expanded the discussion of this system in the Introduction.

https://doi.org/10.7554/eLife.62574.sa2

Article and author information

Author details

  1. Giulia Malaguti

    AMOLF, Science Park, Amsterdam, Netherlands
    Contribution
    Conceptualization, Investigation, Methodology, Writing - original draft, Writing - review and editing
    Competing interests
    No competing interests declared
  2. Pieter Rein ten Wolde

    AMOLF, Science Park, Amsterdam, Netherlands
    Contribution
    Conceptualization, Resources, Funding acquisition, Investigation, Methodology, Writing - review and editing
    For correspondence
    tenwolde@amolf.nl
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-9933-4016

Funding

Nederlandse Organisatie voor Wetenschappelijk Onderzoek

  • Giulia Malaguti
  • Pieter Rein ten Wolde

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We wish to acknowledge Bela Mulder, Tom Shimizu, and Tom Ouldridge for many fruitful discussions and a careful reading of the manuscript. This work is part of the research program of the Netherlands Organisation for Scientific Research (NWO) and was performed at the research institute AMOLF.

Senior Editor

  1. Detlef Weigel, Max Planck Institute for Developmental Biology, Germany

Reviewing Editor

  1. Raymond E Goldstein, University of Cambridge, United Kingdom

Publication history

  1. Received: September 24, 2020
  2. Accepted: February 12, 2021
  3. Accepted Manuscript published: February 17, 2021 (version 1)
  4. Version of Record published: March 10, 2021 (version 2)

Copyright

© 2021, Malaguti and ten Wolde

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 712
    Page views
  • 145
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Developmental Biology
    2. Physics of Living Systems
    Yonghyun Song, Changbong Hyeon
    Research Article Updated

    Spatial boundaries formed during animal development originate from the pre-patterning of tissues by signaling molecules, called morphogens. The accuracy of boundary location is limited by the fluctuations of morphogen concentration that thresholds the expression level of target gene. Producing more morphogen molecules, which gives rise to smaller relative fluctuations, would better serve to shape more precise target boundaries; however, it incurs more thermodynamic cost. In the classical diffusion-depletion model of morphogen profile formation, the morphogen molecules synthesized from a local source display an exponentially decaying concentration profile with a characteristic length λ. Our theory suggests that in order to attain a precise profile with the minimal cost, λ should be roughly half the distance to the target boundary position from the source. Remarkably, we find that the profiles of morphogens that pattern the Drosophila embryo and wing imaginal disk are formed with nearly optimal λ. Our finding underscores the cost-effectiveness of precise morphogen profile formation in Drosophila development.

    1. Physics of Living Systems
    Tianfa Xie et al.
    Research Article

    A monolayer of highly motile cells can establish long-range orientational order, which can be explained by hydrodynamic theory of active gels and fluids. However, it is less clear how cell shape changes and rearrangement are governed when the monolayer is in mechanical equilibrium states when cell motility diminishes. In this work, we report that rat embryonic fibroblasts (REF), when confined in circular mesoscale patterns on rigid substrates, can transition from the spindle shapes to more compact morphologies. Cells align radially only at the pattern boundary when they are in the mechanical equilibrium. This radial alignment disappears when cell contractility or cell-cell adhesion is reduced. Unlike monolayers of spindle-like cells such as NIH-3T3 fibroblasts with minimal intercellular interactions or epithelial cells like Madin-Darby canine kidney (MDCK) with strong cortical actin network, confined REF monolayers present an actin gradient with isotropic meshwork, suggesting the existence of a stiffness gradient. In addition, the REF cells tend to condense on soft substrates, a collective cell behavior we refer to as the 'condensation tendency'. This condensation tendency, together with geometrical confinement, induces tensile prestretch (i.e., an isotropic stretch that causes tissue to contract when released) to the confined monolayer. By developing a Voronoi-cell model, we demonstrate that the combined global tissue prestretch and cell stiffness differential between the inner and boundary cells can sufficiently define the cell radial alignment at the pattern boundary.