Bayesian gradient sensing in the presence of rotational diffusion

Many motile biological cells navigate in concentration gradients of signaling molecules in a process termed chemotaxis [1–6]. At cellular scales, the stochastic binding of signaling molecules results in molecular shot noise and renders concentration measurements inherently noisy [7]. This sensing noise imposes physical limits on the precision of chemotaxis [7–13]. Experimental work suggests that biological cells indeed operate at these limits in shallow concentration gradients [14–18]. Temporal averaging over extended measurement intervals is a common strategy to reduce sensing noise [7, 19–21]. Yet, temporal averaging may be limited in time-varying environments [18, 22–25], or more directly by the stochastic motion of chemotactic agents themselves [26]. How should a chemotactic agent integrate previous and more recent measurements to most accurately estimate the relative direction of a concentration gradient if this direction changes stochastically in time?

Classical work by Berg and Purcell showed that the relative accuracy ⟨(δc/c)²⟩ for sensing absolute concentration c in the presence of measurement noise scales as 1/T with total measurement time T [7]. Intriguingly, Mora et al recently showed that in time-varying environments, this accuracy scales instead as (τr)^−1/2, where τ is the time-scale on which the concentration changes, and r is a rate of information gain [25]. Correspondingly, the optimal effective measurement time for temporal averaging reads T_eff = (τ/r)^1/2 [25]. To the best of our knowledge, it is not known if similar statements hold also for the nonlinear problem of estimating the direction of a concentration gradient. Here, we develop a theory of optimal estimation of gradient direction in the presence of rotational diffusion with rotational diffusion coefficient D, and hence rotational diffusion time τ = (2D)⁻¹. We derive results that are analogous to Mora et al, provided r is now interpreted as a signal-to-noise ratio of gradient sensing normalized by unit time.

Review of previous results. Previous work suggested that bacteria use Kalman filters to track a time-dependent concentration signal [27], providing an optimal weighting of past and recent measurements, see also [26]. Yet, Kalman filters address linear problems [28], while sensing a direction is a nonlinear problem of circular statistics, which prompts Bayesian updating of an angular likelihood distribution at each time step [29]. Bayesian estimation had been successfully applied, e.g. for monitoring time-varying environments with two states [30], estimating absolute concentration [25, 31–33], or temporal changes thereof [24, 34]. In contrast, previous theories of Bayesian sensing of direction [11, 12, 15, 35], or target position [36], did not account for motility noise as we do. In the engineering literature, estimation algorithms known as 'bearing tracking' do account for motility noise [37, 38], yet lack an analytical theory on the theoretical limits of gradient sensing.

Structure of manuscript. Our manuscript is structured as follows: we first introduce our minimal model of gradient sensing and formulate sequential Bayesian estimation, consisting of alternating update and prediction steps for a likelihood distribution of relative gradient direction. We use a suitably discretized time, and the calculus of von-Mises distributions; where von-Mises distributions can be considered the analogue of the normal distribution for direction angles. As our main result, we derive an analytic result for the ultimate precision of gradient estimation, see equation (19). We show that the estimated precision of an individual agent faithfully reflects the true accuracy of estimation errors in an ensemble of identical agents, as expected for optimal Bayesian estimation. Estimation becomes sub-optimal if the agents do not know their rotational diffusion coefficient exactly. If the agents, however, do not account for their rotational diffusion at all, gradient sensing would fail in a characteristic manner. We emphasize that our minimal model is representative of a larger class of problems, e.g. agents exploring time-dependent concentration gradients that likewise change stochastically on a characteristic time-scale τ.

We consider a chemotactic agent that seeks to estimate the direction of an external concentration gradient by detecting binding events of signaling molecules that diffuse to its surface.

As a minimal model, we consider a stationary agent that is positioned at the origin of a two-dimensional plane. The agent is disklike with radius a. To leading order, the external concentration gradient will be linear

$$\begin{equation}c\left(\mathbf{r}\right)={c}_{0}+\nabla c\cdot \mathbf{r}={c}_{0}\left(1+\alpha {\;\mathbf{e}}_{c}\cdot \mathbf{r}/a\right).\end{equation} \tag{ 1 }$$

Here, c₀ denotes a base concentration, α = |∇c|a/c₀ a normalized strength of the concentration gradient, and e_c = ∇c/|∇c| is the unit vector pointing in the direction of the gradient. We introduce a material frame of the agent with orthonormal unit vectors h₁ and h₂, see figure 1(A). The agent seeks to estimate the direction angle ψ* of the gradient relative to the material frame, such that e_c = cos ψ* h₁ + sin ψ* h₂. The circumference r(φ) of the agent is parameterized by r(φ) = a  cos φ  h₁ + a  sin φ  h₂. The rate Λ(φ) at which molecules are detected at the circumference of the agent at position r(φ) is proportional to the local concentration

$$\begin{equation}{\Lambda}\left(\varphi \right)={\left(2\pi \right)}^{-1}\lambda c\left[\mathbf{r}\left(\varphi \right)\right],\end{equation} \tag{ 2 }$$

where λ is a proportionality constant that depends on the diffusion constant of the signaling molecules.

Zoom In Zoom Out Reset image size

The sequence of binding events at respective times θ_j and polar direction angles φ_j can be conveniently encoded in a complex signal

$$\begin{equation}s\left(t\right)=\sum\limits _{j}\delta \left(t-{\theta }_{j}\right){\mathrm{e}}^{\mathrm{i}{\varphi }_{j}}.\end{equation} \tag{ 3 }$$

By definition, the expectation value of this signal equals the Fourier transform of the angle-dependent binding rate, ⟨s(t)⟩ = ∮dφ Λ(φ) exp(iφ). Generally, it is much harder to estimate a gradient direction than absolute concentration c₀. Correspondingly, we focus on the problem of estimating the relative gradient direction angle ψ* and assume that the agent a priori knows c₀ as well as gradient strength α.

We introduce the signal-to-noise ratio of gradient sensing for a given time interval Δt, following [39]

$$\begin{equation}\mathrm{S}\mathrm{N}\mathrm{R}={\alpha }^{2}\lambda {c}_{0}\;\;{\Delta}t/2.\end{equation} \tag{ 4 }$$

This signal-to-noise ratio scales linearly with Δt and defines an information rate with units of an inverse time

$$\begin{equation}r=\frac{\mathrm{S}\mathrm{N}\mathrm{R}}{{\Delta}t}={\alpha }^{2}\lambda {c}_{0}/2.\end{equation} \tag{ 5 }$$

The information rate defines a characteristic time scale of the problem as 1/r. On times shorter than 1/r, the agent cannot gather much useful information on the gradient direction.

We now turn to a time-dependent problem and consider rotational diffusion of the agent with rotational diffusion coefficient D = D_rot. The rotational diffusion coefficient defines a second time-scale τ = 1/(2D). We will show that these two time-scales, 1/r and τ, together determine the ultimate precision of gradient sensing. In the following, we consider the limit of weak rotational diffusion, D ≪ r, which is the biologically relevant case. In the opposite limit, D ≳ r, the relative gradient direction would change so fast that the agent could simply not track it.

To make analytical progress, we assume a discrete time dynamics with equidistant time points t_j = jΔt. We choose the interval Δt of time discretization such that

$$\begin{equation}\frac{1}{r}\ll {\Delta}t\ll {T}_{\mathrm{e}\mathrm{ff}},\end{equation} \tag{ 6 }$$

where

$$\begin{equation}{T}_{\mathrm{e}\mathrm{ff}}=\sqrt{\tau /r}.\end{equation} \tag{ 7 }$$

Here, we introduced a time-scale T_eff as the geometric mean of 1/r and τ; the time-scale T_eff will later characterize the time-scale of convergence of the estimated precision of gradient sensing to its limit value (and, in fact, represent an optimal effective measurement time). Thus, the second inequality in equation (6) ensures that the time discretization is finer than the dynamics of the problem. Our choice of Δt formally corresponds to a large signal-to-noise ratio SNR = rΔt ≫ 1 for a time interval of duration Δt. At the same time, our choice ensures that the effect of rotational diffusion is small compared to sensing noise during a time interval of duration Δt, i.e. DΔt ≪ 1/SNR. Note that for a time step that is too coarse, alternating update and prediction steps would cause the estimated precision to oscillate; a small time step ensures that the amplitude of these oscillations is small.

At the times t_j , the agent shall undergo discrete rotational diffusion events, see figure 1(B). Correspondingly, the relative gradient direction angle ψ* will be constant in each time interval ${\psi }^{{\ast}}\left(t\right)={\psi }_{j}^{{\ast}}$ for $t\in {\mathcal{T}}_{j}$ with ${\mathcal{T}}_{j}=\left({t}_{j-1},{t}_{j}\right)$, but change by a discrete amount ${{\Delta}}_{j}={\psi }_{j+1}^{{\ast}}-{\psi }_{j}^{{\ast}}$ at t_j . The stochastic increments Δ_j are uncorrelated and satisfy ⟨Δ_j Δ_k ⟩ = 2DΔtδ_jk . While estimating the direction of the concentration gradient, the agent will account for its rotational diffusion, assuming that its rotational diffusion coefficient is $\hat{D}$.

During the time intervals ${\mathcal{T}}_{j}=\left({t}_{j-1},{t}_{j}\right)$, the agent does not rotate and integrates the stochastic signal s(t) for a time span Δt

$$\begin{equation}{\mathcal{M}}_{j}={\int }_{{t}_{j-1}}^{{t}_{j}} \mathrm{d}t\enspace s\left(t\right).\end{equation} \tag{ 8 }$$

The angle ${\psi }_{j}=\mathrm{arg}\enspace {\mathcal{M}}_{j}$ of this complex number ${\mathcal{M}}_{j}$ represents a noisy measurement of the true current gradient direction angle ${\psi }_{j}^{{\ast}}$.

3.1. Statistics of single measurements

What is the statistics of the measurements ${\mathcal{M}}_{j}$? In total, the agent will detect on average ν = ∫dφ Λ(φ)Δt = λc₀Δt molecules during a time interval of duration Δt. A short calculation yields for the mean and variance of ${\mathcal{M}}_{j}$

$$\begin{equation}\langle {\mathcal{M}}_{j}\rangle =\frac{\alpha \nu }{2}\enspace {\mathrm{e}}^{\mathrm{i}{\psi }_{j}^{{\ast}}}\end{equation} \tag{ 9 }$$

$$\begin{equation}\langle \vert {\mathcal{M}}_{j}{\vert }^{2}\rangle -\vert \langle {\mathcal{M}}_{j}\rangle {\vert }^{2}=\nu .\end{equation} \tag{ 10 }$$

For this calculation, we considered all binding events in the time interval ${\mathcal{T}}_{j}$, enumerated by an index k with respective angles φ_k , where the total number ν_j of these events is a Poisson random variable with mean ν and variance ν. We then find $\langle {\mathcal{M}}_{j}\rangle =\langle {\sum }_{k}\mathrm{exp}\enspace i{\varphi }_{k}\rangle =\int \mathrm{d}\varphi \enspace {\Lambda}\left(\varphi \right)\mathrm{exp}\left(i\varphi \right){\Delta}t$, and $\langle \vert {\mathcal{M}}_{j}{\vert }^{2}\rangle =\langle {\sum }_{k\ne l}\mathrm{exp}\enspace i\left({\varphi }_{k}-{\varphi }_{l}\right)\rangle +\langle {\sum }_{k=l}\mathrm{exp}\enspace i\left({\varphi }_{k}-{\varphi }_{l}\right)\rangle ={\left(\alpha /2\right)}^{2}\langle {\nu }_{j}\left({\nu }_{j}-1\right)\rangle +\langle {\nu }_{j}\rangle ={\left(\alpha \nu /2\right)}^{2}+\nu$, where we used (Δt/ν)² ∫∫dφ dφ' Λ(φ)Λ(φ') exp i(φ − φ') = (α/2)². Similarly, $\langle {\mathcal{M}}_{j}^{2}\rangle ={\langle {\mathcal{M}}_{j}\rangle }^{2}$, which implies that fluctuations of ${\mathcal{M}}_{j}$ in the complex plane are isotropic [40].

Our choice of Δt in equation (6) implies in particular ν ≫ 1. In this limit of large molecule counts, we can approximate the probability distribution of the measurements ${\mathcal{M}}_{j}$ by a normal distribution with the same mean $\mu =\langle {\mathcal{M}}_{j}\rangle$ and the same variance Σ = ν/2 of the real and imaginary part of ${\mathcal{M}}_{j}$ (diffusion approximation)

$$\begin{align}\hfill p\left({\mathcal{M}}_{j}\vert {\psi }_{j}^{{\ast}}\right)& ={\left(2\pi {\Sigma}\right)}^{-1}\enspace \mathrm{exp}\left(-\frac{\vert {\mathcal{M}}_{j}-\mu {\vert }^{2}}{\nu }\right)\hfill \\ \hfill & \sim \mathrm{exp}\left(\underset{{\zeta }_{j}}{\underbrace{\alpha \vert {\mathcal{M}}_{j}\vert }}\enspace \mathrm{cos}\left({\psi }_{j}-{\psi }_{j}^{{\ast}}\right)\right).\hfill \end{align} \tag{ 11 }$$

In the last step, we simply used the Pythagorean trigonometric identity. The factor of proportionality depends only on $\vert {\mathcal{M}}_{j}\vert$. If we interpret the last expression as a function of ${\psi }_{j}^{{\ast}}$, it provides a first example of a von-Mises distribution with so-called measure-of-concentration ζ_j . von-Mises distributions can be regarded as a generalization of the normal distribution for angular variables and are further discussed in section 3.2 and appendix C. The expectation value ⟨ζ_j ⟩ of the measure-of-concentration ζ_j for a single measurement is just the signal-to-noise ratio defined in equation (4)

$$\begin{equation}\langle {\zeta }_{j}\rangle \approx \mathrm{S}\mathrm{N}\mathrm{R}.\end{equation} \tag{ 12 }$$

A cursory explanation for equation (12) reads $\langle {\zeta }_{j}\rangle =\alpha \langle \vert {\mathcal{M}}_{j}\vert \rangle \;\;\left({\ast}\right){\approx }\;\;\alpha \vert \langle {\mathcal{M}}_{j}\rangle \vert ={\alpha }^{2}\nu /2=\mathrm{S}\mathrm{N}\mathrm{R}$, where we used equation (6) in step (*); a more rigorous derivation can be found in appendix B.

3.2. Estimation of directions using von-Mises distributions

The agent uses the sequence of measurements ${\mathcal{M}}_{j}$, j = 1, ..., m, to compute a likelihood distribution ${\mathcal{L}}_{m}\left(\psi \right)$ for the unknown gradient direction angle ${\psi }_{m}^{{\ast}}$ at time t_m . In each time step, the agent performs a prediction step (accounting for rotational diffusion), and an update step (accounting for the new measurement ${\mathcal{M}}_{m}$), see figure 1(C). For the prediction step, the agents starts with a prior ${\mathcal{L}}_{m-1}\left(\psi \right)$ that incorporates all previous measurements, ${\mathcal{M}}_{j}$ with j ⩽ m − 1, and computes from it a modified likelihood distribution ${\mathcal{L}}_{m-1}^{\prime }\left(\psi \right)$ that accounts for the loss of knowledge due to the rotational diffusion event at t_m−1. In the update step, the agent computes the new likelihood distribution ${\mathcal{L}}_{m}\left(\psi \right)$ using Bayes' rule, incorporating the new measurement ${\mathcal{M}}_{m}$, with the predicted likelihood ${\mathcal{L}}_{m-1}^{\prime }\left(\psi \right)$ as prior.

A fundamental difficulty of this estimating problem is that it concerns an angular variable, i.e. all likelihood distributions must be 2π-periodic. This sets our problem apart from other problems previously considered in the literature [25, 30].

We will approximate the different likelihood distributions $\mathcal{L}\left(\psi \right)$ by von-Mises distributions

$$\begin{equation}\mathcal{VM}\left(\psi ;\mu ,\kappa \right)\sim \mathrm{exp}\left[\kappa \enspace \mathrm{cos}\left(\psi -\mu \right)\right],\end{equation} \tag{ 13 }$$

where μ denotes the mean and κ is the so-called measure-of-concentration. The von-Mises distribution can be regarded as a generalization of the normal distribution for angular variables, where κ plays a role similar to an inverse variance. The von-Mises distribution is the maximum-entropy distribution p(ψ) that maximizes information entropy ∫dψ p(ψ)ln p(ψ) if the first complex moment ⟨exp(iψ)⟩ is specified. Importantly, the product of two von-Mises distributions is again a von-Mises distribution (up to normalization), $\mathcal{VM}\left(\psi ;{\mu }_{1},{\kappa }_{1}\right)\mathcal{VM}\left(\psi ;{\mu }_{2},{\kappa }_{2}\right)\sim \mathcal{VM}\left(\psi ;\mu ,\kappa \right)$. Here, the new mean μ and new measure-of-concentration κ satisfy κ  exp(iμ) = κ₁ exp(iμ₁) + κ₂ exp(iμ₂), i.e. the new mean is a weighted circular average of the old mean values. If the old mean values are close, we have κ ≈ κ₁ + κ₂, which generalizes an analogous, exact statement for normal distributions. Likewise, μ ≈ (κ₁/κ) μ₁ + (κ₂/κ) μ₂. The convolution of two von-Mises distributions is only approximately a von-Mises distribution, $\mathcal{VM}\left(\psi ;{\mu }_{1},{\kappa }_{1}\right){\ast}\mathcal{VM}\left(\psi ;{\mu }_{2},{\kappa }_{2}\right)\approx \mathcal{VM}\left(\psi ;\mu ,\kappa \right)$. Here, the new mean is the sum of the old mean values, μ = μ₁ + μ₂, and the new measure-of-concentration is approximately the harmonic sum of the old measures-of-concentration, ${\kappa }^{-1}\approx {\kappa }_{1}^{-1}+{\kappa }_{2}^{-2}$. The sharper the distributions are (κ₁, κ₂ ≫ 1), the more accurate this approximation for convolutions becomes. Appendix C provides further information on von-Mises distributions.

We will use products and convolutions of von-Mises distributions in the update and the prediction step, respectively. The mean values of these von-Mises distribution represent the current maximum-likelihood estimates ${\hat{\psi }}_{j}$ of the gradient direction angle.

3.3. Bayesian updating

We assume that just before the updating step at time t_m , the agent starts with a Bayesian prior ${\mathcal{L}}_{m-1}^{\prime }\left(\psi \right)$ for the unknown gradient direction angle ψ*, given in terms of a von-Mises distribution with mean ${\hat{\psi }}_{m-1}$ and measure-of-concentration κ'_m−1

$$\begin{equation}{\mathcal{L}}_{m-1}^{\prime }\left(\psi \right)=\mathcal{VM}\left(\psi ;{\hat{\psi }}_{m-1},{\kappa }_{m-1}^{\prime }\right).\end{equation} \tag{ 14 }$$

In light of the new measurement ${\mathcal{M}}_{m}$, the agent updates this prior according to Bayes' rule

$$\begin{equation}\text{update}\;\text{step:}\;\;\;\;{\mathcal{L}}_{m}\left(\psi \right)=\frac{P\left({\mathcal{M}}_{m}\vert \psi \right)}{P\left({\mathcal{M}}_{m}\vert {\mathcal{M}}_{1:m-1}\right)}{\mathcal{L}}_{m-1}^{\prime }\left(\psi \right).\end{equation} \tag{ 15 }$$

Here, $P\left({\mathcal{M}}_{m}\vert \psi \right)$ is the probability to obtain the specific measurement ${\mathcal{M}}_{m}$ if the true gradient angle were ψ, i.e. the measurement model. Similarly, the normalization factor $P\left({\mathcal{M}}_{m}\vert {\mathcal{M}}_{1:m-1}\right)$ is the probability to obtain this specific measurement ${\mathcal{M}}_{m}$ if the true gradient direction were distributed according to the Bayesian prior ${\mathcal{L}}_{m-1}^{\prime }\left(\psi \right)$; explicitly, $P\left({\mathcal{M}}_{m}\vert {\mathcal{M}}_{1:m-1}\right)=\oint \mathrm{d}\psi \enspace P\left(\mathcal{M}\vert \psi \right){\mathcal{L}}_{m-1}^{\prime }\left(\psi \right)$.

The right-hand side of equation (15) represents a product of von-Mises distributions, which yields again a von-Mises distribution, whose measure-of-concentration shall be called κ_m . Our choice of time step Δt in equation (6) implies that both the mean ${\mu }_{1}={\psi }_{m}=\mathrm{arg}\enspace {\mathcal{M}}_{m}$ of the first factor $P\left({\mathcal{M}}_{m}\vert \psi \right)$ and the mean ${\mu }_{2}={\hat{\psi }}_{m-1}$ of the second factor ${\mathcal{L}}_{m-1}^{\prime }\left(\psi \right)$ are close to the true angle ${\psi }_{j}^{{\ast}}$, provided true and assumed rotational diffusion coefficient are equal, $D=\hat{D}$, or at least of similar magnitude, see below. Consequently, also |μ₁ − μ₂| ≪ 1. The calculus of von-Mises distributions thus yields the iteration rule κ_m ≈ κ'_m−1 + ζ_m , or, using equation (12)

$$\begin{equation}\text{update}\;\text{rule:}\;\;\;\;{\kappa }_{m}\approx {\kappa }_{m-1}^{\prime }+\mathrm{S}\mathrm{N}\mathrm{R}.\end{equation} \tag{ 16 }$$

Fluctuations of ζ_m cause small fluctuations of κ_m not considered here.

3.4. Prediction step: rotational diffusion

For rotational diffusion, the stochastic increments Δ_j of the relative direction angle are distributed according to a wrapped normal distribution with variance parameter 2DΔt. In the limit of weak rotational diffusion considered here, this distribution is well approximated by a von-Mises distribution with zero mean and measure-of-concentration (2DΔt)⁻¹. The agent accounts for this stochastic dynamics of the relative direction angle in a prediction step, assuming its rotational diffusion coefficient to be $\hat{D}$ (which may or may not be equal to the true value D). The prediction step thus corresponds to a convolution of two von-Mises distributions

$$\begin{equation}\text{prediction}\;\text{step:}\;\;\;\;{\mathcal{L}}_{m}^{\prime }\approx {\mathcal{L}}_{m}{\ast}\mathcal{VM}\left({\Delta};0,{\left(2\hat{D}{\Delta}t\right)}^{-1}\right).\end{equation} \tag{ 17 }$$

Equation (17) is a special case of a Chapman–Kolmogorov equation, ${\mathcal{L}}_{m}^{\prime }={\mathcal{L}}_{m}^{\prime }\left(\psi \vert {\mathcal{M}}_{1:m}\right)=\oint P\left({\psi }^{\prime }\vert \psi \right)\mathcal{L}\left(\psi \vert {\mathcal{M}}_{1:m}\right)\mathrm{d}\psi$ with assumed transition probability P(ψ'|ψ) from ψ to ψ' given by $P\left({\psi }^{\prime }\vert \psi \right)\approx \mathcal{VM}\left({\psi }^{\prime };\psi ,{\left(2\hat{D}{\Delta}t\right)}^{-1}\right)$.

By the calculus of von-Mises distributions, we obtain again an iteration rule, now for the measure-of-concentration κ'_m of the likelihood distribution after the prediction step

$$\begin{equation}\text{prediction}\;\text{rule:}\;\;\;\;\frac{1}{{\kappa }_{m}^{\prime }}\approx \frac{1}{{\kappa }_{m}}+2\hat{D}{\Delta}t.\end{equation} \tag{ 18 }$$

The prediction step does not change the maximum-likelihood angle, i.e. ${\hat{\psi }}_{m}^{\prime }={\hat{\psi }}_{m}$, since the diffusion kernel is centered at zero.

3.5. Theoretical limit of gradient sensing in the presence of rotational diffusion

Together, the nested iteration rules equations (16) and (18) define a monotonically increasing sequence of measures-of-concentration κ_m , which converges to a limit value κ_∞ for $D=\hat{D}$. The measure-of-concentration κ_∞ represents a theoretical limit for the estimated precision of gradient sensing at steady-state. We compute κ_∞ as a fixed point of the nested iteration rules equations (16) and (18). Specifically, the limit values κ_∞ = lim_m→∞ κ_m and κ'_∞ = lim_m→∞ κ'_m satisfy κ_∞ = κ'_∞ + SNR and ${\left({\kappa }_{\infty }^{\prime }\right)}^{-1}={\left({\kappa }_{\infty }\right)}^{-1}+2D{\Delta}t$; eliminating κ'_∞ yields ${\kappa }_{\infty }={\left[{\left({\kappa }_{\infty }\right)}^{-1}+2D{\Delta}t\right]}^{-1}+\mathrm{S}\mathrm{N}\mathrm{R}$, or ${\left({\kappa }_{\infty }-\mathrm{S}\mathrm{N}\mathrm{R}\right)}^{-1}={\left({\kappa }_{\infty }\right)}^{-1}+2D{\Delta}t$. Multiplying with the common denominator and rearranging yields the quadratic equation κ_∞(κ_∞ − SNR) = SNR/(2DΔt), hence κ_∞ ≈ [SNR/(2DΔt)]^1/2 since κ_∞ ≫ SNR by equation (6). We thus have

$$\begin{equation}\text{Case}\;D=\hat{D}:{\vert c\vert } {\kappa }_{\infty }=\sqrt{\frac{r}{2D}}=\sqrt{\frac{\lambda }{4\enspace {c}_{0}D}}\vert \nabla c\vert a.\end{equation} \tag{ 19 }$$

This result for κ_∞ highlights the competition between information gain with rate r = SNR/Δt and information loss by rotational diffusion with rotational diffusion coefficient D. We can rewrite equation (19) as κ_∞ = rT_eff, where we introduced the effective measurement time T_eff = (τ/r)^1/2.

This theoretical limit is strict and faithful, i.e. it reflects the distribution of estimation errors in an ensemble of identical agents, as we show next.

3.6. Distribution of estimation errors

The mean ${\hat{\psi }}_{m}$ of ${\mathcal{L}}_{m}\left(\psi \right)$ represents the maximum-likelihood estimate of the agent at time t_m . By the update step, equation (15), which expresses ${\mathcal{L}}_{m}\left(\psi \right)$ as a product of von-Mises distributions [and the rules for products of von-Mises distributions, see appendix C], ${\hat{\psi }}_{m}$ is a weighted circular average of the previous estimate ${\hat{\psi }}_{m-1}$ and the current measurement ψ_m

$$\begin{equation}{\kappa }_{m}\enspace {\mathrm{e}}^{\mathrm{i}{\hat{\psi }}_{m}}={\kappa }_{m-1}^{\prime }\enspace {\mathrm{e}}^{\mathrm{i}{\hat{\psi }}_{m-1}}+\mathrm{S}\mathrm{N}\mathrm{R}\enspace {\mathrm{e}}^{\mathrm{i}{\psi }_{m}}.\end{equation} \tag{ 20 }$$

We consider the distribution p(ɛ_m ) of estimation errors in an ensemble of agents with ${\varepsilon }_{m}={\hat{\psi }}_{m}-{\psi }_{m}^{{\ast}}$. We are interested in the true accuracy of estimation, defined as the circular variance of this distribution, CV[p(ɛ_m )] = 1 − |⟨exp iɛ_m ⟩|. We have $\mathrm{C}\mathrm{V}\left[p\left({\varepsilon }_{m}\right)\right]=1-\langle \mathrm{exp}\enspace i{\varepsilon }_{m}\rangle \approx \langle {\varepsilon }_{m}^{2}/2\rangle$. By multiplying equation (20) with $\mathrm{exp}\left(-{\psi }_{m}^{{\ast}}\right)$, and linearizing the imaginary part, we obtain an iteration rule for the estimation error ɛ_m as an affine weighted sum

$$\begin{equation}{\varepsilon }_{m}\approx \frac{{\kappa }_{m-1}^{\prime }}{{\kappa }_{m}}\underset{\begin{matrix}\hfill \;\text{previous}\;\text{estimate}\hfill \\ \hfill \;+\;\text{motility}\;\text{noise}\hfill \end{matrix}}{\underbrace{\left({\varepsilon }_{m-1}-{{\Delta}}_{m-1}\right)}}+\frac{\mathrm{S}\mathrm{N}\mathrm{R}}{{\kappa }_{m}}\underset{\begin{matrix}\hfill \;\text{measurement}\hfill \\ \hfill \;\text{noise}\hfill \end{matrix}}{\underbrace{\left({\psi }_{m}-{\psi }_{m}^{{\ast}}\right)}}.\end{equation} \tag{ 21 }$$

Recall that ${{\Delta}}_{m-1}={\psi }_{m}^{{\ast}}\hspace{-1pt}-\hspace{-1pt}{\psi }_{m-1}^{{\ast}}$ is the random change of the true gradient direction angle ψ*(t) at time t_m−1.

Case $D=\hat{D}$. We first assume that the agent knows the exact value of its own rotational diffusion coefficient, i.e. $D=\hat{D}$. Equation (21) represents a linear superposition of two independent random angles. In the limit of sufficiently sharp distributions, valid by equation (6), the circular variance of the superposition is a weighted sum of the circular variances of these random angles, with coefficients squared, analogous to the case of normal distributions. We thus conclude from equation (21) an iteration rule for the true accuracy

$$\begin{equation}\mathrm{C}\mathrm{V}\left[p\left({\varepsilon }_{m}\right)\right]\approx {\left(\frac{{\kappa }_{m-1}^{\prime }}{{\kappa }_{m}}\right)}^{2}\left(\mathrm{C}\mathrm{V}\left[p\left({\varepsilon }_{m-1}\right)\right]+D{\Delta}t\right)+{\left(\frac{\mathrm{S}\mathrm{N}\mathrm{R}}{{\kappa }_{m}}\right)}^{2}\mathrm{C}\mathrm{V}\left[p\left({\psi }_{m}-{\psi }_{m}^{{\ast}}\right)\right].\end{equation} \tag{ 22 }$$

Note that $\mathrm{C}\mathrm{V}\left[p\left({\psi }_{m}-{\psi }_{m}^{{\ast}}\right)\right]={\left(2\mathrm{S}\mathrm{N}\mathrm{R}\right)}^{-1}$ is just the sensing noise of a single measurement. For agents that are aware of their rotational diffusion with $\hat{D}=D$, this iteration rule for the true accuracy CV[p(ɛ_m )] in an ensemble of identical agents, equation (22), is equivalent to the combination of update and prediction rule for the estimated precision κ_m of an individual agent, equations (16) and (18). Indeed, if CV[p(ɛ_m−1)] ≈ 1/(2κ_m−1) at t = t_m−1, then CV[p(ɛ_m )] ≈ 1/(2κ_m ) at t = t_m by a short calculation. Specifically, from the initial assumption it follows that CV[p(ɛ_m−1)] + DΔt ≈ 1/(2κ'_m−1), hence the right-hand side of equation (22) becomes ${\left[{\kappa }_{m-1}^{\prime }/\left({\kappa }_{m-1}^{\prime }+\mathrm{S}\mathrm{N}\mathrm{R}\right)\right]}^{2}{\left(2{\kappa }_{m-1}^{\prime }\right)}^{-1}+{\left[\mathrm{S}\mathrm{N}\mathrm{R}/\left({\kappa }_{m-1}^{\prime }+\mathrm{S}\mathrm{N}\mathrm{R}\right)\right]}^{2}{\left(2\mathrm{S}\mathrm{N}\mathrm{R}\right)}^{-1}=\left({\kappa }_{m-1}^{\prime }+\mathrm{S}\mathrm{N}\mathrm{R}\right)/\left[2{\left({\kappa }_{m-1}^{\prime }+\mathrm{S}\mathrm{N}\mathrm{R}\right)}^{2}\right]=1/\left[2\left({\kappa }_{m-1}^{\prime }+\mathrm{S}\mathrm{N}\mathrm{R}\right)\right]=1/\left(2{\kappa }_{m}\right)$. This confirms that Bayesian estimation of gradient direction is faithful, i.e. the estimated precision of an individual agent follows, on average, the true accuracy in an ensemble of agents, as expected for optimal Bayesian estimation. As a corollary, the true accuracy converges to the same limit value as the (circular variance corresponding to the) estimated precision

$$\begin{equation}\underset{m\to \infty }{\mathrm{lim}}\enspace \mathrm{C}\mathrm{V}\left[p\left({\varepsilon }_{m}\right)\right]\approx {\left(2{\kappa }_{\infty }\right)}^{-1}.\end{equation} \tag{ 23 }$$

Figure 2(A) shows simulation results for a case where estimated precision and true accuracy are initially different, but quickly converge to a common limit value. For details on numerical methods, see appendix A.

Zoom In Zoom Out Reset image size

Case $\hat{D}\ne D$. More generally, we can consider agents that assume a value $\hat{D}$ for their rotational diffusion coefficient when performing the prediction step equation (18), while the true rotational diffusion coefficient is D. In this case, the estimated precision κ_m converges to ${\kappa }_{\infty }\approx {\left[\mathrm{S}\mathrm{N}\mathrm{R}/\left(2\hat{D}{\Delta}t\right)\right]}^{1/2}$, i.e. equation (19) with D replaced by $\hat{D}$, while the true accuracy CV[p(ɛ_m )] converges to a limit value CV_∞ different from 1/(2κ_∞).

To compute the true accuracy CV_∞ in the long term limit for the case $\hat{D}\ne D$, ${\mathrm{C}\mathrm{V}}_{\infty }=\underset{m\to \infty }{\mathrm{lim}}\enspace \mathrm{C}\mathrm{V}\left[p\left({\varepsilon }_{m}\right)\right]$, we note that the iteration rule for CV[p(ɛ_m )], equation (22), still holds. Taking the limit m → ∞ yields a new self-consistency equation for CV_∞. A short calculation yields the unique solution, see appendix D

$$\begin{equation}\text{Case}\;D\ne \hat{D}:{\vert c\vert } {\mathrm{C}\mathrm{V}}_{\infty }\approx \frac{1}{2}\sqrt{\frac{2\enspace D{\Delta}t}{\mathrm{S}\mathrm{N}\mathrm{R}}}\cdot \frac{D+\hat{D}}{2\sqrt{D\hat{D}}}.\end{equation} \tag{ 24 }$$

This expression attains the minimal value CV_∞ = [DΔt/(2SNR)]^1/2 exactly for $\hat{D}=D$. Equation (24) compares favorably to simulation results, see figure 3. For details on numerical methods, see appendix A. We conclude that Bayesian gradient estimation can still operate close to its optimum even if the agent has only inaccurate knowledge of its true rotational diffusion coefficient, provided the true and the assumed rotational diffusion coefficient, D and $\hat{D}$, are of similar magnitude.

Zoom In Zoom Out Reset image size

However, the estimation of gradient direction fails if the agent would completely neglect its rotational diffusion, corresponding to $\hat{D}=0$, as we show next.

Agents unaware of their rotational diffusion: D > 0, $\hat{D}=0$. If agents were unaware of their own rotational diffusion with $\hat{D}=0$, they would perform only the update step, equation (15). Correspondingly, the likelihood distribution ${\mathcal{L}}_{m}\left(\psi \right)$ for the gradient direction angle at time t_m becomes a product of m + 1 von-Mises distributions. The calculus of von-Mises distribution yields an expression for the measure-of-concentration κ_m of this likelihood distribution as a circular average of all previous measurements

$$\begin{equation}{\kappa }_{m}\enspace {\mathrm{e}}^{\mathrm{i}{\hat{\psi }}_{m}}={\kappa }_{0}^{\prime }\enspace {\mathrm{e}}^{\mathrm{i}{\hat{\psi }}_{0}}+\sum\limits _{j=1}^{m}{\zeta }_{j}\enspace {\mathrm{e}}^{\mathrm{i}{\psi }_{j}}.\end{equation} \tag{ 25 }$$

Here, the first term represents the contribution from the Bayesian prior, whereas the sum represents the contribution from the sequential measurements. For simplicity, we assume a flat prior with κ'₀ = 0. We can rewrite equation (25) in terms of the complex measurement ${\mathcal{M}}_{j}$ as ${\kappa }_{m}\enspace \mathrm{exp}\left(i{\hat{\psi }}_{m}\right)=\alpha {\sum }_{j}{\;\;\mathcal{M}}_{j}$. To evaluate this sum, we note the correlations between subsequent measurements, stemming from the rotational diffusion of the agent

$$\begin{equation}\langle {\mathcal{M}}_{j}{\mathcal{M}}_{k}^{{\ast}}\rangle ={\left(\frac{\alpha \nu }{2}\right)}^{2}\enspace {\mathrm{e}}^{-D\vert {t}_{k}-{t}_{j}\vert }\quad \;\text{for}\;\enspace j\ne k\quad .\end{equation} \tag{ 26 }$$

We can now compute $\langle {\kappa }_{m}^{2}\rangle$ in a straight-forward manner, noting that the term ${\sum }_{j\ne k}\langle {\mathcal{M}}_{j}{\mathcal{M}}_{k}^{{\ast}}\rangle$ yields a double-geometric series, see appendix E for details and the (lengthy) result. Using the approximation ${\kappa }_{m}\approx {\langle {\kappa }_{m}^{2}\rangle }^{1/2}$, we obtain the asymptotic scaling

$$\begin{equation}\text{Case}\;\hat{D}=0:{\vert c\vert } {\kappa }_{m}\approx \begin{cases}r{t}_{m}\quad \;\text{for}\;\;\enspace {t}_{m}\ll \tau \quad \hfill \\ 2r\sqrt{{t}_{m}\tau }\quad \;\text{for}\;{t}_{m}\gg \tau \quad \hfill \end{cases}.\end{equation} \tag{ 27 }$$

Hence, for D > 0 but $\hat{D}=0$, the estimated precision κ_m will initially increase linearly with m, but then cross-over to the asymptotic scaling κ_m ∼ m^1/2 after a characteristic cross-over time on the order of τ = (2D)⁻¹. Simultaneously, CV[p(ɛ_m )] → 1 for t_m ≫ τ. These simple scaling relations are corroborated by numerical simulations, see figure 2(B). For details on numerical methods, see appendix A.

Thus, agents not aware of their own rotational diffusion will arrive at erroneous gradient estimates, because past measurements will have become partially invalidated by rotational diffusion, yet are nonetheless incorporated by the agent in its gradient estimates with full weight. Concomitantly, the estimated precisions κ_m that individual agents assign to their estimates of gradient direction do not reflect the true accuracy, i.e. the dispersion of estimation errors in an ensemble of agents. Individual agents are 'over-confident' of their own estimates.

Case $D=\hat{D}=0$. Finally, in the absence of rotational diffusion with $D=\hat{D}=0$, the estimated precision κ_m grows linearly with time as ⟨κ_m ⟩ ≈ mSNR = rt_m , as expected. Correspondingly, the true accuracy of estimation errors in an ensemble of agents decreases inversely with total measurement time, CV_m ≈ 1/(2rt_m ). This is Berg-and-Purcell's classical result for the estimation of absolute concentration [7], generalized to the problem of estimating gradient direction.

Summary of results. We considered a minimal model of gradient sensing in the presence of both sensing and motility noise. We analytically derived the theoretical limit of gradient sensing for a chemotactic agent that undergoes rotational diffusion. The chemotactic agent is placed in an external concentration gradient and detects stochastic binding events of signaling molecules on its circumference, while it undergoes rotational diffusion with rotational diffusion coefficient D.

We find that the accuracy by which this agent can estimate the time-varying relative direction of the gradient depends on two fundamental time-scales: the rotational diffusion time τ = 1/(2D) on which the relative gradient direction changes, and a second time-scale set by a rate of information gain r = SNR/Δt defined as a signal-to-noise ratio of gradient sensing normalized by unit time. The error ɛ between true and estimated gradient direction is distributed at steady-state with a circular variance ${\mathrm{C}\mathrm{V}}_{\infty }\left[p\left(\varepsilon \right)\right]\approx {\left(r{T}_{\mathrm{e}\mathrm{ff}}\right)}^{-1}$ that scales inversely with an effective measurement time ${T}_{\mathrm{e}\mathrm{ff}}=\sqrt{\tau /r}$ given by the geometric mean of the two time-scales of the problem. Numerical simulations corroborate this scaling relation, see figure 4. For details on numerical methods, see appendix A. Intuitively, the agent averages noisy measurements only over the recent past defined by a time interval of duration T_eff. The effective measurement time T_eff reflects a trade-off between the requirement to average noisy measurements over long time intervals to reduce sensing noise, and a second requirement to use only recent measurements because the relative gradients direction already changes during the measurement interval due to motility noise.

Zoom In Zoom Out Reset image size

For the case of optimal estimation, we assume that the agent uses sequential Bayesian estimation and has perfect knowledge of its own rotational diffusion coefficient. In this case, the precision that individual agents assign to their own gradient estimates matches the accuracy of gradient estimates in an ensemble of identical agents, i.e. gradient estimation is faithful. If the agent assumes a rotational diffusion coefficient $\hat{D}$ that is either too small or too large, this is no longer the case, and the circular variance CV_∞[p(ɛ)] of estimation errors will be larger than in the optimal case $D=\hat{D}$. If, however, the agent would fully ignore its owns rotational diffusion, corresponding to $\hat{D}=0$, gradient sensing will fail in a characteristic way: agents not aware of their rotational diffusion will extend their temporal averaging infinitely into the past. Concomitantly, the estimated gradient direction decorrelates from the true direction on the time-scale τ of rotational diffusion, while agents erroneously believe that their estimates become more and more precise as function of total measurement time t. We find an abnormal asymptotic scaling of the variance of gradient estimates, which decreases with time as $\mathrm{C}\mathrm{V}\left[\mathcal{L}\left(\psi \right)\right]\sim {t}^{-1/2}$ (estimated precision), while the errors between true and estimated gradient direction converges to one CV[p(ɛ)] → 1 (true accuracy). The abnormal scaling of the estimated precision represents a signature of erroneous state estimation, and could be tested in real-world applications, e.g. of bearing tracking.

Cause of the square-root scaling of optimal measurement time. Our theory shows that the optimal measurement time for the estimation of a time-varying gradient direction scales as T_eff = (τ/r)^1/2, while the variance of the expected estimation error scales as 1/(rT_eff) = (rτ)^−1/2. This result parallels a recent result by Mora et al on the optimal estimation of an absolute concentration that varies in time [25]. In both Mora et al and our work, τ and r denote a characteristic fluctuation time-scale, and a rate of information gain of the problem, respectively.

We propose that the inverse square-root scaling of the optimal measurement time could be a general feature of Bayesian estimation of time-varying signals. The following heuristic argument highlights the common idea: if an agent averages past measurements of a time-varying signal s*(t) with an effective measurement time T_eff, this agent effectively estimates $\bar{{s}^{{\ast}}}\left(t\right)={T}_{\mathrm{e}\mathrm{ff}}^{-1}{\int }_{t-{T}_{\mathrm{e}\mathrm{f}\mathrm{f}}}^{t} \mathrm{d}{t}^{\prime }{s}^{{\ast}}\left({t}^{\prime }\right)\approx {s}^{{\ast}}\left(t-{T}_{\mathrm{e}\mathrm{ff}}/2\right)$, i.e. the estimate lags behind the current value of the signal. The variance σ² of the expected estimation error $\varepsilon \left(t\right)=\hat{s}\left(t\right)-{s}^{{\ast}}\left(t\right)$ between the estimated value $\hat{s}\left(t\right)$ of the signal and its true value s*(t) can be approximately decomposed into two contributions: ${\sigma }_{\text{sensing}}^{2}=\langle {\left[\hat{s}-\bar{{s}^{{\ast}}}\enspace \right]}^{2}\rangle$ and ${\sigma }_{\text{time}}^{2}=\langle {\left[{s}^{{\ast}}-\bar{{s}^{{\ast}}}\right]}^{2}\rangle$. The first contribution predominantly reflects sensing noise, and approximately scales as ${\sigma }_{\text{sensing}}^{2}\sim 1/\left(r{T}_{\mathrm{e}\mathrm{ff}}\right)$ as function of effective measurement time T_eff. The second contribution accounts for the stochastic change of the signal during the time T_eff and approximately scales as ${\sigma }_{\text{time}}^{2}\sim {T}_{\mathrm{e}\mathrm{ff}}/\tau$. While the first contribution decreases with T_eff, the second summand increases, thus posing an optimization problem for the optimal T_eff that minimizes the total variance of the error ${\sigma }^{2}\approx {\sigma }_{\text{sensing}}^{2}+{\sigma }_{\text{time}}^{2}$. Solving for dσ²/dT_eff = 0 yields T_eff ∼ (τ/r)^1/2. A detailed version of this argument can be found in appendix F.

Difference to bacterial chemotaxis. We emphasize that our and Mora's recent results on the optimal effective measurement time for estimation of time-varying signals are different from previous work [26]. Previously, Strong et al extended the seminal work by Berg and Purcell [7] to compute the optimal filter for bacterial chemotaxis, i.e. the optimal sensori-motor transfer function for bacteria with run-and-tumble motility in the presence of sensing and motility noise. This filter extended temporal comparison over a time span set by the rotational diffusion time, i.e. T_eff ∼ D⁻¹, (see equation (73) in loc. cit.), and additionally had an extremely long tail (see equation (70) in loc. cit.). This difference may be related to the fact that bacteria like E. coli use a chemotaxis strategy based on temporal comparison, which is different from the strategy of spatial comparison considered here: these bacteria estimate a time derivative of a concentration signal traced along their swimming path [19]. In contrast, most eukaryotic cells use spatial comparison, i.e. estimate spatial concentration gradients across the cell diameter [4, 5, 41–43], as considered here. A third chemotaxis strategy represented by marine sperm cells navigating along helical paths [3], can be mapped on the case of spatial comparison: while moving along helical paths, these cells 'visit' different sensor positions during one helical turn.

Directional memory and rotational diffusion. Biological cells performing chemotaxis possess internal memory of recent measurements of gradient direction. In eukaryotic cells with crawling cell motility, this memory manifests itself, e.g. in a persistent polarization of cell shape and the acto-myosin cytoskeleton, as well as anisotropic distributions of chemotactic effector molecules at the cell boundary [44, 45]. Stochastic remodeling of the cytoskeleton, as well as stochastic formation of pseudopods, would cause a decorrelation of cell polarization in the absence of chemotactic guidance cues. For marine sperm cells,the axis of helical swimming paths represents a consolidated memory of previous noisy steering responses [46]. In our minimal, coarse-grained model, we capture all effects that randomize direction in terms of a single effective rotational diffusion coefficient.

Biochemical implementation. Storing the likelihood distribution of estimated gradient directions, or just a proxy thereof, requires internal memory. Recent mathematical models of cellular gradient sensing such as LEGI (local excitation, global inhibition) [47], or balanced-inactivation models [48], consider distributions of chemotactic effector molecules on the cell boundary, which could serve as a proxy for a likelihood distribution of gradient direction: while the position of a concentration peak in such a distribution can represent a maximum likelihood estimate, its amplitude could encode a level of certainty. Slow surface diffusion of effector molecules would relax their distribution toward a uniform distribution in a way that is formally equivalent to a prediction step that accounts for the loss of knowledge due to rotational diffusion.

Typical parameters. Typical rotational diffusion coefficients for the bacterium E. coli are D ∼ 0.1 s⁻¹, close to the theoretical lower limit of a passive particle of same size and shape. For ten-fold larger sperm cells, active fluctuations dominate [49], resulting in an estimate D ∼ 0.01–0.1 s⁻¹ [40, 46]. The motility of crawling Dictyostelium cells was characterized by a persistence time of ∼5 min in the absence of chemoactractant [50], which sets an effective rotational diffusion coefficient D ∼ 0.003 s⁻¹. For immune T cells in three-dimensional tissue, a persistence time of ∼1 min was found (displaying a characteristic speed dependence) [51]. If binding of signaling molecules to receptors on the cell surface is diffusion-limited, we can estimate the rate of binding by λc₀ = 4πD_c ac₀ for a perfectly absorbing spherical cell of radius a, where D_c denotes the translational diffusion coefficient of signaling molecules [52]. For typical values (D_c ∼ 300 μm² s⁻¹, a ∼ 10 μm, c₀ ∼ 1 nM), we estimate λc₀ ∼ 10⁴ s⁻¹. Thus, for a concentration gradient of either α₀ = 1%, 0.5%, or 0.1% across the diameter of a cell, we estimate a rate of information gain of r ≈ 0.5 s⁻¹, ≈0.12 s⁻¹ and ≈0.005 s⁻¹, respectively, or, equivalently, a signal-to-noise ratio of gradient sensing SNR ≈ 0.5, ≈0.12, and ≈0.005 for a time interval Δt = 1 s. Assuming D ∼ 0.003 s⁻¹ [50], our main result, equation (19), predicts for the ultimate precision of gradient sensing CV_∞ ≈ 0.06, ≈0.11, and ≈0.55 for these three gradient strengths, respectively. Reversible binding of signaling molecules effectively increases sensing noise, thus decreasing the signal-to-noise ratio by a constant prefactor [8, 9, 12]. Some cells, such as sperm cells, respond chemotactically even at pico-molar concentrations [53], corresponding to even lower rates of information gain [20, 39]. For the slime mold dictyostelium, the precision of chemotactic migration is set by sensing noise at low rates of information gain (corresponding to low signal-to-noise ratios), while noise of downstream intracellular signaling [54] becomes relevant at high rates of information gain [12, 16, 17]. Our theory of Bayesian gradient estimation provides a theoretical limit for the precision of gradient sensing, even if some cells do not achieve this limit, e.g. because these cells use shorter effective measurement times.

Time-varying gradients. The minimal model of a stationary agent subject to rotational diffusion considered here represents a prototypical example for a broader class of problems: if the agent actively moves in a concentration landscape, also the absolute concentration of chemoattractant and the strength of the gradient will change as the agent moves. Thus, the agent faces two estimation problems at the same time: estimating absolute concentration, which is comparatively easy, and estimating the direction of a weak concentration gradient, which is considerably more difficult. Our minimal model decouples these two problems and studies the second, more difficult problem in isolation. Similarly, in a dynamic concentration landscape, the direction of the concentration gradient will change in time with a characteristic fluctuation time-scale τ. For an idealized case of random changes of gradient direction, this case of time-varying concentration gradients is formally equivalent to the case of rotational diffusion with rotational diffusion time τ considered here.

No new data were created or analysed in this study.

MN and BMF are supported by the German National Science Foundation (DFG) through the Excellence Initiative by the German Federal and State Governments (Clusters of Excellence cfaed EXC-1056 and PoL EXC-2068), as well as DFG grant FR3429/3-1 and FR3429/4-1 (Heisenberg program) to BMF. We thank all members of the 'Biological Algorithms Group', in particular Andrea Auconi and Steffen Lange, for stimulating discussions.

For numerical simulations, we represented likelihood distributions ${\mathcal{L}}_{m}\left(\psi \right)$ using an equidistant grid for ψ with step size of Δψ = 0.01. For the case of an agent unaware of its own rotational diffusion (D > 0, $\hat{D}=0$), where likelihood distributions become very sharp, we used a locally refined grid: we first identified a region of interest (ROI), where the likelihood distribution exceeded a relative threshold of $1{0}^{-7}\mathrm{max}\left({\mathcal{L}}_{m}\left(\psi \right)\right)$, and used two equidistant grids, inside and outside of this ROI, with 3000 grid points each.

In each time step, the likelihood distributions were updated according to the update step, equation (15), and the prediction step, equation (17). We used the diffusion approximation, equation (11), to sample stochastic measurements ${\mathcal{M}}_{j}$. Drawing stochastic measurements from a Poisson distribution instead gave virtually identical results in initial simulations (same circular variance and Shannon information of measurements even for small SNR). Likewise, we confirmed that results do not change if we use a finer grid spacing for ψ.

To numerically compute the true accuracy CV_∞ at steady state shown in figure 3, we simulated agents for m = 100 time steps and averaged CV[p(ɛ_m )] over the last 25 time steps. In all cases, we confirmed that CV[p(ɛ_m )] had converged and was fluctuating around a limit value.

To compute the average estimated precision ${\bar{\kappa }}_{\infty }$ in figure 4, we first determined the circular variances $\mathrm{C}\mathrm{V}\left[{\mathcal{L}}_{m}\left(\psi \right)\right]$ and $\mathrm{C}\mathrm{V}\left[{\mathcal{L}}_{m}^{\prime }\left(\psi \right)\right]$ of simulated likelihood distributions after update and prediction step, respectively, and then converted these circular variances to equivalent measures-of-concentration κ_m and κ'_m by inverting CV = 1 − I₁(κ)/I₀(κ). Finally, ${\bar{\kappa }}_{\infty }$ is obtained as limit value of the geometric mean $\sqrt{{\kappa }_{m}{\kappa }_{m}^{\prime }}$ (determined as average over the last 25 of a total of 100 time steps, when all measures-of-concentration have converged). Similarly, to compute the true accuracy, we first determined the circular variances CV[p(ɛ_m )] and CV[p(ɛ'_m )] of the distribution of estimation errors in an ensemble of agents after update and prediction step, respectively, converted these again to equivalent measures-of-concentration, and finally determined the limit value of their geometric mean.

Error bars determined by bootstrapping are smaller than symbol size.

To prove equation (12) stating $\langle {\zeta }_{j}\rangle =\alpha \langle \vert {\mathcal{M}}_{j}\vert \rangle \approx \mathrm{S}\mathrm{N}\mathrm{R}$, we compute the expectation value $\langle \vert {\mathcal{M}}_{j}\vert \rangle$ using the law of cosines. We write $\rho ={\mathcal{M}}_{j}-\langle {\mathcal{M}}_{j}\rangle$ for the deviation of ${\mathcal{M}}_{j}$ from its mean $\mu =\langle {\mathcal{M}}_{j}\rangle =\alpha \nu /2\enspace \mathrm{exp}\enspace i{\psi }_{j}^{{\ast}}$ with polar representation $\rho =b\enspace \mathrm{exp}\enspace i\left({\psi }_{j}^{{\ast}}+\gamma \right)$. We thus have a triangle in the complex plane with vertices $0,\mu ,{\mathcal{M}}_{j}\in \mathbb{C}$ and side lengths $a=\vert \langle {\mathcal{M}}_{j}\rangle \vert$, b = |ρ|, $c=\vert {\mathcal{M}}_{j}\vert$, where γ is the angle enclosed by the sides of length a and b. The law of cosines for this triangle reads

$$\begin{equation}c=\sqrt{{a}^{2}+{b}^{2}-2ab\enspace \mathrm{cos}\enspace \gamma }.\end{equation} \tag{ S1 }$$

According to equation (11), the deviations ρ are isotropic, hence the angle γ is a uniformly distributed with p(γ) = 1/(2π). The squared magnitude of the deviation b² = |ρ|² follows a χ²-distribution for 2 degrees of freedom (namely, real and imaginary part of ρ, respectively); hence, p(b) = 2b/ν exp(−b²/ν).

Thus,

$$\begin{equation}\langle \vert {\mathcal{M}}_{j}\vert \rangle =\langle c\rangle ={\int }_{0}^{\infty } \mathrm{d}b\enspace p{\left(b\right)\oint }_{0}^{2\pi } \mathrm{d}\varphi \enspace p\left(\varphi \right)\enspace c\left(b,\varphi \right).\end{equation} \tag{ S2 }$$

The first integration, $\mathcal{I}\left(b\right){=\oint }_{0}^{2\pi } \mathrm{d}\varphi \enspace p\left(\varphi \right)c\left(b,\varphi \right)$, results in an elliptic integral (explicitly, $2\pi \mathcal{I}\left(b\right)={\pm}\left(\alpha \nu -2b\right)E\left[-8\alpha \nu b/{\left(\alpha \nu -2b\right)}^{2}\right]+\left(\alpha \nu +2b\right)E\left[8\alpha \nu b/{\left(\alpha \nu +2b\right)}^{2}\right]$, where E denotes the complete elliptic integral of the second kind, and the '+' sign applies for b < αν/2 and '−' otherwise). The second integration in equation (S2) with respect to b cannot be performed analytically. However, $\mathcal{I}\left(b\right)$ can be well approximated by

$$\begin{equation}\mathcal{I}\left(b\right)\approx \begin{cases}\frac{\alpha \nu }{2}\left[1+{\left(\frac{b}{\alpha \nu }\right)}^{2}\right]\quad b{\leqslant}\alpha \nu /2\quad \hfill \\ b\left[1+{\left(\frac{\alpha \nu }{4b}\right)}^{2}\right]\quad b{ >}\alpha \nu /2\quad \hfill \end{cases}.\end{equation} \tag{ S3 }$$

The second integration with respect to b can now be easily done, yielding

$$\begin{equation}\alpha \langle \vert {\mathcal{M}}_{j}\vert \rangle =\alpha \langle c\rangle \approx \mathrm{S}\mathrm{N}\mathrm{R}+\frac{1}{2}+\cdots \enspace ,\end{equation} \tag{ S4 }$$

where the ellipses represent terms that decay exponentially fast for SNR ≫ 1 (explicitly, $\alpha \langle c\rangle =\mathrm{S}\mathrm{N}\mathrm{R}+\left(1/2\right)-\left[\left(2+\mathrm{S}\mathrm{N}\mathrm{R}\right)/4\right]\enspace \mathrm{exp}\left(-\mathrm{S}\mathrm{N}\mathrm{R}/2\right)+\sqrt{2\pi \mathrm{S}\mathrm{N}\mathrm{R}}\left[\left(4+\mathrm{S}\mathrm{N}\mathrm{R}\right)/8\right]\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}\left(\sqrt{\mathrm{S}\mathrm{N}\mathrm{R}}\right)$).

A probability distribution p(ψ) of angles should be 2π-periodic, i.e. p(ψ) = p(ψ + 2π), and normalized to one on the unit circle, i.e. ${\oint }_{0}^{2\pi } \mathrm{d}\psi \enspace p\left(\psi \right)=1$. The circular variance of such a circular distribution is defined as

$$\begin{equation}\mathrm{C}\mathrm{V}\left[p\left(\psi \right)\right]=1-\left\vert {\oint }_{0}^{2\pi } \mathrm{d}\psi \enspace {\mathrm{e}}^{\mathrm{i}\psi }p\left(\psi \right)\right\vert .\end{equation} \tag{ S5 }$$

An important circular distribution is the wrapped normal distribution

$$\begin{equation}\mathcal{WN}\left(\psi ;\mu ,{\sigma }^{2}\right)=\sum\limits _{l=-\infty }^{\infty }{\left(2\pi {\sigma }^{2}\right)}^{-1/2}\enspace \mathrm{exp}\left(-\frac{{\left(\psi -\mu +2\pi l\right)}^{2}}{2{\sigma }^{2}}\right)\end{equation} \tag{ S6 }$$

with variance parameter σ², whose circular variance reads $\mathrm{C}\mathrm{V}=1-{\mathrm{e}}^{-{\sigma }^{2}/2}$. In the limit σ² ≪ 1, CV ≈ σ²/2. The wrapped normal distribution is closely approximated by the von-Mises distribution, which is commonly used in directional statistics due to its mathematical tractability [55]

$$\begin{equation}\mathcal{VM}\left(\psi ;\mu ,\kappa \right)={\left(2\pi {I}_{0}\left(\kappa \right)\right)}^{-1}\enspace \mathrm{exp}\left[\kappa \enspace \mathrm{cos}\left(\psi -\mu \right)\right],\end{equation} \tag{ S7 }$$

where κ is the so-called measure-of-concentration or precision, and I_n (κ) is the modified Bessel function of order n. The circular variance reads $\mathrm{C}\mathrm{V}=1-{I}_{1}\left(\kappa \right)/{I}_{0}\left(\kappa \right)=1/\left(2\kappa \right)+1/\left(8{\kappa }^{2}\right)+\mathcal{O}\left({\kappa }^{-3}\right)$.

The normalized product of two von-Mises distributions, say $\mathcal{VM}\left(\psi ;{\mu }_{1},{\kappa }_{1}\right)$ and $\mathcal{VM}\left(\psi ;{\mu }_{2},{\kappa }_{2}\right)$, is again a von-Mises distribution $\mathcal{VM}\left(\psi ;\mu ,\kappa \right)$. Such normalized product appears, e.g. in Bayes formula, equation (15). Specifically, the mean μ and measure-of-concentration κ of the normalized product satisfy

$$\begin{equation}\kappa \enspace {\mathrm{e}}^{\mathrm{i}\mu }={\kappa }_{1}\enspace {\mathrm{e}}^{i{\mu }_{1}}+{\kappa }_{2}\enspace {\mathrm{e}}^{\mathrm{i}{\mu }_{2}}.\end{equation} \tag{ S8 }$$

If mean values are close, |μ₁ − μ₂| ≪ 1, we have the approximate sum rule κ ≈ κ₁ + κ₂.

In contrast, the circular convolution of two von-Mises distributions $\mathcal{VM}\left(\psi ;{\mu }_{1},{\kappa }_{1}\right)$ and $\mathcal{VM}\left(\psi ;{\mu }_{2},{\kappa }_{2}\right)$ is only approximately a von-Mises distribution [55]. To find such approximation, one can first map the two von-Mises distributions onto wrapped normal distributions of same respective mean and circular variance, and compute the convolution of these wrapped normal distributions [38]. The convolution of two wrapped normal distributions, say with variance parameters ${\sigma }_{1}^{2}$ and ${\sigma }_{2}^{2}$, is again a wrapped normal distribution with new variance parameter ${\sigma }^{2}={\sigma }_{1}^{2}+{\sigma }_{2}^{2}$. Finally, this new wrapped normal distribution $\mathcal{WN}\left(\psi ;{\mu }_{1}+{\mu }_{2},{\sigma }^{2}\right)$ is mapped back onto a von-Mises distribution $\mathcal{VM}\left(\psi ;\mu ,\kappa \right)$. Thus, $\mathcal{VM}\left(\psi ;\mu ,\kappa \right)\approx \mathcal{VM}{\left(\psi ;{\mu }_{1},{\kappa }_{1}\right)}^{{\ast}}\mathcal{VM}\left(\psi ;{\mu }_{2},{\kappa }_{2}\right)$ with μ = μ₁ + μ₂ and ${I}_{1}\left(\kappa \right)/{I}_{0}\left(\kappa \right)=\mathrm{exp}\left(-{\sigma }^{2}/2\right)=\mathrm{exp}\left(-{\sigma }_{1}^{2}\right)\mathrm{exp}\left(-{\sigma }_{2}^{2}\right)=\left[{I}_{1}\left({\kappa }_{1}\right)/{I}_{0}\left({\kappa }_{1}\right)\right]\cdot \left[{I}_{1}\left({\kappa }_{2}\right)/{I}_{0}\left({\kappa }_{2}\right)\right]$. In the limit κ₁, κ₂ ≫ 1, we have ${\kappa }^{-1}={\kappa }_{1}^{-1}+{\kappa }_{2}^{-1}$.

For ease of reference, we highlight the sum rules for the measure-of-concentration of either a normalized product or convolution of two von-Mises distributions

$$\begin{equation}\text{convolution:}\;\;\frac{1}{\kappa }\approx \frac{1}{{\kappa }_{1}}+\frac{1}{{\kappa }_{2}},\end{equation} \tag{ S9 }$$

$$\begin{equation}\text{multiplication:}\;\;\kappa \approx {\kappa }_{1}+{\kappa }_{2}.\end{equation} \tag{ S10 }$$

While equation (S9) is valid for κ₁, κ₂ ≫ 1, equation (S10) is valid for |μ₁ − μ₂| ≪ 1. Table C presents a summary of this section. Figure S1 graphically relates von-Mises distributions and wrapped normal distributions of equal circular variance.

Table C. Mathematical properties of normal distribution, von-Mises distribution, and wrapped normal distribution. See appendix C for details.

	Normal distribution	von-Mises distribution	Wrapped normal distribution
Definition	$\mathcal{N}\left(\psi ;\mu ,{\sigma }^{2}\right)=\frac{{\mathrm{e}}^{{\left(\psi -\mu \right)}^{2}/\left(2{\sigma }^{2}\right)}}{\sqrt{2\pi {\sigma }^{2}}}$	$\mathcal{VM}\left(\psi ;\mu ,\kappa \right)=\frac{{\mathrm{e}}^{\kappa \enspace \mathrm{cos}\left(\psi -\mu \right)}}{2\pi {I}_{0}\left(\kappa \right)}$	$\mathcal{WN}\left(\psi ;\mu ,{\sigma }^{2}\right)=\sum _{l}\mathcal{N}\left(\psi ;\mu +2\pi l,{\sigma }^{2}\right)$
Normalized	$\mathcal{N}\left({\mu }_{1},{\sigma }_{1}^{2}\right)\mathcal{N}\left({\mu }_{2},{\sigma }_{2}^{2}\right)\sim \mathcal{N}\left(\mu ,{\sigma }^{2}\right)$	$\mathcal{VM}\left({\mu }_{1},{\kappa }_{1}\right)\mathcal{VM}\left({\mu }_{2},{\kappa }_{2}\right)\sim \mathcal{VM}\left(\mu ,\kappa \right)$	$\mathcal{WN}\left({\mu }_{1},{\sigma }_{1}^{2}\right)\mathcal{WN}\left({\mu }_{2},{\sigma }_{2}^{2}\right) \approx \mathcal{WN}\left(\mu ,{\sigma }^{2}\right)$
	product
	$\frac{\mu }{{\sigma }^{2}}=\frac{{\mu }_{1}}{{\sigma }_{1}^{2}}+\frac{{\mu }_{2}}{{\sigma }_{2}^{2}}$	$\kappa \enspace {\mathrm{e}}^{\mathrm{i}\mu }={\kappa }_{1}\enspace {\mathrm{e}}^{\mathrm{i}{\mu }_{1}}+{\kappa }_{2}\enspace {\mathrm{e}}^{\mathrm{i}{\mu }_{2}}$	$\frac{\mu }{{\sigma }^{2}} = \frac{{\mu }_{1}}{{\sigma }_{1}^{2}}+\frac{{\mu }_{2}}{{\sigma }_{2}^{2}}$ for σ1,2 ≪ 1, |μ1 − μ2|≪ 1
$\frac{1}{{\sigma }^{2}}=\frac{1}{{\sigma }_{1}^{2}}+\frac{1}{{\sigma }_{2}^{2}}$	κ ≈ κ1 + κ2 for |μ1 − μ2| ≪ 1	$\frac{1}{{\sigma }^{2}} = \frac{1}{{\sigma }_{1}^{2}}+\frac{1}{{\sigma }_{2}^{2}}$ for σ1,2 ≪ 1, |μ1 − μ2|≪ 1
		equations (S10) and (S8)
Convolution	$\mathcal{N}\left({\mu }_{1},{\sigma }_{1}^{2}\right)\;{\ast}\;\mathcal{N}\left({\mu }_{1},{\sigma }_{1}^{2}\right)=\mathcal{N}\left(\mu ,{\sigma }^{2}\right)$	$\mathcal{VM}\left({\mu }_{1},{\kappa }_{1}\right)\;{\ast}\;\mathcal{VM}\left({\mu }_{2},{\kappa }_{2}\right)\approx \mathcal{VM}\left(\mu ,\kappa \right)$	$\mathcal{WN}\left({\mu }_{1},{\sigma }_{1}^{2}\right)\;{\ast}\;\mathcal{WN}\left({\mu }_{1},{\sigma }_{1}^{2}\right) = \mathcal{WN}\left(\mu ,{\sigma }^{2}\right)$
	μ = μ1 + μ2	μ = μ1 + μ2	μ = μ1 + μ2
	$\sigma ={\sigma }_{1}^{2}+{\sigma }_{2}^{2}$	$\frac{1}{\kappa }\approx \frac{1}{{\kappa }_{1}}+\frac{1}{{\kappa }_{2}}$ for κ1, κ2 ≫ 1 equation (S9)	$\sigma ={\sigma }_{1}^{2}+{\sigma }_{2}^{2}$
Circular	—	1 − I1(κ)/I0(κ) ≈ 1/(2κ) for κ ≫ 1	$1-{\mathrm{e}}^{-{\sigma }^{2}/2}$ ≈σ2/2 for σ2 ≪ 1
variance CV

Zoom In Zoom Out Reset image size

We provide details on the calculation of the true accuracy of gradient sensing in the long time limit, CV_∞ = lim_m→∞CV[p(ɛ_m )], for the case that agents assume a wrong rotational diffusion coefficient $\hat{D}\ne D$.

First, we note that the iteration rule for CV_m , equation (22), still holds in the case $\hat{D}\ne D$. We take the limit m → ∞ on both sides an obtain a self-consistency equation for CV_∞

$$\begin{align}\hfill {\mathrm{C}\mathrm{V}}_{\infty }& \approx {\left(\frac{{\kappa }_{\infty }^{\prime }}{{\kappa }_{\infty }}\right)}^{2}\left({\mathrm{C}\mathrm{V}}_{\infty }+D{\Delta}t\right)+{\left(\frac{\mathrm{S}\mathrm{N}\mathrm{R}}{{\kappa }_{\infty }}\right)}^{2}\frac{1}{2\;\mathrm{S}\mathrm{N}\mathrm{R}}\hfill \\ \hfill & \approx \frac{{\kappa }_{\infty }^{2}-2\;\mathrm{S}\mathrm{N}\mathrm{R}{\kappa }_{\infty }}{{\kappa }_{\infty }^{2}}\left({\mathrm{C}\mathrm{V}}_{\infty }+D{\Delta}t\right)+{\left(\frac{\mathrm{S}\mathrm{N}\mathrm{R}}{{\kappa }_{\infty }}\right)}^{2}\frac{1}{2\;\mathrm{S}\mathrm{N}\mathrm{R}},\hfill \end{align} \tag{ S11 }$$

where we used κ'_∞ = lim_m→∞ κ'_m ≈ κ_∞ − SNR. Note that κ_∞, κ'_∞ ≫ SNR by equation (6). In the main text, we had shown for the estimated precision κ_∞ in the long time limit, ${\kappa }_{\infty }={\left[\mathrm{S}\mathrm{N}\mathrm{R}/\left(2\hat{D}{\Delta}t\right)\right]}^{1/2}$, hence $\mathrm{S}\mathrm{N}\mathrm{R}/{\kappa }_{\infty }={\left(\mathrm{S}\mathrm{N}\mathrm{R}\cdot 2\hat{D}{\Delta}t\right)}^{1/2}$. Inserting this into equation (S11) gives

$$\begin{equation}{\mathrm{C}\mathrm{V}}_{\infty }\approx \left(1-2\sqrt{\mathrm{S}\mathrm{N}\mathrm{R}\cdot 2\hat{D}{\Delta}t}\right)\left({\mathrm{C}\mathrm{V}}_{\infty }+D{\Delta}t\right)+\hat{D}{\Delta}t,\end{equation} \tag{ S12 }$$

or after rearranging

$$\begin{equation}1-2\sqrt{\mathrm{S}\mathrm{N}\mathrm{R}\cdot 2\hat{D}{\Delta}t}=\frac{{\mathrm{C}\mathrm{V}}_{\infty }-\hat{D}{\Delta}t}{{\mathrm{C}\mathrm{V}}_{\infty }+D{\Delta}t}\approx 1-\frac{\hat{D}+D}{{\mathrm{C}\mathrm{V}}_{\infty }}{\Delta}t.\end{equation} \tag{ S13 }$$

Solving for CV_∞ yields the result stated in equation (24), ${\mathrm{C}\mathrm{V}}_{\infty }\approx \frac{1}{2}\sqrt{2D{\Delta}t/\mathrm{S}\mathrm{N}\mathrm{R}}\cdot \left(D+\hat{D}\right)/\sqrt{4D\hat{D}}$.

We compute the second moment $\langle {\kappa }_{m}^{2}\rangle$ of the measure-of-concentration κ_m of the likelihood distribution ${\mathcal{L}}_{m}\left(\psi \right)$ for the gradient angle ψ after m measurements.

Using the calculus of von-Mises distributions, we can express κ_m as a circular average of all previous measurements, see equation (25), which reads ${\kappa }_{m}\enspace \mathrm{exp}\left(i{\hat{\psi }}_{m}\right)={\kappa }_{0}^{\prime }\enspace \mathrm{exp}\left(i{\hat{\psi }}_{0}\right)+{\sum }_{j=1}^{m}{\zeta }_{j}\enspace \mathrm{exp}\left(i{\psi }_{j}\right)$. For simplicity, we assume a flat initial prior with κ'₀ = 0. Using ζ_j = α|M_j |, we can rewrite this equation in terms of the complex measurement ${\mathcal{M}}_{j}$ as

$$\begin{equation}{\kappa }_{m}\enspace {\mathrm{e}}^{\mathrm{i}{\hat{\psi }}_{m}}=\alpha \sum\limits _{j=1}^{m}{\mathcal{M}}_{j}.\end{equation} \tag{ S14 }$$

Note

$$\begin{equation}{\alpha }^{2}\langle \vert {\mathcal{M}}_{j}{\vert }^{2}\rangle ={\alpha }^{2}\left[{\left(\frac{\alpha \nu }{2}\right)}^{2}+\nu \right]={\mathrm{S}\mathrm{N}\mathrm{R}}^{2}+2\mathrm{S}\mathrm{N}\mathrm{R},\end{equation} \tag{ S15 }$$

and

$$\begin{equation}{\alpha }^{2}\langle {\mathcal{M}}_{j}{M}_{k}^{{\ast}}\rangle ={\alpha }^{2}{\left(\frac{\alpha \nu }{2}\right)}^{2}{\mathrm{e}}^{-D\vert j-k\vert {\Delta}t}.\end{equation} \tag{ S16 }$$

We also note the expression for the double-geometric series with 0 < ξ < 1

$$\begin{align}\hfill \sum\limits _{j\ne k}^{m}{\xi }^{\vert j-k\vert }& =2\sum\limits _{j=1}^{m-1}\sum\limits _{k=j+1}^{m}{\xi }^{k-j}\hfill \\ \hfill & =2\sum\limits _{j=1}^{m-1}\xi \frac{1-{\xi }^{m-j}}{1-\xi }\hfill \\ \hfill & =2\left(m-1\right)\frac{\xi }{1-\xi }-2\frac{{\xi }^{2}}{{\left(1-\xi \right)}^{2}}\left(1-{\xi }^{m-1}\right).\hfill \end{align} \tag{ S17 }$$

Thus,

$$\begin{align}\hfill \langle {\kappa }_{m}^{2}\rangle & ={\alpha }^{2}\sum\limits _{j=1}^{m}\langle \vert {\mathcal{M}}_{j}{\vert }^{2}\rangle +{\alpha }^{2}\sum\limits _{j\ne k}\langle {\mathcal{M}}_{j}{\mathcal{M}}_{k}^{{\ast}}\rangle \hfill \\ \hfill & =m\left({\mathrm{S}\mathrm{N}\mathrm{R}}^{2}+2\;\mathrm{S}\mathrm{N}\mathrm{R}\right)+{\mathrm{S}\mathrm{N}\mathrm{R}}^{2}\sum\limits _{j\ne k}{\xi }^{\vert j-k\vert },\quad \xi ={\mathrm{e}}^{-D{\Delta}t}\hfill \\ \hfill & =2m\;\mathrm{S}\mathrm{N}\mathrm{R}+\left(m\frac{{\mathrm{e}}^{D{\Delta}t}+1}{{\mathrm{e}}^{D{\Delta}t}-1}-2\enspace {\mathrm{e}}^{D{\Delta}t}\frac{1-{\mathrm{e}}^{-mD{\Delta}t}}{{\left({\mathrm{e}}^{D{\Delta}t}-1\right)}^{2}}\right){\mathrm{S}\mathrm{N}\mathrm{R}}^{2}.\hfill \end{align} \tag{ S18 }$$

For long times, t_m ≫ τ, this expression scales with m, as $\langle {\kappa }_{m}^{2}\rangle \approx 2{t}_{m}\tau {r}^{2}$. This follows from exp(DΔt) ≈ 1 + DΔt and exp(−mDΔt) ≈ 0 to leading order in DΔt. For short times, t_m ≪ τ, however, $\langle {\kappa }_{m}^{2}\rangle \approx {m}^{2}{\mathrm{S}\mathrm{N}\mathrm{R}}^{2}$. Specifically, a Fourier expansion of exp(−mDΔt) ≈ 1 − mDΔt + (mDΔt)²/2 to second order yields $\langle {\kappa }_{m}^{2}\rangle =2m\;\mathrm{S}\mathrm{N}\mathrm{R}+\left[{m}^{2}+\left(m-{m}^{3}\right)D{\Delta}t/3\right]{\mathrm{S}\mathrm{N}\mathrm{R}}^{2}+\mathcal{O}{\left(D{\Delta}t\right)}^{2}$.

We consider a minimal example for the estimation of a time-varying signal to illustrate the origin of the square-root scaling of the optimal effective measurement time with the fluctuation time-scale of the signal. Suppose a signal s*(t) varies stochastically with characteristic fluctuation time-scale τ, e.g. as a diffusion process with ds*/dt = ξ(t), where ξ(t) denotes Gaussian white noise with ⟨ξ(t)ξ(t')⟩ = τ⁻¹ δ(t − t'). The agent estimates $\hat{s}\left(t\right)={T}_{\mathrm{e}\mathrm{ff}}^{-1}{\int }_{t-{T}_{\mathrm{e}\mathrm{f}\mathrm{f}}}^{t} \enspace \mathrm{d}{t}^{\prime }\enspace \left({s}^{{\ast}}\left({t}^{\prime }\right)+\eta \left({t}^{\prime }\right)\right)$ by averaging noisy measurements over a time interval T_eff with measurement noise η(t) satisfying ⟨η(t)η(t')⟩ = r⁻¹ δ(t − t'). The variance σ² of the expected estimation error $\delta \left(t\right)=\hat{s}\left(t\right)-{s}^{{\ast}}\left(t\right)$ between the estimated value $\hat{s}\left(t\right)$ of the signal and its true value s*(t) can then be composed into two parts, where $\bar{{s}^{{\ast}}}={T}_{\mathrm{e}\mathrm{ff}}^{-1}{\int }_{t-{T}_{\mathrm{e}\mathrm{f}\mathrm{f}}}^{t} \enspace \mathrm{d}{t}^{\prime }\enspace {s}^{{\ast}}\left({t}^{\prime }\right)$

$$\begin{equation}{\sigma }^{2}=\langle {\left[\enspace \hat{s}\left(t\right)-{s}^{{\ast}}\left(t\right)\right]}^{2}\rangle \approx \underset{\begin{matrix}\hfill {\sigma }_{\text{sensing}}^{2}=1/\left(r{T}_{\mathrm{e}\mathrm{ff}}\right)\hfill \\ \hfill \;\text{sensing}\;\text{noise}\hfill \end{matrix}}{\underbrace{\left\langle {\left[\hat{s}\left(t\right)-\bar{{s}^{{\ast}}}\right]}^{2}\right\rangle }}+\underset{\begin{matrix}\hfill {\sigma }_{\text{time}}^{2}={T}_{\mathrm{e}\mathrm{ff}}/\left(3\tau \right)\hfill \\ \hfill \;\text{time}-\text{varying}\;\text{signal}\hfill \end{matrix}}{\underbrace{\left\langle {\left[{s}^{{\ast}}\left(t\right)-\bar{{s}^{{\ast}}}\left(t\right)\right]}^{2}\right\rangle }}.\end{equation} \tag{ S19 }$$

An explicit calculation yields

$$\begin{equation}{\sigma }^{2}={\left(r{T}_{\mathrm{e}\mathrm{ff}}\right)}^{-1}+{T}_{\mathrm{e}\mathrm{ff}}/\left(3\tau \right).\end{equation} \tag{ S20 }$$

Solving for dσ²/dT_eff = 0 yields T_eff = (3τ/r)^1/2.