Physical limit to concentration sensing in a changing environment

TL;DR

This paper establishes a new physical limit on the accuracy of cellular concentration sensing in fluctuating environments, showing it scales differently than previous bounds and can be achieved by adaptive biochemical networks.

Contribution

It formulates a non-linear field-theoretic model for sensing in fluctuating environments and derives a novel bound on sensing accuracy, contrasting with classical bounds.

Findings

Derived a new physical bound: δc/c ∼ (Dacτ)^(-1/4)

Showed the bound can be achieved by adaptive biochemical networks

Contrasted the new bound with the classical Berg-Purcell limit

Abstract

Cells adapt to changing environments by sensing ligand concentrations using specific receptors. The accuracy of sensing is ultimately limited by the finite number of ligand molecules bound by receptors. Previously derived physical limits to sensing accuracy have assumed that the concentration was constant and ignored its temporal fluctuations. We formulate the problem of concentration sensing in a strongly fluctuating environment as a non-linear field-theoretic problem, for which we find an excellent approximate Gaussian solution. We derive a new physical bound on the relative error in concentration $c$ which scales as $\delta c/c \sim (Dac\tau)^{-1/4}$ with ligand diffusivity $D$, receptor cross-section $a$, and characteristic fluctuation time scale $\tau$, in stark contrast with the usual Berg and Purcell bound $\delta c/c \sim (DacT)^{-1/2}$ for a perfect receptor sensing concentration during time $T$. We show how the bound can be achieved by a simple biochemical network downstream the receptor that adapts the kinetics of signaling as a function of the square root of the sensed concentration.

Equations

$$P_{\mathrm{prior}}(\{\phi (t)\})=\frac{1}{Z_{\mathrm{prior}}}\exp \left[-\frac{\tau }{2}{\int }_0^{T}dt{\left(\frac{d\phi }{dt}\right)}^2\right].$$

$$\begin{gathered} P(\{\phi (t)\})=\frac{P(t_1,\dots ,t_n|\{\phi (t)\})P_{\mathrm{prior}}(\{\phi (t)\})}{P(t_1,\dots ,t_n)} \\ =\frac{1}{Z}\exp \left\{-{\int }_0^{T}dt\left[\frac{\tau }{2}{\left(\frac{d\phi }{dt}\right)}^2+r_0{e}^{-\phi (t)}\right]-\sum _{i=1}^n\phi (t_i)\right\}, \end{gathered}$$

$$\begin{gathered} Z(\phi ,t)=\int \mathcal{D}\phi (t)\delta (\phi (t)-\phi )\exp [-\frac{\tau }{2}{\int }_0^{t}d{t}^{\prime }{\left(\frac{d\phi }{dt}\right)}^2 \\ -{\int }_0^{t}d{t}^{\prime }\left(r_0{e}^{-\phi (t^{\prime })}+\phi (t^{\prime })\sum _{i=1}^n\delta (t^{\prime }-t_i)\right)]. \end{gathered}$$

$$P(\phi ,t)=\frac{Z(\phi ,t)}{Z(t)},\text{with}Z(t)={\int }_{-\infty }^{\infty }d{\phi }^{\prime }Z({\phi }^{\prime },t).$$

$$\frac{\partial P(\phi ,t)}{\partial t}=-r_0(e^{-\phi }-⟨e^{-\phi }⟩)P(\phi ,t)+\frac{1}{2\tau }\frac{{\partial }^{2}P}{\partial {\phi }^2},$$

$$P(\phi ,t_i^{+})=\frac{e^{-\phi }P(\phi ,t_i^{-})}{⟨e^{-\phi }⟩},$$

$$\frac{d\hat{\phi }}{dt}$$

$$\frac{d{\sigma }^2}{dt}$$

$$d\left(\sum _{i=1}^n\delta (t-t_i)\right)\approx r_0{e}^{-{\phi }^{\ast }}dt+(r_0{e}^{-{\phi }^{\ast }}{)}^{1/2}d{W}^{\prime },$$

$$d\hat{\phi }\approx (r_0{e}^{-\hat{\phi }}/\tau {)}^{1/2}({\phi }^{\ast }-\hat{\phi })+{\tau }^{-1/2}d{W}^{\prime },$$

$$dϵ=-(r/\tau {)}^{1/2}ϵdt+{\tau }^{-1/2}d{W}^{\prime }-({\tau }^{\ast }{)}^{-1/2}dW.$$

$$⟨{ϵ}^2⟩=\frac{1}{2\sqrt{r}}\left(\frac{1}{\sqrt{\tau }}+\frac{\sqrt{\tau }}{{\tau }^{\ast }}\right),$$

$$\frac{⟨(\hat{c}-c^{\ast }{)}^2⟩}{c^2}\approx ⟨{ϵ}^2⟩=\frac{1}{\sqrt{r\tau }}=\frac{1}{\sqrt{4Dac\tau }}.$$

$$\frac{\delta c^2}{c^2}=\frac{1}{4DacT}$$

$$\frac{d\hat{c}}{dt}=\sqrt{4Da\hat{c}/\tau }\left(\frac{1}{4Da}\sum _{i=1}^n\delta (t-t_i)-\hat{c}\right).$$

$$\Delta [A^{\ast }]=k_A^{+}[A][B^{\ast }],$$

$$\frac{d[A^{\ast }]}{dt}=-k_A^{-}[B^{\ast }][A^{\ast }],$$

$$\frac{d[B^{\ast }]}{dt}=k_B^{+}[B][A^{\ast }]-k_B^{-}[B^{\ast }{]}^2,$$

$$\begin{gathered} \frac{d[A^{\ast }]}{dt}=\alpha \sqrt{[A^{\ast }]} \\ \left(\beta \sum _{i=1}\delta (t-t_i)-[A^{\ast }]\right). \end{gathered}$$

$$d\phi (t)={\tau }^{-1/2}dW.$$

$$P_{\mathrm{prior}}(\{\phi (t)\})=\frac{1}{Z_{\mathrm{prior}}}\exp \left[-\frac{\tau }{2}{\int }_0^{T}dt{\left(\frac{d\phi }{dt}\right)}^2\right],$$

$$P(\{\phi (t)\})=\frac{P(t_1,\dots ,t_n|\{\phi (t)\})P_{\mathrm{prior}}(\{\phi (t)\})}{P(t_1,\dots ,t_n)}=\frac{1}{Z}\exp \left\{-{\int }_0^{T}dt\left[\frac{\tau }{2}{\left(\frac{d\phi }{dt}\right)}^2+r_0{e}^{-\phi (t)}\right]-\sum _{i=1}^n\phi (t_i)\right\},$$

$$Z(\phi ,t)=\int \mathcal{D}\phi \delta (\phi (t)-\phi )\exp \left\{{\int }_0^{t}d{t}^{\prime }\left[-\frac{\tau }{2}{\left(\frac{d\phi }{dt}\right)}^2-r_0{e}^{-\phi (t^{\prime })}-\phi (t^{\prime })\sum _{i=1}^n\delta (t^{\prime }-t_i)\right]\right\}.$$

$$P(\phi ,t)=\frac{Z(\phi ,t)}{Z(t)},Z(t)=\int d{\phi }^{\prime }Z({\phi }^{\prime },t),$$

$$Z(\phi ,t+dt)=\int d{\phi }^{\prime }Z({\phi }^{\prime },t)\exp \left[-\frac{\tau (\phi -{\phi }^{\prime }{)}^2}{2dt}-dt{r}_0{e}^{-\phi }\right].$$

$$P(\phi ,t+dt)=\frac{1}{z(t)}\int d{\phi }^{\prime }P({\phi }^{\prime },t)\exp \left[-\frac{\tau (\phi -{\phi }^{\prime }{)}^2}{2dt}-dt{r}_0{e}^{-\phi }\right],$$

$$z(t)=\int d\phi d{\phi }^{\prime }P({\phi }^{\prime },t)\exp \left[-\frac{\tau (\phi -{\phi }^{\prime }{)}^2}{2dt}-dt{r}_0{e}^{-\phi }\right]$$

$$\frac{\partial P(\phi ,t)}{\partial t}=-r_0(e^{-\phi }-⟨{\mathrm{e}}^{-\phi }⟩)P(\phi ,t)+\frac{1}{2\tau }\frac{{\partial }^{2}P}{\partial {\phi }^2},$$

$$z(t)=(1-dt⟨e^{-\phi }⟩)\sqrt{2\pi dt/\tau }.$$

$$\frac{\tilde{Z}(t+dt)}{\tilde{Z}(t)}=(1-dt⟨e^{-\phi }⟩),\text{or}\frac{d\ln \tilde{Z}(t)}{dt}=-⟨e^{-\phi }⟩.$$

$$P(\phi ,t_i+dt)=\frac{e^{-\phi }P(\phi ,t_i)}{⟨e^{-\phi }⟩},$$

Full text

Cells must respond to extracellular signals to guide their actions in the world. The signals typically come in the form of changing concentrations of various molecular ligands, which are conveyed to the cell through ligand binding to cell surface receptors. A lot of ink has been expended on deriving the fundamental limits to the precision with which a cell can measure the concentrations from the activity of its receptors, constrained by the stochasticity of ligand binding and unbinding Berg and Purcell (1977); Bialek and Setayeshgar (2005); Kaizu et al. (2014); Aquino et al. (2016). In particular, it has become clear that the temporal sequence of binding-unbinding events carries more information about the underlying ligand concentration than just the mean receptor occupancy, typically used in deterministic chemical kinetics models of this problem Endres and Wingreen (2009). In particular, such precise temporal information allows cells to estimate the concentration of a cognate ligand even in a sea of weak spurious ligands Siggia and Vergassola (2013); Lalanne and François (2015); Mora (2015), as well as to estimate concentrations of multiple ligands from fewer receptor types Singh and Nemenman (2017, 2019), and molecular network motifs able to perform such complex estimation exist in the real world, even potentially taking advantage of cross-talk between receptor-ligand pairs Carballo-Pacheco et al. (2019).

Importantly, concentrations of ligands are worth measuring only when they are a priori unknown; or, in other words, if they change with time, allowing for instance cells to adapt their behaviour accordingly and maximize their long-term growth Kussell and Leibler (2005). However, all of the preceding analyses have focused on the regime with a clear time scale separation, where the concentration is constant or constantly changing Mora and Wingreen (2010) during the period over which it is estimated. In this article, we will fill in this gap by calculating the accuracy with which a temporally varying ligand concentration may be estimated from a sequence of binding and unbinding events. This requires making assumptions about the time scale over which significant changes of the concentration are possible. In our formulation, the optimal sensor performs a Bayesian computation, formalized mathematically as a stochastic field theory. Crucially, we show how simple biochemical circuits allow one to perform the relevant complex computations.

Field theory of concentration sensing. We associate to the ligand concentration $c(t)$ a field $\varphi(t)$ through $c(t)=c_{0}e^{-\varphi(t)}$, where $c_{0}$ is an irrelevant reference concentration. Ligand concentration controls the ligand-receptor binding rate $r(t)=4Dac(t)=4Dac_{0}e^{-\varphi(t)}\equiv r_{0}e^{-\varphi(t)}$, where $4Da$ is the diffusion-limited binding rate per molecule of the ligand to its target receptor, modeled as a circle of diameter $a$ on the cell’s surface, and $D$ is the ligand diffusivity. This binding rate can be readily generalized to $N$ receptors by using instead $r(t)=4NDac(t)$. All our results will then hold with this additional $N$ factor. We assume that the concentration follows a geometric random walk, with characteristic time scale $\tau$: $d\varphi=\tau^{-1/2}dW$, with $W$ a Wiener process. This choice is justified by the fact that in many biological contexts, such as bacterial chemotaxis, concentrations may vary over many orders of magnitude.

The probability of the concentration temporal evolution over the time interval $(0,T)$ is given by

[TABLE]

The receptor sees binding events at times $t_{1},t_{2},\ldots,t_{n}$, each occuring with rate $4Dac(t_{i})=r_{0}e^{-\varphi(t_{i})}$. To simplify, let us assume that unbinding is instantaneous (generalization to finite binding times is discussed later). The posterior distribution of the concentration profile then follows Bayes’ rule:

[TABLE]

where $Z$ is a normalization constant independent of $\varphi$. The term $r_{0}e^{-\varphi}dt$ in the integral corresponds the probability of not binding a ligand between $t$ and $t+dt$ (except at times $t_{i}$). The binding events at $t=t_{i}$ are generated by the true temporal trace of ligand concentration, $c^{*}(t)=c_{0}e^{-\varphi^{*}(t)}$. In the following the true trace $\varphi^{*}(t)$ will be distinguished from the field $\varphi$, which refers to our observation-based belief.

The one-dimensional field-theoretic problem (2) is a particular case of Bayesian filtering Chen (2003). When collecting information from binding events, cells do not have access to the future and cannot use the full span $[0,T]$ of observations to infer the concentration at time $t$. Instead, they must infer it solely based on past observation in the interval $[0,t]$, which distinguishes our problem from the mathematically similar inference of a continuous probability density Bialek et al. (1996); Nemenman and Bialek (2002); Kinney (2014, 2015); Chen et al. (2018). This inference can be performed recursively by the rules of Bayesian sequential forecasting, similar to the transfer matrix technique, and also known as the forward algorithm Chen (2003). To do this recursion, we first define:

[TABLE]

Considering past observations during the interval $[0,t]$, the posterior distribution of $\varphi$ at time $t$ reads:

[TABLE]

When considering periods during which no binding event was observed, we can write a recursion for $Z(\varphi,t)$ between $t$ and $t+dt$. Taking the $\delta t\to 0$ limit yields, for $t\neq t_{i}$ (App. A SI ):

[TABLE]

where $\langle\cdot\rangle$ denotes an average over $P(\varphi)$. When a binding even does occur at time $t_{i}$, the posterior distribution is updated using Bayes’ rule:

[TABLE]

where $t_{i}^{\pm}$ refer to the values right before and after the observation. The partition function $Z(t)$ can be similarly calculated (App. A SI ) and could in principle be used to infer the correct timescale $\tau$ by maximizing $P(\tau|\{t_{1},\ldots,t_{N}\})\propto Z$ (App. C SI ).

Gaussian solution. Because of the $P(\varphi)$ dependence in $\langle e^{-\varphi}\rangle$, the equations for the evolution of the posterior probability (5)-(6) are nonlinear. However, assuming a Gaussian Ansatz $P(\varphi,t)=(2\pi\sigma(t)^{2})^{-1/2}\exp[-(\varphi-\hat{\varphi}(t))^{2}/2\sigma(t)^{2}]$, which is accurate in the limit of long measurement times (see below), gives a closed-form solution (App. B SI ), with:

[TABLE]

The maximum a posterior estimator for the concentration is then simply given by $\hat{c}(t)=c_{0}e^{-\hat{\varphi}(t)}$, while $\sigma(t)^{2}$ defines the Bayesian uncertainty on the estimator.

To check the validity of the Gaussian solution, we simulated (5)-(6) numerically, starting from a uniform distribution ($P(\varphi,0)=1/2$ for $\varphi\in[-1,1]$ and [math] otherwise), with $r_{0}\tau=50$ and a true $\varphi^{*}(t)$ starting at $\varphi^{*}(0)=0$. The numerical solution quickly approaches the Gaussian solution given by (7)-(8) starting with $\hat{\varphi}(0)=\langle\varphi\rangle_{t=0}$ and $\sigma(0)^{2}=\mathrm{Var}(\varphi)_{t=0}$. The Kullback-Leibler divergence between the numerical and analytical solutions falls rapidly (Fig. 1A) and the numerical solution approaches the predicted Gaussian very closely (Fig. 1A, inset). Thus, the Gaussian solution provides an excellent approximation.

Error estimate. To study the typical behaviour of (7)-(8), we now assume that the rate of binding events is large compared to the rate of change of the concentration, $4Dac\tau=r\tau\gg 1$. This regime is the biologically relevant one: to sense concentration, cells need to record many binding events over the time scale on which the concentration fluctuates. In that limit the estimator $\hat{\varphi}$ is close to the true value $\varphi^{*}$, and the Bayesian uncertainty $\sigma^{2}$ is small, allowing for two simplifications. First, (8) relaxes over time scale $r(t)^{-1}$ to a quasi-steady state value $\sigma^{2}\approx 1/\sqrt{r_{0}e^{-\hat{\varphi}}\tau}\ll 1$. Second, we can make a small noise approximation for binding events: over some time interval $\Delta t$, with $r^{*}(t)^{-1}\ll\Delta t\ll\tau$, the number of binding events has both mean and variance equal to $r^{*}(t)\Delta t$, allowing us to replace discrete jumps in (7) by:

[TABLE]

where $W^{\prime}$ is a Wiener process. As a result, the estimator $\hat{\varphi}$ tracks the true value $\varphi^{*}$ according to:

[TABLE]

where we have expanded at first order in $\hat{\varphi}-\varphi^{*}$. In the general case, the true field may evolve according to a different characteristic time scale, $\tau^{*}$, than the one assumed by the Bayesian filter, $\tau$, so that $d\varphi^{*}=(\tau^{*})^{-1/2}dW$. The estimation error $\epsilon=\hat{\varphi}-\varphi^{*}$ then evolves according to:

[TABLE]

Intrigingy, the noises $dW^{\prime}$ and $dW$ have very different interpretations, one being due to the random arrival of binding events, and the other to the geometric diffusion of the concentration. Yet they come in the same form in this equation. Relying again on the assumption that $r\tau\gg 1$, we get an estimate of the error:

[TABLE]

which has a minimum as a function of $\tau$, reached for the true value of the characteristic fluctuation time $\tau=\tau^{*}$:

[TABLE]

This error is equal to the Bayesian uncertainty $\sigma^{2}=1/\sqrt{r_{0}\tau e^{-\hat{\varphi}}}\approx 1/\sqrt{4Dac\tau}$ and is consistent with the error found using the saddle-point approximation in the related problem of probability density estimate Bialek et al. (1996).

We checked the validity of our small-noise approximation by comparing the prediction from (12) with the results of a numerical simulation of (7)-(8), in which we averaged the error $\langle(\hat{c}-c^{*})^{2}\rangle$ as a function of $c$ for many realization of the process. The agreement is found to be excellent, and gets better as $r\tau=4Dac\tau$ becomes larger (Fig. 1B).

The error in (13) sets a fundamental physical limit on any concentration sensing device, biological or artificial, in a concentration profile that follows a geometric random walk. This bound is radically different from that obtained by Berg and Purcell for the concentration sensing error by a single receptor integrating binding events over time $T$ Berg and Purcell (1977); Endres and Wingreen (2009):

[TABLE]

(in the limit where binding events are short so that the receptor is always free).

The major difference is that Berg and Purcell, as well as most of the literature on concentration sensing, assume that the sensed concentration does not change with time. Our result can be reconciled with Berg and Purcell by defining an effective measurement time $T\sim\sqrt{\tau/4Dac}$ — the geometric mean between the mean time between binding events and the time scale of variation. This $T$ realizes the optimal tradeoff between the requirement to integrate over many binding events, $T\gg 1/(4Dac)$, but over a relatively constant concentration, $T\ll\tau$ Kolmogoroff (1939).

Plausible biological implementation. Can cells implement the optimal Bayesian filtering scheme and reach the bound set by (13)? To gain intuition, it is useful to rewrite (7)-(8) in term of the concentration estimator $\hat{c}$, in the limit $4Dac\tau\gg 1$ where $\sigma^{2}$ can be eliminated:

[TABLE]

Each binding event should lead to an increment of $\hat{c}$, followed by a continuous, exponential decay, with a rate given by $T^{-1}=\sqrt{4Da\hat{c}/\tau}$.

This scheme can be implemented by a simple biochemical network schematized in Fig. 2A. The concentration readout $\hat{c}_{A}$ may be represented by the “active” (for instance phosphorylated) form $A^{*}$ of a chemical species. Binding events cause the receptor to activate $A$ into $A^{*}$, which gets subsequently deactivated. Both the activation and deactivation of $A$ are catylized by a second chemical species in its active form, $B^{*}$. Thus, upon a binding event, the concentration of active $A^{*}$ is increased by:

[TABLE]

and it decays between binding events according to:

[TABLE]

where $k_{A}^{\pm}$ are biochemical parameters.

To implement (15), the concentration of $B^{*}$ must be controled by the square root of $A^{*}$. This dependence can be achieved by assuming that $B$ is activated into $B^{*}$ through the catalityc activity of $A^{*}$, and that $B^{*}$ gets deactivated cooperatively as a dimer:

[TABLE]

where $k_{B}^{\pm}$ are biochemical reaction rates.

Assuming that the kinetics of $B$ are fast compared to $A$, we obtain $B^{*}=(Bk^{+}_{B}/k^{-}_{B})^{1/2}\sqrt{A^{*}}$ and

[TABLE]

with $\alpha=k_{A}^{-}({[B]k^{+}_{B}/k^{-}_{B}})^{1/2}$ and $\beta=(k_{A}^{+}[A]/k_{A}^{-})$. If $A$ and $B$ are in excess, and thus approximately constant, then this biochemical network exactly implements (15), with $4Da\hat{c}_{A}\equiv k_{A}^{-}[A^{*}]/k_{A}^{+}[A]$, and $\tau=\tau_{\rm net}\equiv 1/(\alpha^{2}\beta)=k^{-}_{B}/(k^{+}_{B}k_{A}^{+}k_{A}^{-}[A][B])$.

Interestingly, the amount of inactive ($\approx$ total) $B$ controls the time scale of concentration fluctuations, and could be tuned through gene regulation to adapt to different speeds of environmental fluctuations. A biochemical network might be able to find the optimal $\tau$ and then adjust $[B]$ accordingly by empirically measuring the fold-change of $r(t)$ (which can be done by biochemical networks, see e.g. Goentoro et al. (2009)) but with a delay, $\langle r(t+\Delta t)/r(t)\rangle=e^{\Delta t/2\tau}$, and then inverting the relationship to extract $\tau$.

We tested the performance of the biochemical network for sensing concentration by simulating (16)-(18) with a fluctuating ligand concentration $c(t)$ with characteristic time scale $\tau$. For concreteness, we set $c^{*}(0)=10$nM, $\tau^{*}=10$s, $k_{A}^{+}[A]=0.01$, $k_{A}^{-}=k_{B}^{+}=k_{B}^{+}=1\mu{\rm M}^{-1}{\rm s}^{-1}$ and $[B]=10\mu$M, so that $\tau_{\rm net}=\tau^{*}$. Fig. 2B shows the network estimate $\hat{c}_{A}(t)$ along with the true value $c^{*}(t)$. The empirical error $\langle(\hat{c}_{A}-c^{*})^{2}\rangle$ as a function of $c^{*}$ averaged over $10^{4}$s (Fig. 2B, inset), again shows an excellent agreement with the theoretical bound $1/\sqrt{4Dac\tau}$.

Discussion. For the sake of clarity our analysis made simplifying assumptions which can be easily relaxed. Our proposed biochemical implementation assumed a constant burst of activity following each binding event, consistent with the optimal estimation strategy. However, in real receptors, stochasticity in the bound time is known to double the variance in the estimate Endres and Wingreen (2009) (App. D SI ). Treating this effect simply adds a factor $\sqrt{2}$ in the noise term of (9) as well as in (13), $\langle\delta c^{2}\rangle/c^{2}\approx 1/\sqrt{2Dac\tau}$. We also ignored periods during which the receptor was bound. During that time the receptor is blind to the external world, and the posterior evolves according to the prior: $\partial_{t}P=(1/2\tau)\partial^{2}_{\varphi}P$, $\partial_{t}\hat{\varphi}=0$ and $\partial_{t}\sigma^{2}=1/\tau$. In our results, these “down times” renormalize the effective observation time by the fraction of time the receptor is free, $p_{\rm free}=(1+4Dacu)^{-1}$, where $u$ is the average bound time, $\langle\delta c^{2}\rangle/c^{2}\approx 1/\sqrt{4Dacp_{\rm free}\tau}$ (App. D SI ). Combining the two effects (stochasticity in bound time and receptor availability) would yield $\langle\delta c^{2}\rangle/c^{2}\approx 1/\sqrt{2Dacp_{\rm free}\tau}$.

The field theory of (2) is mathematically similar to the problem of estimating a density function from a small sample set with a smoothing prior Bialek et al. (1996); Nemenman and Bialek (2002); Kinney (2014); Chen et al. (2018). The main difference lies in the domain of observations. In density estimation the whole function $\{\varphi(t)\}_{t\in[0,T]}$ is infered together on the whole domain of $t$, while sensors can only learn from past observations, i.e. the $t^{\prime}<t$ half-plane. However, our solution can easily be generalized to deal with the entire time domain using the forward-backward algorithm (App. E SI ). Eqs. (5)-(6) and (7)-(8) can be solved both forward (from [math] to $t$) and backward (from $T$ to $t$, with time reversal) in time, giving $P_{\to}(\varphi)$, $\hat{\varphi}_{\to}$, $\sigma^{2}_{\to}$ for the forward solution (the one treated in this article), and $P_{\leftarrow}(\varphi)$, $\hat{\varphi}_{\leftarrow}$, $\sigma^{2}_{\leftarrow}$ for the backward solution. The Bayesian posterior at any given time is then given by $\propto P_{\to}(\varphi)P_{\leftarrow}(\varphi)$, of mean $(\sigma^{2}_{\leftarrow}\hat{\varphi}_{\to}+\sigma^{2}_{\to}\hat{\varphi}_{\leftarrow})/(\sigma^{2}_{\leftarrow}+\sigma^{2}_{\to})$ and variance $\sigma^{2}_{\leftarrow}\sigma^{2}_{\to}/(\sigma^{2}_{\leftarrow}+\sigma^{2}_{\to})$ in the Gaussian approximation. While this situation is not relevant for concentration sensing, our general solution should be applicable to problems of density estimation. The saddle-point approximation usually made in that context Bialek et al. (1996); Nemenman and Bialek (2002); Kinney (2014) is expected to work in the same limit as our Gaussian Ansatz; however, recent work has emphasized the importance of non-Gaussian fluctuations for small datasets Chen et al. (2018).

The biological implementation we propose is speculative. An interesting direction would be to identify square-root or similar control of receptor signaling in real biological systems, and interpret them in terms of optimal Bayesian filtering. Signaling pathways dealing with concentration changes over several orders of magnitude, such as bacterial chemotaxis, typically use adaptation mechanisms to increase the dynamic range of sensing Lazova et al. (2011)—a feature that is absent from our approach as we neglect noise in the signaling output. Combining adaptation design with ideas from Bayesian estimation could help us gain insight into the fundamental bounds and resource allocation tradeoffs that limit biological information processing.

Acknowledgments. The authors would like to thank the Casa Matemática Oaxaca from the Banff International Research Station where this work was initiated. TM was partially supported by Agence National pour la Recherche (ANR) grant No. ANR-17-ERC2-0025-01 “IRREVERSIBLE” and IN by NSF Grants No. PHY-1410978 and IOS-1822677.

Appendix A Field theory for concentration sensing

We first recall the problem outlined in the main text for self-consistency. Receptor binding happens with rate $r(t)=4Dac(t)=r_{0}e^{-\varphi(t)}$. The field $\varphi(t)$ follows a random walk with characteristic time $\tau^{*}$:

[TABLE]

In the following $\varphi^{*}$ will refer to the actual realization of the random walk, with characteristic time $\tau^{*}$, and $\varphi$ will refer to the guess made based on the observation of binding events, while $\tau$ denotes the assumes characteristic time scale.

The probability of the time trace of $\varphi$ in the absence of any observation is given by

[TABLE]

with $Z_{\rm prior}=(2\pi dt/\tau)^{T/2dt}$, where $dt$ is an infinitesimal discretization scale.

During each interval $[t,t+dt]$ without a binding event, the likelihood reads $e^{-dtr(t)}=e^{-dtr_{0}e^{-\varphi(t)}}$. A binding event in interval $[t,t+dt]$ has likelihood $dtr_{0}e^{-\varphi(t)}$. Thus the posterior probability thus reads:

[TABLE]

We define:

[TABLE]

The marginal of $\varphi$ at time $t$, $P(\varphi,t|\{t_{1},\ldots,t_{n^{\prime}}\})$, where $n^{\prime}$ is the last binding event before $t$, is then:

[TABLE]

and $Z=Z(T)$.

The partial partition function $Z(\varphi,t)$ can be computed recursively. Let us start with the case where no binding occurs between $t$ and $t+dt$. Then:

[TABLE]

Equivalently,

[TABLE]

where

[TABLE]

is a normalization constant, which can be used to calculate $Z(t)$ recursively: $Z(t+dt)=Z(t)z(t)$. Let us now take the limit $dt\to 0$. The Gaussian integral becomes infinitely peaked. Expanding $P(\varphi^{\prime},t)=P(\varphi,t)+\partial P(\varphi,t)/\partial\varphi(\varphi^{\prime}-\varphi)+\partial^{2}P(\varphi,t)/\partial\varphi^{2}(\varphi^{\prime}-\varphi)^{2}/2$, we obtain:

[TABLE]

where $\langle f(\varphi)\rangle=\int d\varphi\,P(\varphi,t)f(\varphi)$, and

[TABLE]

Defining $Z(t)=(2\pi dt/\tau)^{t/2dt}\tilde{Z}(t)$, we get:

[TABLE]

Now assume that there is a binding event between $t=t_{i}$ and $t_{i}+dt$. Then $P(\varphi,t)$ is discontinuous at $t_{i}$ and terms of order $dt$ can be ignored:

[TABLE]

and $\tilde{Z}(t_{i}+dt)/\tilde{Z}(t_{i})=\langle e^{-\varphi}\rangle$.

In summary, the evolution equation for $P$ reads:

[TABLE]

and the partition function is given by

[TABLE]

Appendix B Gaussian approximation

Because of the term $\langle e^{-\varphi}\rangle$, (32) is non-linear and cannot be solved analytically. However, if we assume that $P(\varphi,t)$ is Gaussian,

[TABLE]

then closed equations can be obtained for the mean $\hat{\varphi}(t)$ and variance $\sigma^{2}(t)$:

[TABLE]

where we have used the Gaussian integral rules: $\langle e^{-\varphi}\rangle=e^{-\hat{\varphi}+\sigma^{2}/2}$, $\langle\varphi e^{-\varphi}\rangle=(\hat{\varphi}-\sigma^{2})e^{-\hat{\varphi}+\sigma^{2}/2}$, $\langle\varphi^{2}e^{-\varphi}\rangle=(\sigma^{2}+(\hat{\varphi}-\sigma^{2})^{2})e^{-\hat{\varphi}+\sigma^{2}/2}$, and $\langle\varphi^{2}\rangle=\hat{\varphi}^{2}+\sigma^{2}$, and we have used integration by parts to calculate the integrals.

Appendix C Partition function and time scale inference

The most likely timescale $\tau$ can be inferred from the observations as well by using Bayes’s rule again:

[TABLE]

where we have used $Z_{\rm prior}(\tau)=(2\pi dt/\tau)^{T/2dt}$.

We can calculate $\tilde{Z}$ from the Gaussian approximation (37):

[TABLE]

This expression looks like the log-likelihood of a sequence of binding events, up to the $\sigma^{2}$ corretions. Bear in mind that $\hat{r}(t)$ is the estimated rate, not the true one $r^{*}(t)$. We have $\hat{r}(t)=r^{*}(t)e^{-\epsilon(t)}$, where we $\epsilon=\hat{\varphi}-\varphi^{*}$. Expanding in $\epsilon$ and $\sigma^{2}$, we obtain:

[TABLE]

Both $\epsilon(t)$ and $\sum_{i}\delta(t-t_{i})-r^{*}$ are stochastic processes of mean 0. They are also uncorrelated with each other, so that the third term is sub-linear in $T$. The last term, which scales with $T$, thus dominates the $\tau$-dependent part of the likelihood. It is maximized for minimum mean squared error, that is for $\tau=\tau^{*}$, as shown in the main text. At large $T$, the $\tilde{Z}(\tau)$ term exponentially dominates the prior $P_{\rm prior}(\tau)$, so that $P(\tau)$ is peaked around the maximum of $\tilde{Z}(\tau)$.

Eq. (39) is an example of the usual bias-variance tradeoff (where “bias” refers to errors made from overfitting the data, and “variance” to errors due to limited data). At small $\tau$, $\hat{r}$ changes rapidly, jumping when a new binding happens, and then rapidly decreases. Thus the term $\sum\ln\hat{r}(t_{i})$ increases, indicating increase in the goodness of fit. At the same time $\sigma^{2}(t)$ increases, so that at small $\tau$ the integral term in (39) becomes large and negative, exploding exponentially for $\tau\to 0$. In contrast, for $\tau\to\infty$, $\hat{r}(t_{i})=n/T$, and $\sigma^{2}\to 0$, so that now the goodness of fit is small.

Overall, there’s an optimal $\tau$ that maximizes $\tilde{Z}$. The three terms in (39) parallel the three terms (fluctuation determinant, goodness of fit, and the kinetic term) in the field-theoretic formulation of continuous probability density estimation from samples Bialek et al. (1996); Nemenman and Bialek (2002), and thus we expect that maximizing $\tilde{Z}$ will result in an optimal $\tau$ not only when the true concentration undergoes a geometric random walk, but also when it undergoes various anomalous walks Nemenman and Bialek (2002).

Appendix D Bound time

We have so far neglected the time the receptor remained bound the ligand Let us denote by $t_{i,\rm off}$ the unbinding time following the binding time at $t_{i}$. When the receptor is bound, $t_{i}<t\leq t_{i,\rm off}$, no information can be obtained from the environment, and the evolution equation for the posterior simply follows the diffusion law:

[TABLE]

The rest of the time, $t_{i,\rm off}<t\leq t_{i}$, (32) holds. Similarly, in the Gaussian approximation, we get

[TABLE]

for bound receptors, $t_{i}<t\leq t_{i,\rm off}$, and (35)-(36) for unbound receptors. In the limit where binding and unbinding events are frequent compared to $\tau$, $r\tau\gg 1$, we have:

[TABLE]

where

[TABLE]

where $u=\langle t_{i,\rm off}-t_{i}\rangle$ is the average bound time, and $r(t)^{-1}=\langle t_{i}-t_{i-1,\rm off}\rangle$ the average unbound time. The uncertainty then reads:

[TABLE]

The time the receptor remains bound has another impact on the ability to sense concentration. In the biochemical scheme proposed in the main text, each binding event causes a fixed burst of activity $\delta(t-t_{i})$, regardless of the bound time. In the simplest receptors however, signaling occurs during the time the receptor is bound, which is itself stochastic. We can model this by replacing the Dirac delta by a random burst of activity, $b_{i}\delta(t-t_{i})$, with $b_{i}$ proportional to the bound time, $b_{i}=(t_{i,\rm off}-t_{i})/u$, so that $\langle b_{i}\rangle=1$ and $\mathrm{Var}(b_{i})=1$, since the bound time is distributed exponentially according to $(1/u)e^{-(t_{i,\rm off}-t_{i})/u}$. More generally we can consider $\langle b_{i}\rangle=1$ and $\mathrm{Var}(b_{i})=CV$, where $0\leq CV\leq$ denotes the coefficient of variation. The general base $\langle b_{i}\rangle$ can be renormalized away into the biochemical parameters. The special case $CV=0$ gives back the results of the main text. When $CV>0$, the variance of $\int_{t_{0}}^{t_{0}+\Delta t}\sum_{i}b_{i}\delta(t-t_{i})$ over an interval of duraction $\Delta t$ instead reads:

[TABLE]

where $n^{\prime}$ is the number of binding events in the interval. As the result, the noise $dW^{\prime}$ in the main text gains a factor $\sqrt{1+CV}$, and the error becomes:

[TABLE]

and minimal error reached for $\tau=(1+CV)\tau^{*}$:

[TABLE]

Taking into account both receptor occupancy and stochasticity in bound times finally yields:

[TABLE]

or, in the case of complete stochastic unbinding, $CV=1$:

[TABLE]

Appendix E Beyond concentration sensing – using the future

Information about future binding events can be exploited by using the backward equation for $P_{\leftarrow}=P(\varphi,t|\{t_{n^{\prime}+1},\ldots,t_{n}\})$:

[TABLE]

where $m$ is the last binding event before $t$. We denote by $P_{\to}=P(\varphi,t|\{t_{1},\ldots,t_{n^{\prime}}\})$ the solution of the forward equation (32) discussed before. The distribution of $\varphi$ at time $t$ is then given by:

[TABLE]

where we have used the fact that the past and future where conditionally independent given $\varphi$ at time $t$, since the process is Markovian, and a uniform prior. We thus have:

[TABLE]

Using the Gaussian solution, and denoting $\hat{\varphi}_{\to},\sigma^{2}_{\to}$ the parameters of the forward solution, and $\hat{\varphi}_{\leftarrow},\sigma^{2}_{\leftarrow}$ those of the backward solution, we obtain:

[TABLE]

with

[TABLE]

These formulas could be used in estimates of density $r(t)$ from sparse observations $t_{i}$, where $r(t)$ is interpreted as a density of events, and $t$ is the variable whose density we want to infer.