This paper was converted on www.awesomepapers.org from LaTeX by an anonymous user.
Want to know more? Visit the Converter page.

Also at ]Universal Biology Institute, The University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, Tokyo 113-8654, Japan. Also at ]1Department of Mathematical Informatics, the Graduate School of Information Science and Technology, the University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, Tokyo, 113-8654, Japan

Gradient sensing limit of a cell when controlling the elongating direction

Kento Nakamura RIKEN Center for Brain Science, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan    Tetsuya J. Kobayashi tetsuya@mail.crmind.net Institute of Industrial Science, The University of Tokyo, 4-6-1, Komaba, Meguro-ku, Tokyo 153-8505 Japan [ [
Abstract

Eukaryotic cells perform chemotaxis by determining the direction of chemical gradients based on stochastic sensing of concentrations at the cell surface. To examine the efficiency of this process, previous studies have investigated the limit of estimation accuracy for gradients. However, these studies assume that the cell shape and gradient are constant, and do not consider how adaptive regulation of cell shape affects the estimation limit. Dynamics of cell shape during gradient sensing is biologically ubiquitous and can influence the estimation by altering the way the concentration is measured, and cells may strategically regulate their shape to improve estimation accuracy. To address this gap, we investigate the estimation limits in dynamic situations where cells change shape adaptively depending on the sensed signal. We approach this problem by analyzing the stationary solution of the Bayesian nonlinear filtering equation. By applying diffusion approximation to the ligand-receptor binding process and the Laplace method for the posterior expectation under a high signal-to-noise ratio regime, we obtain an analytical expression for the estimation limit. This expression indicates that estimation accuracy can be improved by elongating perpendicular to the estimated direction, which is also confirmed by numerical simulations. Our analysis provides a basis for clarifying the interplay between estimation and control in gradient sensing and sheds light on how cells optimize their shape to enhance chemotactic efficiency.

preprint: APS/123-QED

I Introduction

Cells have evolved sophisticated mechanisms to sense and adapt to their fluctuating environment. An example of such adaptation is eukaryotic chemotaxis, manifested by the directed movement of cells along chemical gradients [1]. In eukaryotic cells, chemotaxis relies on the spatial sensing of chemical gradients using an array of receptors distributed across the cell surface. The reliable detection of chemical gradients is challenging due to the inherent stochasticity of receptor-ligand binding [2, 3]. This stochasticity introduces measurement noise in the local ligand concentration and leads to errors in estimating spatial difference of the concentrations. Despite this uncertainty, eukaryotic cells exhibit a remarkable ability to detect and navigate shallow gradients [4, 5, 6, 7], raising the question of how sensory systems are optimized to process noisy chemical information. To address this question, theoretical studies have sought to elucidate the limits on the accuracy of chemical gradient sensing in eukaryotic cells set by stochastic uncertainty [3, 8, 9, 10, 11, 12, 13]. These studies have employed tools from estimation theory, information theory, and statistical physics to derive bounds on the precision of gradient sensing.

One notable factor that can influence the accuracy of gradient sensing is the shape of the cell. During chemotaxis, eukaryotic cells undergo significant shape changes, such as elongation, pseudopod extension, and polarization [14, 15, 16]. These shape changes modulate the spatial distribution of receptors on the cell surface, altering the way in which the cell samples the chemical gradients. For instance, an elongated cell may be more sensitive to concentration differences along its long axis but less sensitive to differences along its short axis. By dynamically regulating its shape, a eukaryotic cell may enhance its gradient sensing performance by sampling gradient information more effectively.

However, the theoretical understanding of how cell shape affects gradient sensing remains limited. Previous studies have investigated the impact of cell elongation on the static gradient sensing based on the optimal estimation theory [17, 18]. These studies have revealed that the estimation accuracy of gradient depends on the spatial configuration between the cell and the true gradient, such as the angle between the cell’s elongation axis and the gradient direction. However, these studies assume that the elongation direction is temporally fixed and externally determined, not actively adjusted by the cell in response to its environment. This assumption limits the applicability of their findings to actual scenarios where cells continuously adjust their elongation direction based on sensory inputs. Another line of studies has addressed this limitation by developing a phenomenological model that couples cell shape dynamics with an internal polarity mechanism that represents the cell’s estimate of the gradient direction [19, 14]. While these models provide valuable insights into the interplay between shape and sensing during chemotaxis, it does not provide a framework for deriving estimation limit attained by dynamic cell shape control.

In this work, we bridge this gap by developing optimality theory to derive the gradient sensing limit under the feedback control of cell shape. We frame the gradient sensing problem as a Bayesian nonlinear filtering problem [20, 21], where the cell estimates the gradient direction from noisy receptor activity measurements while controlling its shape. The Bayesian nonlinear filtering theory has been applied to derive the optimal cellular sensing strategy under fluctuating environments [22, 23, 24, 25]. The theory has recently been applied to the gradient sensing by a eukaryotic cell with a circular shape [13, 26], but has not been used under feedback control of sensory apparatus. We show that the Bayesian nonlinear filtering theory can be combined with feedback control strategy of sensory processes.

Our paper is organized as follows. In Section II, we introduce a mathematical model of gradient sensing by a cell with a dynamically controlled shape. This model incorporates the stochastic dynamics of receptor-ligand binding and the feedback control of cell shape based on the sensory information. In Section III, we formulate the gradient sensing problem as a Bayesian nonlinear filtering problem and derive the optimal filtering equation for the posterior distribution of the gradient direction, which sets a limit on the accuracy of gradient sensing. In Section IV, we specialize our analysis to the case where the cell can control its elongation direction to be either parallel or perpendicular to the estimated gradient direction. In this case, we derive an analytical expression for the asymptotic estimation error in the limit of high signal-to-noise ratio, and we validate our theory using numerical simulations. Finally, in section V, we conclude this work by discussing the possible directions for future research.

II Model

II.1 Setting

We formulate a spatial gradient sensing model that incorporates cell-shape modulation (Figure 1).

Refer to caption
Figure 1: (A)Schematic diagram of gradient sensing under feedback control of cell shape. The cell shape is controlled based on the estimate of the gradient direction and, in turn, affects the gradient estimation by biasing the ligand binding process. (B)Time lapse of a sample process of the gradient sensing. Blue and red arrows indicate the true and estimated gradient direction, respectively. The length of the red arrows indicate the confidence of estimate |𝐞^t||\hat{\bm{\mathrm{e}}}_{t}|. Orange dots show the locations of receptors at which ligand bound during time [t,t+Δt)[t,t+\Delta t) with Δt\Delta t being the discretization width. (C)Time evolution of variables in the sample process shown in (B). The first row shows the true gradient direction θt\theta_{t} (blue curve) and the angular locations of receptors at which ligand binds dYt\mathrm{d}Y_{t} (orange dots). The second row shows the posterior density of the gradient direction ZtZ_{t} and the estimated gradient direction θ^t\hat{\theta}_{t} (red line). The third row shows the true (blue curve) and estimated (red curve) gradient direction. The final row shows the error between the true and estimated direction measured by (θtθ^t)2/2(\theta_{t}-\hat{\theta}_{t})^{2}/2 (black curve) and the uncertainty of estimate represented by 1|𝐞^t|1-|\hat{\bm{\mathrm{e}}}_{t}| (red curve). The trajectories were simulated under the signal-to-noise ratio Λ=101\Lambda=10^{1} with the following dimensionless parameter values: λ=2103,α=101\lambda=2\cdot 10^{3},\alpha=10^{-1}, and σϵ=0.4\sigma\epsilon=-0.4.

We consider a cell placed at the origin of a 2-dimensional space with a linear concentration gradient. The gradient direction is denoted by θt[0,2π)\theta_{t}\in[0,2\pi) and is assumed to fluctuate over time following a circular diffusion process

dθt=2DdWt,\displaystyle\mathrm{d}\theta_{t}=\sqrt{2D}\mathrm{d}W_{t}, (1)

where WtW_{t} is a standard Wiener process and DD is the rotational diffusion coefficient. The ligand concentration c(𝒙,θ)c(\bm{x},\theta) at a given position 𝒙2\bm{x}\in\mathbb{R}^{2} when the gradient direction is θt=θ\theta_{t}=\theta is given by

c(𝒙,θ)=c0+c𝒙,𝐞(θ).\displaystyle c(\bm{x},\theta)=c_{0}+\nabla c\left\langle\bm{x},\bm{\mathrm{e}}(\theta)\right\rangle. (2)

Here, c0c_{0} and c\nabla c are constants representing the concentration at the origin and the gradient steepness, respectively, and 𝐞(θ)=[cos(θ),sin(θ)]𝕋\bm{\mathrm{e}}(\theta)=[\cos(\theta),\sin(\theta)]^{\mathbb{T}} is a unit vector pointing in the direction θ\theta.

The cell senses the gradient using receptors on its surface, represented as a closed curve rt={rtϕϕ[0,2π)}r_{t}=\{r_{t}^{\phi}\mid\phi\in[0,2\pi)\} based on the Fourier series:

rtϕ=r0(1+k=2ϵtkcos(k(ϕϕtk)))\displaystyle r_{t}^{\phi}=r_{0}\left(1+\sum_{k=2}^{\infty}\epsilon_{t}^{k}\cos(k(\phi-\phi_{t}^{k}))\right) (3)

where ϵtk\epsilon_{t}^{k}\in\mathbb{R} and ϕtk[0,2π)\phi_{t}^{k}\in[0,2\pi) are the amplitude and phase of the kk-th Fourier component, respectively. The first Fourier component is removed to ensure that the cell center 𝒟ϕrtϕ\oint\mathcal{D}\phi r_{t}^{\phi} coincides with the coordinate origin. Given the cell radius rtϕr_{t}^{\phi}, we assume that NN\in\mathbb{N} receptors are located on the cell surface at equal arc length intervals. The direction of the ii-th receptor is denoted by ϕi\phi_{i}, and the ligand binds to the receptor with a rate proportional to the local ligand concentration c(rtϕi𝐞(ϕi),θ)c(r_{t}^{\phi_{i}}\bm{\mathrm{e}}(\phi_{i}),\theta) at the receptor’s position rtϕi𝐞(ϕi)r_{t}^{\phi_{i}}\bm{\mathrm{e}}(\phi_{i}). We assume that the unbinding time is negligible compared to the time binding time. The number of ligand binding events at the ii-th receptor up to time tt is described by the Poisson process YtϕiY_{t}^{\phi_{i}} with rate process

λtϕi=λϕi(rt,θt),λϕi(r,θ)\displaystyle\lambda_{t}^{\phi_{i}}=\lambda^{\phi_{i}}(r_{t},\theta_{t}),\ \lambda^{\phi_{i}}(r,\theta) :=λ0Nc(rϕi𝐞(ϕi),θ).\displaystyle:=\frac{\lambda_{0}}{N}c(r^{\phi_{i}}\bm{\mathrm{e}}(\phi_{i}),\theta). (4)

where λ0/N\lambda_{0}/N is the binding rate of a single receptor under unit concentration.

Finally, we assume that the cell controls its shape rtr_{t} based on the ligand binding history measured by all receptors. By introducing the combined notation for the ligand binding process over receptors, Yt:={Ytϕii=1,,N}Y_{t}:=\{Y_{t}^{\phi_{i}}\mid i=1,\ldots,N\}, the feedback control assumption is expressed by setting ϵtk\epsilon_{t}^{k} and ϕtk\phi_{t}^{k} as functionals of the history of YtY_{t}:

ϵtk=ϵtk[Y0:t],ϕtk=ϕtk[Y0:t].\displaystyle\epsilon_{t}^{k}=\epsilon_{t}^{k}[Y_{0:t^{-}}],\ \phi_{t}^{k}=\phi_{t}^{k}[Y_{0:t^{-}}]. (5)

Here, Y0:t={Yss[0,t)}Y_{0:t^{-}}=\{Y_{s}\mid s\in[0,t)\} represents the history of observation process up to, but not including, time tt. The exclusion of YtY_{t} ensures causality between sensing and control; otherwise, the value of YY after a jump at time tt would affect whether the jump occurs at that time. In Section IV, we will consider a more concrete setting for feedback control of the cell shape to derive the analytical result.

II.2 Dimensionless system

To narrow down the parameter space without losing generality, we consider a dimensionless system on the above setting. By rescaling time as t~:=t/τD\tilde{t}:=t/\tau_{D} with τD:=1/D\tau_{D}:=1/D, we obtain the dimensionless dynamics of the gradient direction:

dθt~\displaystyle\mathrm{d}\theta_{\tilde{t}} =2dWt~,\displaystyle=\sqrt{2}\mathrm{d}W_{\tilde{t}}, (6)

By introducing several dimensionless parameters, the ligand-receptor binding process is obtained as the Poisson process Yt~ϕiY_{\tilde{t}}^{\phi_{i}} with the following dimensionless rate process:

λ~t~ϕi(θ)\displaystyle\tilde{\lambda}_{\tilde{t}}^{\phi_{i}}(\theta) :=λt~τDϕi(θ)τD=λ¯N(1+αr~t~ϕi𝐞(ϕi),𝐞(θ)).\displaystyle:=\lambda_{\tilde{t}\tau_{D}}^{\phi_{i}}(\theta)\tau_{D}=\frac{\bar{\lambda}}{N}\left(1+\alpha\tilde{r}_{\tilde{t}}^{\phi_{i}}\left\langle\bm{\mathrm{e}}(\phi_{i}),\bm{\mathrm{e}}(\theta)\right\rangle\right). (7)

Here, r~t~ϕi:=rt~τDϕi/r0\tilde{r}_{\tilde{t}}^{\phi_{i}}:=r_{\tilde{t}\tau_{D}}^{\phi_{i}}/r_{0} is the normalized cell radius, and λ¯:=λ0c0τD\bar{\lambda}:=\lambda_{0}c_{0}\tau_{D} and α:=r0c/c0\alpha:=r_{0}\nabla c/c_{0} are dimensionless parameters. λ¯\bar{\lambda} represents the binding number during the diffusion time τD\tau_{D} in the absence of a gradient, and α\alpha represents the relative difference in ligand concentration across the cell radius.

In the following sections, we consider the dimensionless system and omit the tildes from t~,λ~,\tilde{t},\tilde{\lambda}, and r~\tilde{r} for visual clarity.

III Estimation and control by Bayesian nonlinear filtering

We formulate a sensing problem for the gradient direction θt\theta_{t} under the feedback control of the cell shape rtr_{t}, and describe the optimal estimator derived based on the Bayesian nonlinear filtering theory. We assume that the cell estimates the gradient direction θt\theta_{t} based on the sensing history Y0:t:={Yss[0,t]}Y_{0:t}:=\{Y_{s}\mid s\in[0,t]\}, i.e., the estimated direction is expressed as ϑt[Y0:t]\vartheta_{t}[Y_{0:t}] with a functional ϑt\vartheta_{t}.

The goal is to find the functional ϑt\vartheta_{t} that optimizes an estimation performance. As a measure of estimation performance for directional statistics, we consider the circular variance of the difference between the true and estimated gradient direction, as in the previous study [13]:

CV(θtϑt[Y0:t]):=1|𝔼θt,Y0:t[exp(i(θtϑt[Y0:t]))]|,\displaystyle\mathrm{CV}(\theta_{t}-\vartheta_{t}[Y_{0:t}]):=1-\left|\mathbb{E}^{\theta_{t},Y_{0:t}}\left[\exp\left(\mathrm{i}(\theta_{t}-\vartheta_{t}[Y_{0:t}])\right)\right]\right|, (8)

where i\mathrm{i} is an imaginary unit and 𝔼θt,Y0:t\mathbb{E}^{\theta_{t},Y_{0:t}} represents the expectation with respect to the joint distribution of θt\theta_{t} and Y0:tY_{0:t}. As described in Ref. [13], the circular variance relates to the square error as CV(θtϑt[Y0:t])𝔼[(θtϑt[Y0:t])2]/2\mathrm{CV}(\theta_{t}-\vartheta_{t}[Y_{0:t}])\approx\mathbb{E}[(\theta_{t}-\vartheta_{t}[Y_{0:t}])^{2}]/2 when the estimated direction concentrates around the true gradient direction. The optimal estimator θ^t\hat{\theta}_{t} that minimizes the estimation error is represented based on the posterior density function ZtZ_{t} of the gradient direction θt\theta_{t} conditioned on the sensing history Y0:tY_{0:t} (see Appendix A for derivation):

θ^t\displaystyle\hat{\theta}_{t} =ϑt[Y0:t]:=argminϑt[Y0:t]CV(θtϑt[Y0:t])=arg𝐞^t,\displaystyle=\vartheta_{t}^{\ast}[Y_{0:t}]:=\mathop{\rm arg~{}min}\limits_{\vartheta_{t}[Y_{0:t}]}\mathrm{CV}(\theta_{t}-\vartheta_{t}[Y_{0:t}])=\arg\hat{\bm{\mathrm{e}}}_{t}, (9)
𝐞^t\displaystyle\hat{\bm{\mathrm{e}}}_{t} :=𝔼θt[𝐞(θt)Y0:t]=𝒟θ𝐞(θ)Zt(θ),\displaystyle:=\mathbb{E}^{\theta_{t}}[\bm{\mathrm{e}}(\theta_{t})\mid Y_{0:t}]=\oint\mathcal{D}\theta\bm{\mathrm{e}}(\theta)Z_{t}(\theta), (10)

where Zt(θ):=(θt=θY0:t)Z_{t}(\theta):=\mathbb{P}\left(\theta_{t}=\theta\mid Y_{0:t}\right) is the posterior density function of θt\theta_{t} given Y0:tY_{0:t}, and 𝔼θt[Y0:t]\mathbb{E}^{\theta_{t}}[\cdot\mid Y_{0:t}] represents the expectation with respect to Zt(θ)Z_{t}(\theta). 𝐞^t\hat{\bm{\mathrm{e}}}_{t} is a vector whose argument represents the estimated direction θ^t\hat{\theta}_{t}, and its magnitude 𝐞^t1\|\hat{\bm{\mathrm{e}}}_{t}\|\leq 1, where \|\cdot\| denotes a Euclidean norm, represents the certainty of the estimate. The magnitude takes the minimum 𝐞^t=0\|\hat{\bm{\mathrm{e}}}_{t}\|=0 when the posterior density function is uniform, Zt(θ)=(2π)1Z_{t}(\theta)=(2\pi)^{-1}, and the maximum 𝐞^t=1\|\hat{\bm{\mathrm{e}}}_{t}\|=1 when the posterior density concentrates completely on the estimated direction, Zt(θ)=δ(θθ^t)Z_{t}(\theta)=\delta(\theta-\hat{\theta}_{t}).

Using the vector 𝐞^t\hat{\bm{\mathrm{e}}}_{t}, we can express the limit of estimation error attained by the optimal estimator:

minϑtCV(θtϑt[Y0:t])=1𝔼Y0:t[|𝐞^t|],\displaystyle\min_{\vartheta_{t}}\mathrm{CV}(\theta_{t}-\vartheta_{t}[Y_{0:t}])=1-\mathbb{E}^{Y_{0:t}}[|\hat{\bm{\mathrm{e}}}_{t}|], (11)

where 𝔼Y0:t\mathbb{E}^{Y_{0:t}} represents the expectation with respect to the marginal distribution of Y0:tY_{0:t} with θt\theta_{t} being averaged out.

To calculate the optimal estimator based on the posterior density function Zt(θ)Z_{t}(\theta), we can use the following stochastic partial differential equation known as the Kushner equation in the nonlinear filtering theory (for a pedagogical introduction, see Ref. [21]; for detailed derivations, see the references therein and Ref. [27]):

dZt(θ)\displaystyle\mathrm{d}Z_{t}(\theta) =dpredZt(θ)+dobsZt(θ),\displaystyle=\mathrm{d}^{\mathrm{pred}}Z_{t}(\theta)+\mathrm{d}^{\mathrm{obs}}Z_{t}(\theta), (12)
dpredZt(θ)\displaystyle\mathrm{d}^{\mathrm{pred}}Z_{t}(\theta) =2Ztθ2(θ)dt,\displaystyle=\frac{\partial^{2}Z_{t}}{\partial\theta^{2}}(\theta)\mathrm{d}t, (13)
dobsZt(θ)=Zt(θ)i=1Nλtϕi(θ)λ^tϕiλ^tϕi(dYtϕiλ^tϕidt).\displaystyle\mathrm{d}^{\mathrm{obs}}Z_{t}(\theta)=Z_{t^{-}}(\theta)\sum_{i=1}^{N}\frac{\lambda_{t^{-}}^{\phi_{i}}(\theta)-\hat{\lambda}_{t^{-}}^{\phi_{i}}}{\hat{\lambda}_{t^{-}}^{\phi_{i}}}\left(\mathrm{d}Y_{t}^{\phi_{i}}-\hat{\lambda}_{t^{-}}^{\phi_{i}}\mathrm{d}t\right). (14)

Here, ZtZ_{t^{-}} and λ^tϕi:=𝔼θt[λtϕi(θt)Y0:t]\hat{\lambda}_{t^{-}}^{\phi_{i}}:=\mathbb{E}^{\theta_{t}}[\lambda_{t^{-}}^{\phi_{i}}(\theta_{t^{-}})\mid Y_{0:t^{-}}] with tt^{-} representing the time infinitesimally before tt are statistics conditioned on Y0:tY_{0:t^{-}}. This conditioning ensures that the coefficient of the jump term dYtϕi\mathrm{d}Y_{t}^{\phi_{i}} in the equation is independent of the value of YY at time tt. The first term, dpredZt(θ)\mathrm{d}^{\mathrm{pred}}Z_{t}(\theta), represents the prediction based on the prior knowledge about the dynamics of θt\theta_{t} and the second term, dobsZt(θ)\mathrm{d}^{\mathrm{obs}}Z_{t}(\theta), represents the correction based on the sensed signal dYt\mathrm{d}Y_{t}.

IV Analytical derivation of estimation limit

In this section, we derive an analytical expression for the estimation limit by approximating the stationary solution of the filtering equation. We introduce additional assumptions to make the problem analytically tractable.

IV.1 Additional assumptions for an analytical treatment

Refer to caption
Figure 2: Parameters of controlled cell elongation. (A) σ\sigma dependency of the relative angle between elongation direction (gray curve indicates an elongating cell) and the estimated direction (red arrow). A cell elongates along the estimated direction if σ=+1\sigma=+1 (upper panel) and elongates perpendicularly to the estimated direction if σ=1\sigma=-1 (lower panel). (B) ϵ\epsilon dependency of the extent of elongation for an elliptic shape (left) and a shape described by the second Fourier component (right).

To derive the estimation limit analytically, we focus on the control pattern of cell shape such that the cell elongates either parallel or perpendicular to the estimated direction θ^t\hat{\theta}_{t} (Figure 2A). We model this control pattern by imposing a condition on the second mode of the cell shape:

ϵt,2[Y0:t]=σϵ,ϕt,2[Y0:t]=θ^t=ϑt[Y0:t],\displaystyle\epsilon_{t,2}[Y_{0:t^{-}}]=\sigma\epsilon,\ \phi_{t,2}[Y_{0:t^{-}}]=\hat{\theta}_{t^{-}}=\vartheta_{t^{-}}^{\ast}[Y_{0:t^{-}}], (15)

where ϵ0\epsilon\geq 0 and σ{±1}\sigma\in\{\pm 1\} are constants. When σ=+1\sigma=+1 and σ=1\sigma=-1, the cell elongates parallel and perpendicular to the estimated direction, respectively (Figure 2A). In addition, we consider the first order perturbation of the cell shape around a circular one by assuming ϵt,k=O(ϵ)\epsilon_{t,k}=O(\epsilon) for every kk and ignoring the second and higher order terms, O(ϵ2)O(\epsilon^{2}).

We further restrict the parameter region of interest for analytical tractability. We assume that the ligand gradient is shallow so that the difference in ligand concentration across the spatial scale of the cell αrtϕi\alpha r_{t}^{\phi_{i}} is sufficiently small. Second, we assume that there are a sufficiently large number of receptors, N1N\gg 1. Finally, we assume that Λ:=λ¯α2/2\Lambda:=\bar{\lambda}\alpha^{2}/2 is sufficiently large. This assumption indicates that the cell is provided with sufficient amount of gradient information because the normalized ligand binding rate λ¯\bar{\lambda} and gradient steepness α\alpha are high. We note that Λ\Lambda is a dimensionless version of the “information rate” defined in the previous study [13].

IV.2 Derivation of estimation limit

The filtering equation is generally difficult to analyze because the posterior density function Zt(θ)Z_{t}(\theta) is not closed as a finite-dimensional model except in some special cases. To simplify the problem, we first discuss the behavior of the observation term and then combine the effect of the prediction term.

IV.2.1 Coordinate transform of the observation term

We can see from Eq. (14) that the observation term is proportional to Zt(θ)Z_{t}(\theta) except for λ^tϕi\hat{\lambda}_{t}^{\phi_{i}}, suggesting that Zt(θ)Z_{t}(\theta) may grow exponentially. For analyzing an exponentially growing variable, we employ the coordinate transformation to the logarithm of the original variable.

To easily transform the dynamics of Zt(θ)Z_{t}(\theta) to that of the log-posterior density logZt(θ)\log Z_{t}(\theta), we apply a diffusion approximation of dYtϕi\mathrm{d}Y_{t}^{\phi_{i}} in advance based on the weak gradient assumption. By ignoring o((αrtϕi)2)o((\alpha r_{t}^{\phi_{i}})^{2}) terms in the corresponding Fokker-Planck equation, we obtain:

dobsZt(θ)\displaystyle\mathrm{d}^{\mathrm{obs}}Z_{t}(\theta) Zt(θ)i=1N2ΛNrtϕi𝐞(ϕi),𝐞(θ)𝐞^tdWti,\displaystyle\approx Z_{t}(\theta)\sum_{i=1}^{N}\sqrt{\frac{2\Lambda}{N}}r_{t}^{\phi_{i}}\left\langle\bm{\mathrm{e}}(\phi_{i}),\bm{\mathrm{e}}(\theta)-\hat{\bm{\mathrm{e}}}_{t}\right\rangle\mathrm{d}W_{t}^{i},\quad (16)

where dWti\mathrm{d}W_{t}^{i} with i=1,,Ni=1,\ldots,N are the mutually independent standard Wiener processes (see Appendix B for derivation). Here, we have also taken the expectation over the true gradient direction θt\theta_{t} because we only need the statistics of 𝐞^t\hat{\bm{\mathrm{e}}}_{t} marginalized over θt\theta_{t} to obtain the estimation limit (Eq. (11)).

By applying the Ito’s formula to Eq. (16), we get the dynamics of the log-posterior density:

dobslogZt(θ)\displaystyle\mathrm{d}^{\mathrm{obs}}\log Z_{t}(\theta) Zt(θ)1dobsZt(θ)12Zt(θ)2dobs[Z(θ)]t\displaystyle\approx Z_{t}(\theta)^{-1}\mathrm{d}^{\mathrm{obs}}Z_{t}(\theta)-\frac{1}{2}Z_{t}(\theta)^{-2}\mathrm{d}^{\mathrm{obs}}[Z(\theta)]_{t} (17)
i=1N2ΛNrtϕi𝐞(ϕi),𝐞(θ)𝐞^tdWtii=1NΛN(rtϕi)2𝐞(ϕi),𝐞(θ)𝐞^t2dt\displaystyle\approx\sum_{i=1}^{N}\sqrt{\frac{2\Lambda}{N}}r_{t}^{\phi_{i}}\left\langle\bm{\mathrm{e}}(\phi_{i}),\bm{\mathrm{e}}(\theta)-\hat{\bm{\mathrm{e}}}_{t}\right\rangle\mathrm{d}W_{t}^{i}-\sum_{i=1}^{N}\frac{\Lambda}{N}(r_{t}^{\phi_{i}})^{2}\left\langle\bm{\mathrm{e}}(\phi_{i}),\bm{\mathrm{e}}(\theta)-\hat{\bm{\mathrm{e}}}_{t}\right\rangle^{2}\mathrm{d}t (18)
=i=1NΛN(rtϕi)2𝐞(ϕi),𝐞(θ)𝐞^t2dt+O(Λ1/2),\displaystyle=-\sum_{i=1}^{N}\frac{\Lambda}{N}(r_{t}^{\phi_{i}})^{2}\left\langle\bm{\mathrm{e}}(\phi_{i}),\bm{\mathrm{e}}(\theta)-\hat{\bm{\mathrm{e}}}_{t}\right\rangle^{2}\mathrm{d}t+O(\Lambda^{1/2}), (19)

where dobs[Z(θ)]t\mathrm{d}^{\mathrm{obs}}[Z(\theta)]_{t} is the quadratic variation process. We neglected the fluctuating part of the observation term (the first term on the second line of Eq. (19)) based on the large Λ\Lambda assumption. We note that the variance of the fluctuating part does not diverge in the limit of NN\rightarrow\infty because dWti\mathrm{d}W_{t}^{i} is independent for each ii.

When viewed in the logarithmic coordinate (Eq. (19)), Zt(θ)Z_{t}(\theta) dependent factors disappear from the observation term. The dependence on θ\theta is now simply represented by 𝐞(ϕi),𝐞(θ)𝐞^t2\left\langle\bm{\mathrm{e}}(\phi_{i}),\bm{\mathrm{e}}(\theta)-\hat{\bm{\mathrm{e}}}_{t}\right\rangle^{2}, which is further expanded up to a θ\theta-independent constant:

𝐞(ϕi),𝐞(θ)𝐞^t2\displaystyle\left\langle\bm{\mathrm{e}}(\phi_{i}),\bm{\mathrm{e}}(\theta)-\hat{\bm{\mathrm{e}}}_{t}\right\rangle^{2} =(𝐞(ϕi),𝐞(θ)|𝐞^t|𝐞(ϕi),𝐞(θ^t))2\displaystyle=\left(\left\langle\bm{\mathrm{e}}(\phi_{i}),\bm{\mathrm{e}}(\theta)\right\rangle-|\hat{\bm{\mathrm{e}}}_{t}|\left\langle\bm{\mathrm{e}}(\phi_{i}),\bm{\mathrm{e}}(\hat{\theta}_{t})\right\rangle\right)^{2} (20)
=12𝐞(2ϕi),𝐞(2θ)+|𝐞^t|(𝐞(θ^t),𝐞(θ)𝐞(2ϕi),𝐞(θ+θ^t))+const.\displaystyle=\frac{1}{2}\left\langle\bm{\mathrm{e}}(2\phi_{i}),\bm{\mathrm{e}}(2\theta)\right\rangle+|\hat{\bm{\mathrm{e}}}_{t}|\left(\left\langle\bm{\mathrm{e}}(\hat{\theta}_{t}),\bm{\mathrm{e}}(\theta)\right\rangle-\left\langle\bm{\mathrm{e}}(2\phi_{i}),\bm{\mathrm{e}}(\theta+\hat{\theta}_{t})\right\rangle\right)+const. (21)
=12𝐞(2(ϕiθ^t)),𝐞(2(θθ^t))|𝐞^t|(𝐞(0)+𝐞(2(ϕiθ^t)),𝐞(θθ^t))+const.\displaystyle=\frac{1}{2}\left\langle\bm{\mathrm{e}}(2(\phi_{i}-\hat{\theta}_{t})),\bm{\mathrm{e}}(2(\theta-\hat{\theta}_{t}))\right\rangle-|\hat{\bm{\mathrm{e}}}_{t}|\left(\left\langle\bm{\mathrm{e}}(0)+\bm{\mathrm{e}}(2(\phi_{i}-\hat{\theta}_{t})),\bm{\mathrm{e}}(\theta-\hat{\theta}_{t})\right\rangle\right)+const. (22)

IV.2.2 Cell shape dependence of the observation term

We can summarize the dependence of the observation term on cell shape by substituting Eq. (22) into Eq. (19) and taking the sum over ii:

dobslogZt(θ)\displaystyle\mathrm{d}^{\mathrm{obs}}\log Z_{t}(\theta) Λ(12𝑹2,𝐞(2(θθ^t))|𝐞^t|𝑹0+𝑹2,𝐞(θθ^t))dt+O(Λ1/2).\displaystyle\approx-\Lambda\left(\frac{1}{2}\left\langle\bm{R}_{2},\bm{\mathrm{e}}(2(\theta-\hat{\theta}_{t}))\right\rangle-|\hat{\bm{\mathrm{e}}}_{t}|\left\langle\bm{R}_{0}+\bm{R}_{2},\bm{\mathrm{e}}(\theta-\hat{\theta}_{t})\right\rangle\right)\mathrm{d}t+O(\Lambda^{1/2}). (23)

Here, 𝑹k:=N1i=1N(rtϕi)2𝐞(k(ϕiθ^t))\bm{R}_{k}:=N^{-1}\sum_{i=1}^{N}(r_{t}^{\phi_{i}})^{2}\bm{\mathrm{e}}(k(\phi_{i}-\hat{\theta}_{t})) is a quantity that summarizes the cell shape and its configuration with respect to the estimated direction θ^t\hat{\theta}_{t}. Specifically, 𝑹0\bm{R}_{0} and 𝑹2\bm{R}_{2} represent the overall radius and elongation of the cell, respectively. In the limit of large NN and small ϵ\epsilon, 𝑹0\bm{R}_{0} and 𝑹2\bm{R}_{2} are approximated as follows:

𝑹0\displaystyle\bm{R}_{0} (𝒟ϕρtϕ)1𝒟ϕρtϕ(rtϕ)2𝐞(0)𝐞(0)+o(ϵ),\displaystyle\approx\left(\oint\mathcal{D}\phi\rho_{t}^{\phi}\right)^{-1}\oint\mathcal{D}\phi\rho_{t}^{\phi}(r_{t}^{\phi})^{2}\bm{\mathrm{e}}(0)\approx\bm{\mathrm{e}}(0)+o(\epsilon), (24)
𝑹2\displaystyle\bm{R}_{2} (𝒟ϕρtϕ)1𝒟ϕρtϕ(rtϕ)2𝐞(2(ϕθ^t))32σϵ𝐞(0)+o(ϵ).\displaystyle\approx\left(\oint\mathcal{D}\phi\rho_{t}^{\phi}\right)^{-1}\oint\mathcal{D}\phi\rho_{t}^{\phi}(r_{t}^{\phi})^{2}\bm{\mathrm{e}}(2(\phi-\hat{\theta}_{t}))\approx\frac{3}{2}\sigma\epsilon\bm{\mathrm{e}}(0)+o(\epsilon). (25)

Here, we defined ρtϕ\rho_{t}^{\phi} such that ρtϕ𝒟ϕ\rho_{t}^{\phi}\mathcal{D}\phi indicates the arc length of the cell circumference spanned by the angle interval [ϕ,ϕ+𝒟ϕ)[\phi,\phi+\mathcal{D}\phi) and used the approximation ρtϕ(rtϕ)21+3(rtϕ1)\rho_{t}^{\phi}(r_{t}^{\phi})^{2}\approx 1+3(r_{t}^{\phi}-1), which holds for sufficiently small ϵ\epsilon. The term ρtϕ\rho_{t}^{\phi} is introduced due to the assumption that receptors are distributed equally with respect to arc length intervals, rather than angular intervals. By substituting Eqs. (24) and (25) into Eq. (23), we get the following approximation of the observation term:

dobslogZt(θ)\displaystyle\mathrm{d}^{\mathrm{obs}}\log Z_{t}(\theta) Λ(34σϵcos(2(θθ^t))|𝐞^t|(1+32σϵ)cos(θθ^t))dt+O(Λ1/2)\displaystyle\approx-\Lambda\left(\frac{3}{4}\sigma\epsilon\cos(2(\theta-\hat{\theta}_{t}))-|\hat{\bm{\mathrm{e}}}_{t}|\left(1+\frac{3}{2}\sigma\epsilon\right)\cos(\theta-\hat{\theta}_{t})\right)\mathrm{d}t+O(\Lambda^{1/2}) (26)
=Λ(|𝐞^t|cos(θθ^t)+32σϵ(|𝐞^t|cos(θθ^t)12cos(2(θθ^t))))dt+O(Λ1/2).\displaystyle=\Lambda\left(|\hat{\bm{\mathrm{e}}}_{t}|\cos(\theta-\hat{\theta}_{t})+\frac{3}{2}\sigma\epsilon\left(|\hat{\bm{\mathrm{e}}}_{t}|\cos(\theta-\hat{\theta}_{t})-\frac{1}{2}\cos(2(\theta-\hat{\theta}_{t}))\right)\right)\mathrm{d}t+O(\Lambda^{1/2}). (27)

Eq. (27) indicates that, when ϵ\epsilon is sufficiently small, the cell shape affects the estimation through its second Fourier mode σϵ\sigma\epsilon and that ZtZ_{t} would concentrate around the estimated direction θ^t\hat{\theta}_{t} under large Λ\Lambda.

.

IV.2.3 Integration of the effect of the prediction term

We incorporate the effect of the prediction term by using the concentration property of Zt(θ)Z_{t}(\theta) around θ^t\hat{\theta}_{t} suggested by Eq. (27). By assuming that the posterior density is unimodal with a maximum point θtmax\theta_{t}^{\max}, which eventually coincides with θ^t\hat{\theta}_{t}, we expand logZt\log Z_{t} into a quadratic function:

logZt(θ)12κt(θθtmax)2+o((θθtmax)2)+const.\displaystyle\log Z_{t}(\theta)\approx-\frac{1}{2}\kappa_{t}(\theta-\theta_{t}^{\max})^{2}+o((\theta-\theta_{t}^{\max})^{2})+\mathrm{const.} (28)

where κt:=2θ2logZt(θtmax)>0\kappa_{t}:=-\frac{\partial^{2}}{\partial\theta^{2}}\log Z_{t}(\theta_{t}^{\max})>0 indicates the concentration of the posterior density. When Λ\Lambda is high, the posterior density would concentrate on a maximal point and κt\kappa_{t} would be large. For a large value of κt\kappa_{t}, we can use Laplace’s method to approximate the posterior expectation:

𝔼[F(θt)Y0:t]\displaystyle\mathbb{E}\left[F(\theta_{t})\mid Y_{0:t}\right] =𝒟θF(θ)exp(logZt(θ))\displaystyle=\oint\mathcal{D}\theta F(\theta)\exp(\log Z_{t}(\theta)) (29)
𝒟θF(θ)exp(12κt(θθtmax)2)\displaystyle\approx\oint\mathcal{D}\theta F(\theta)\exp\left(-\frac{1}{2}\kappa_{t}(\theta-\theta_{t}^{\max})^{2}\right) (30)
F(θtmax)+12F′′(θtmax)κt1+o(κt1).\displaystyle\approx F(\theta_{t}^{\max})+\frac{1}{2}F^{\prime\prime}(\theta_{t}^{\max})\kappa_{t}^{-1}+o(\kappa_{t}^{-1}). (31)

By substituting F(θ)=𝐞(θ)F(\theta)=\bm{\mathrm{e}}(\theta), we obtain the approximation of the circular variance in terms of κt\kappa_{t}:

𝐞^t(112κt1)𝐞(θtmax),CV[θtθ^t]12κt1.\displaystyle\hat{\bm{\mathrm{e}}}_{t}\approx\left(1-\frac{1}{2}\kappa_{t}^{-1}\right)\bm{\mathrm{e}}(\theta_{t}^{\max}),\ \mathrm{CV}[\theta_{t}-\hat{\theta}_{t}]\approx\frac{1}{2}\kappa_{t}^{-1}. (32)

Eq. (32) shows that the maximum point θtmax\theta_{t}^{\max} approximately coincides with the estimated direction θ^t:=arg(𝐞^t)\hat{\theta}_{t}:=\arg(\hat{\bm{\mathrm{e}}}_{t}). Therefore, the concentration κt\kappa_{t} of the posterior density around the estimated direction θ^t\hat{\theta}_{t} is sufficient for expressing the estimation limit.

We finally get the estimation limit by calculating the behavior of the observation and prediction terms around the estimated direction. The prediction term around the estimated direction is characterized by κt\kappa_{t}:

dpredlogZt(θ)\displaystyle\mathrm{d}^{\mathrm{pred}}\log Z_{t}(\theta) =[2logZt(θ)θ2+(logZt(θ)θ)2]dt\displaystyle=\left[\frac{\partial^{2}\log Z_{t}(\theta)}{\partial\theta^{2}}+\left(\frac{\partial\log Z_{t}(\theta)}{\partial\theta}\right)^{2}\right]\mathrm{d}t (33)
κt2(θθ^t)2dt.\displaystyle\approx\kappa_{t}^{2}(\theta-\hat{\theta}_{t})^{2}\mathrm{d}t. (34)

The observation term (Eq. (27)) is expanded around θ=θ^t\theta=\hat{\theta}_{t} as follows:

dobslogZt(θ)\displaystyle\mathrm{d}^{\mathrm{obs}}\log Z_{t}(\theta) 122θ2dobslogZt(θ^t)(θθ^t)2\displaystyle\approx-\frac{1}{2}\frac{\partial^{2}}{\partial\theta^{2}}\mathrm{d}^{\mathrm{obs}}\log Z_{t}(\hat{\theta}_{t})\cdot(\theta-\hat{\theta}_{t})^{2} (35)
12Λ(132σϵ)(θθ^t)2dt,\displaystyle\approx-\frac{1}{2}\Lambda\left(1-\frac{3}{2}\sigma\epsilon\right)(\theta-\hat{\theta}_{t})^{2}\mathrm{d}t, (36)

where we used the approximation |𝐞^t|1|\hat{\bm{\mathrm{e}}}_{t}|\approx 1 which holds under large κ\kappa with Eq. (32). By combining Eqs. (34) and (36), we obtain the stationary value of κt\kappa_{t} satisfying dlogZt(θ)=0\mathrm{d}\log Z_{t}(\theta)=0:

κtΛ2(134σϵ).\displaystyle\kappa_{t}\approx\sqrt{\frac{\Lambda}{2}}\left(1-\frac{3}{4}\sigma\epsilon\right). (37)

We note that this relation between κt\kappa_{t} and Λ\Lambda implies that the assumption of large κt\kappa_{t} under large Λ\Lambda is self-consistent. By substituting Eq. (37) into Eq. (32), we obtain the approximate estimation limit at the stationary state:

CV(θtθ^t)12Λ(1+34σϵ).\displaystyle\mathrm{CV}(\theta_{t}-\hat{\theta}_{t})\approx\frac{1}{\sqrt{2\Lambda}}\left(1+\frac{3}{4}\sigma\epsilon\right). (38)

This result indicates that the estimation is improved when σ<0\sigma<0, i.e., the cell elongates perpendicular to the estimated direction, and impaired when σ>0\sigma>0, i.e., the cell elongates along the estimated direction. We note that this estimation limit reproduces previous result derived for a circular cell under a time-discrete setting when ϵ=0\epsilon=0 [13].

IV.3 Numerical validation of estimation limit

Refer to caption
Figure 3: Estimation error and normalized estimation error as functions of the shape parameter σϵ\sigma\epsilon for different values of the signal-to-noise ratio Λ\Lambda. (A1, A2) Estimation error for an elliptic shape (A1) and a shape described by the second Fourier component (A2). Different curves represent different values of Λ\Lambda, which are obtained by varying the ligand binding rate λ¯\bar{\lambda} while keeping the gradient steepness α\alpha fixed. (B1, B2) Estimation error for an elliptic shape (B1) and a shape described by the second Fourier component (B2). Different curves represent different values of Λ\Lambda, which are obtained by varying the gradient steepness α\alpha while keeping the ligand binding rate λ¯\bar{\lambda} fixed. In all panels, solid lines represent theoretical predictions, while markers indicate simulation results.

To validate the theoretical prediction that the control of elongation direction can improve or impair the estimation performance, as indicated by Eq. (38), we perform numerical simulations of the gradient sensing process under different cell shape configurations (Figure 2(B)). We consider two types of cell shapes: a shape described by the second Fourier mode and an elliptic shape. The shape described by the second Fourier mode serves as a simple setting where only the second mode has non-zero values, and the elliptic shape is a widely used candidate for modeling elongating cells.

To ensure a fair comparison across elongation parameters, we rescale the radius as rtϕrtϕ/S(ϵ)/S(0)r_{t}^{\phi}\rightarrow r_{t}^{\phi}/\sqrt{S(\epsilon)/S(0)}, where S(ϵ)S(\epsilon) represents the area of the elongated shape with elongation parameter ϵ\epsilon, and S(0)S(0) represents the area of the circular shape (i.e., when ϵ=0\epsilon=0). This rescaling factor adjusts the size of the cell so that its area is independent of the elongation parameter ϵ\epsilon. The factor S(ϵ)/S(0)S(\epsilon)/S(0) deviates from unity only by o(ϵ)o(\epsilon), and consequently, the rescaling does not change the analytical expression of the estimation limit given by Eq. (24).

In our simulations, the gradient direction θt\theta_{t} evolves according to Eq. (6), while the ligand-receptor binding is modeled as a Poisson process with the rate given by Eq. (7). The posterior distribution of the gradient direction is computed by discretizing the filtering equation (Eq. (12)) over the direction θ\theta.

To quantify the estimation performance, we calculate the circular variance (Eq. (8)) by approximating the expectation using the Monte Carlo method with 10410^{4} samples. The initial state of the posterior density function is set to be uniform. We observe that the circular variance reaches a nearly stationary value after a time period of 2Λ1/22\Lambda^{-1/2} (see Figure C.4). Therefore, we define the stationary estimation limit as the average of the circular variance over the time interval t[2Λ1/2,3Λ1/2]t\in[2\Lambda^{-1/2},3\Lambda^{-1/2}].

Figure 3 presents the estimation limit for various values of the signal-to-noise ratio Λ\Lambda and the elongation parameter σϵ\sigma\epsilon. In the high signal-to-noise ratio regime, the results confirm our theoretical prediction: the estimation limit improves when the cell elongates perpendicular to the estimated direction (σ<0\sigma<0) and deteriorates when the cell elongates along the estimated direction (σ>0\sigma>0). This trend holds for both the shape described by the second Fourier mode and the elliptic shape. Moreover, the dependence of the estimation limit on the elongation magnitude ϵ\epsilon is well captured by the analytical approximation (Eq. (38)) for both cell shapes. This finding suggests that the second mode of the cell shape plays a dominant role in determining the estimation performance under linear concentration gradient.

Even when the signal-to-noise ratio decreases, the analytical formula explains the simulation results moderately well for Λ10\Lambda\geq 10. While simulation results start to deviate from the analytical prediction as Λ\Lambda decreases below 10, the qualitative result that the estimation improves for σ<0\sigma<0 remains valid even for moderate values of Λ=1\Lambda=1, where estimation becomes increasingly difficult (see Figure D.5 for sample trajectories of θt\theta_{t} and θ^t\hat{\theta}_{t}).

In the ellipse case with Λ=1\Lambda=1, a slight improvement in estimation is also observed for σ>0\sigma>0. This suggests the presence of additional behaviors that are not captured by our current analysis, which is reasonable given the qualitative differences in the low SNR regime.

Despite some deviations in the low SNR regime, these numerical results quantitatively support our theoretical analysis, demonstrating that the strategic control of cell shape elongation can indeed enhance or hinder the accuracy of gradient sensing, depending on the orientation of the elongation relative to the estimated gradient direction. Further investigation into the low SNR regime may reveal additional insights and will be a subject for future work.

V Discussion

In this study, we utilized Bayesian nonlinear filtering theory to derive the estimation limits for fluctuating gradient directions under the feedback control of cell shape. By applying the diffusion approximation to the stochastic ligand-receptor binding process and the Laplace approximation for posterior expectation, we obtained an analytical expression for the estimation limit valid in the high signal-to-noise ratio (SNR) regime. Our results align with previously identified limits for circular cell shapes in discrete-time settings [13].

Our key finding is that incorporating active cell shape control can enhance gradient sensing accuracy. We showed that when cells elongate perpendicular to the estimated gradient direction, the estimation error is reduced compared to the case of a static circular cell shape. This result parallels with previous theoretical findings that elongation perpendicular to the true gradient direction improves estimation performance [17]. Our work extends this previous result by considering the more realistic scenario where the elongation direction is determined by sensory information and fluctuates around the true gradient direction due to sensory noise. We demonstrate that even in this case, dynamic shape control can enhance direction estimation.

An important future direction is to identify scenarios where cell shape dynamics provides biologically more significant enhancement of gradient sensing accuracy. While our focus on the high SNR regime allowed us to obtain the analytical expression, the observed improvement of estimation limit implicates only marginal enhancement in chemotaxis ability. To clarify this point, we note that chemotaxis index, defined as the average migration along the gradient direction, CI=𝔼[cos(θtθ^t)]\mathrm{CI}=\mathbb{E}[\cos(\theta_{t}-\hat{\theta}_{t})], relates to the estimation error as CI=1CV[θtθ^t]\mathrm{CI}=1-\mathrm{CV}[\theta_{t}-\hat{\theta}_{t}], if we assume that cells migrate toward the estimated direction θ^t\hat{\theta}_{t} and the estimated direction fluctuates symmetrically around the true gradient 𝔼[sin(θtθ^t)]=0\mathbb{E}[\sin(\theta_{t}-\hat{\theta}_{t})]=0. Therefore, the same amount of reduction in the estimation error is more pronounced under the lower SNR regime in terms of the increasing rate of the chemotaxis index. We may find biologically significant improvement of gradient sensing ability by extensively investigating low SNR regime.

Another valuable extension of our analysis would be to account for uncertainty in the gradient steepness, especially in the low SNR regime. In our current study, we implicitly assumed that the gradient steepness varies on a slower timescale than the gradient direction, allowing us to treat the steepness as a known constant. However, in reality, the gradient steepness can change significantly during migration, and previous work has shown that uncertainty in the steepness can have a notable impact on the maximum likelihood estimate of the gradient direction [18]. Under the low SNR regime, the estimated direction is shown to be biased perpendicularly to or along with the elongation axis when the steepness is known or unknown, respectively. Our Bayesian filtering framework provides a natural way to incorporate steepness uncertainty by modeling the steepness as a stochastic process with a characteristic fluctuation timescale. By tuning this timescale, we can interpolate the two extreme cases of known and unknown steepness. This approach could reveal interesting transitions in the bias of the estimated direction and, consequently, the optimal sensing strategy as a function of the steepness uncertainty. Characterizing such transitions and their dependence on SNR and other model parameters is an intriguing direction for future research.

Acknowledgements

The first author received a JSPS Research Fellowship (Grant Number 20J21362). This research was supported by JSPS KAKENHI (Grant Number 19H05799 and 24H02148) and JST CREST (Grant Number JPMJCR2011).

References

  • Parent and Devreotes [1999] C. A. Parent and P. N. Devreotes, Science 284, 765 (1999).
  • Ueda et al. [2001] M. Ueda, Y. Sako, T. Tanaka, P. Devreotes, and T. Yanagida, Science 294, 864 (2001).
  • Ueda and Shibata [2007] M. Ueda and T. Shibata, Biophysical journal 93, 11 (2007).
  • Song et al. [2006] L. Song, S. M. Nadkarni, H. U. Bödeker, C. Beta, A. Bae, C. Franck, W.-J. Rappel, W. F. Loomis, and E. Bodenschatz, European journal of cell biology 85, 981 (2006).
  • Van Haastert and Postma [2007] P. J. Van Haastert and M. Postma, Biophysical journal 93, 1787 (2007).
  • Fuller et al. [2010] D. Fuller, W. Chen, M. Adler, A. Groisman, H. Levine, W.-J. Rappel, and W. F. Loomis, Proceedings of the National Academy of Sciences 107, 9656 (2010).
  • Amselem et al. [2012] G. Amselem, M. Theves, A. Bae, C. Beta, and E. Bodenschatz, Physical review letters 109, 108103 (2012).
  • Endres and Wingreen [2008] R. G. Endres and N. S. Wingreen, Proceedings of the National Academy of Sciences 105, 15749 (2008).
  • Rappel and Levine [2008] W.-J. Rappel and H. Levine, Physical review letters 100, 228101 (2008).
  • Hu et al. [2010] B. Hu, W. Chen, W.-J. Rappel, and H. Levine, Physical review letters 105, 048104 (2010).
  • Hu et al. [2011a] B. Hu, W. Chen, H. Levine, and W.-J. Rappel, Journal of statistical physics 142, 1167 (2011a).
  • Aquino et al. [2014] G. Aquino, L. Tweedy, D. Heinrich, and R. G. Endres, Scientific reports 4, 5688 (2014).
  • Novak and Friedrich [2021] M. Novak and B. M. Friedrich, New Journal of Physics 23, 043026 (2021).
  • Tweedy et al. [2013] L. Tweedy, B. Meier, J. Stephan, D. Heinrich, and R. G. Endres, Scientific reports 3, 1 (2013).
  • van Haastert et al. [2017] P. J. van Haastert, I. Keizer-Gunnink, and A. Kortholt, Molecular biology of the cell 28, 922 (2017).
  • Imoto et al. [2021] D. Imoto, N. Saito, A. Nakajima, G. Honda, M. Ishida, T. Sugita, S. Ishihara, K. Katagiri, C. Okimura, Y. Iwadate, et al., PLoS computational biology 17, e1009237 (2021).
  • Hu et al. [2011b] B. Hu, W. Chen, W.-J. Rappel, and H. Levine, Physical Review E 83, 021917 (2011b).
  • Baba et al. [2012] A. Baba, T. Hiraiwa, and T. Shibata, Physical Review E 86, 060901 (2012).
  • Hiraiwa et al. [2013] T. Hiraiwa, A. Baba, and T. Shibata, The European Physical Journal E 36, 1 (2013).
  • Bain and Crisan [2009] A. Bain and D. Crisan, Fundamentals of stochastic filtering, Vol. 3 (Springer, 2009).
  • Kutschireiter et al. [2020] A. Kutschireiter, S. C. Surace, and J.-P. Pfister, Journal of Mathematical Psychology 94, 102307 (2020).
  • Kobayashi [2010] T. J. Kobayashi, Physical review letters 104, 228104 (2010).
  • Zechner et al. [2016] C. Zechner, G. Seelig, M. Rullan, and M. Khammash, Proceedings of the National Academy of Sciences 113, 4729 (2016).
  • Mora and Nemenman [2019] T. Mora and I. Nemenman, Physical review letters 123, 198101 (2019).
  • Nakamura and Kobayashi [2021] K. Nakamura and T. J. Kobayashi, Physical Review Letters 126, 128102 (2021).
  • Auconi et al. [2022] A. Auconi, M. Novak, and B. M. Friedrich, Europhysics Letters 138, 12001 (2022).
  • Venugopal et al. [2015] M. Venugopal, R. M. Vasu, and D. Roy, IEEE Transactions on Automatic Control 61, 823 (2015).
  • Kanazawa and Sornette [2020] K. Kanazawa and D. Sornette, Physical Review Research 2, 033442 (2020).

Appendix A Estimation limit of directional statistics

We derive the equations for estimation limit expressed with Eqs. (9) and (11) by showing that the estimation error CV[θtϑt[Y0:t]]\mathrm{CV}[\theta_{t}-\vartheta_{t}[Y_{0:t}]] is bounded from below by 1𝔼Y0:t[𝐞^t]1-\mathbb{E}^{Y_{0:t}}[\hat{\bm{\mathrm{e}}}_{t}] and the bound is attained by ϑt[Y0:t]=θ^t\vartheta_{t}[Y_{0:t}]=\hat{\theta}_{t}. To do that, we note that 𝔼θt[exp(iθt)Y0:t]\mathbb{E}^{\theta_{t}}[\exp(i\theta_{t})\mid Y_{0:t}] on a complex plane can be idenfied with the real vector 𝐞^t:=𝔼θt[𝐞(θt)Y0:t]\hat{\bm{\mathrm{e}}}_{t}:=\mathbb{E}^{\theta_{t}}[\bm{\mathrm{e}}(\theta_{t})\mid Y_{0:t}] and thus 𝔼θt[exp(iθt)Y0:t]=𝐞^texp(iθ^t)\mathbb{E}^{\theta_{t}}[\exp(\mathrm{i}\theta_{t})\mid Y_{0:t}]=\|\hat{\bm{\mathrm{e}}}_{t}\|\exp(\mathrm{i}\hat{\theta}_{t}) holds. By using the property of conditional expectation, we can evaluate the estimation error as follows:

CV[θtϑt[Y0:t]]\displaystyle\mathrm{CV}\left[\theta_{t}-\vartheta_{t}[Y_{0:t}]\right] :=1|𝔼θt,Y0:t[exp(i(θtϑt[Y0:t]))]|\displaystyle:=1-\left|\mathbb{E}^{\theta_{t},Y_{0:t}}\left[\exp\left(\mathrm{i}(\theta_{t}-\vartheta_{t}[Y_{0:t}])\right)\right]\right| (39)
=1|𝔼Y0:t[𝔼θt[exp(iθt)Y0:t]exp(iϑt[Y0:t])]|\displaystyle=1-\left|\mathbb{E}^{Y_{0:t}}\left[\mathbb{E}^{\theta_{t}}\left[\exp\left(\mathrm{i}\theta_{t}\right)\mid Y_{0:t}\right]\exp\left(-\mathrm{i}\vartheta_{t}[Y_{0:t}]\right)\right]\right| (40)
=1|𝔼Y[𝐞^texp(iθ^t)exp(iϑt[Y0:t])]|\displaystyle=1-\left|\mathbb{E}_{Y}\left[\|\hat{\bm{\mathrm{e}}}_{t}\|\exp(\mathrm{i}\hat{\theta}_{t})\exp\left(-\mathrm{i}\vartheta_{t}[Y_{0:t}]\right)\right]\right| (41)
=1𝔼Y[𝐞^tcos(θ^tϑt[Y0:t])]2+𝔼Y[𝐞^tsin(θ^tϑt[Y0:t])]2\displaystyle=1-\sqrt{\mathbb{E}_{Y}\left[\|\hat{\bm{\mathrm{e}}}_{t}\|\cos(\hat{\theta}_{t}-\vartheta_{t}[Y_{0:t}])\right]^{2}+\mathbb{E}_{Y}\left[\|\hat{\bm{\mathrm{e}}}_{t}\|\sin(\hat{\theta}_{t}-\vartheta_{t}[Y_{0:t}])\right]^{2}} (42)
1𝔼Y[𝐞^t]2(𝔼Y[cos(θ^tϑt[Y0:t])]2+𝔼Y[sin(θ^tϑt[Y0:t])]2)\displaystyle\geq 1-\sqrt{\mathbb{E}_{Y}\left[\|\hat{\bm{\mathrm{e}}}_{t}\|\right]^{2}\left(\mathbb{E}_{Y}\left[\cos(\hat{\theta}_{t}-\vartheta_{t}[Y_{0:t}])\right]^{2}+\mathbb{E}_{Y}\left[\sin(\hat{\theta}_{t}-\vartheta_{t}[Y_{0:t}])\right]^{2}\right)} (43)
1𝔼Y[𝐞^t]2(𝔼Y[cos2(θ^tϑt[Y0:t])+sin2(θ^tϑt[Y0:t])])\displaystyle\geq 1-\sqrt{\mathbb{E}_{Y}\left[\|\hat{\bm{\mathrm{e}}}_{t}\|\right]^{2}\left(\mathbb{E}_{Y}\left[\cos^{2}(\hat{\theta}_{t}-\vartheta_{t}[Y_{0:t}])+\sin^{2}(\hat{\theta}_{t}-\vartheta_{t}[Y_{0:t}])\right]\right)} (44)
=1𝔼Y[𝐞^t].\displaystyle=1-\mathbb{E}_{Y}\left[\|\hat{\bm{\mathrm{e}}}_{t}\|\right]. (45)

Here, we used the Cauchy-Schwartz inequality in the fifth line and the Jensen inequality in the six-th line. We can check that θ^t\hat{\theta}_{t} attains the lower bound by substituting ϑt[Y0:t]=θ^t\vartheta_{t}[Y_{0:t}]=\hat{\theta}_{t} in the third line. Therefore, the optimal estimator and estimation limit are given by the formulae in Eq. (9) and Eq. (11), respectively.

Appendix B Filtering process marginalized over hidden process and diffusion approximation of Poisson observation term

We explain how to obtain the dynamics of the posterior density ZtZ_{t} averaged over the stochastic realization of true gradient θt\theta_{t} dynamics and approximate the jump process dYt\mathrm{d}Y_{t} with a diffusion process. Because the derivation can be applied not only to our setting but also to other filtering process, we explain the result based on a more general setting.

We formulate a filtering problem by modeling hidden and observation processes. We consider a hidden process whose probability density function t(x)=(Xt=x)\mathbb{P}_{t}(x)=\mathbb{P}(X_{t}=x) follows the Fokker-Planck equation

t(x)t=Xt(x)\displaystyle\frac{\partial\mathbb{P}_{t}(x)}{\partial t}=\mathcal{L}_{X}^{\ast}\mathbb{P}_{t}(x) (46)

where \mathcal{L}^{\ast} is the Fokker-Planck operator. We model observation processes YtiY_{t}^{i} with i{1,,nY}i\in\{1,\ldots,n_{Y}\} by the Poisson processes whose rate processes are λti=λi(Xt,Zt)\lambda_{t}^{i}=\lambda^{i}(X_{t},Z_{t^{-}}), i.e., (dYtiXt,Y0:t)=λi(Xt,Zt)dt\mathbb{P}(\mathrm{d}Y_{t}^{i}\mid X_{t},Y_{0:t^{-}})=\lambda^{i}(X_{t},Z_{t^{-}})\mathrm{d}t. Note that Zt(x)=(Xt=xY0:t)Z_{t^{-}}(x)=\mathbb{P}(X_{t^{-}}=x\mid Y_{0:t^{-}}) is a functional of Y0:t.Y_{0:t^{-}}. The posterior density function Zt(x)Z_{t}(x) under this setting follows the filtering equation (see [21] for pedagogical introduction of the filtering equation):

dZt(x)=XZt(x)dt+Zt(x)j=1NYλi(x,Zt)λ^i[Zt]λ^i[Zt](dYtiλ^i[Zt]dt)\displaystyle\mathrm{d}Z_{t}(x)=\mathcal{L}_{X}^{\ast}Z_{t}(x)\mathrm{d}t+Z_{t^{-}}(x)\sum_{j=1}^{N_{Y}}\frac{\lambda^{i}(x,Z_{t^{-}})-\hat{\lambda}^{i}[Z_{t^{-}}]}{\hat{\lambda}^{i}[Z_{t^{-}}]}(\mathrm{d}Y_{t}^{i}-\hat{\lambda}^{i}[Z_{t^{-}}]\mathrm{d}t) (47)

where f^[z]:=𝒟xf(x,z)z(x)\hat{f}[z]:=\int\mathcal{D}xf(x,z)z(x) represents the posterior expectation of ff when the posterior density is Zt=zZ_{t}=z. The problem in the main text corresponds to the case where the hidden process is Xt=θtX_{t}=\theta_{t} with XZt(X)=2ZtX2(X)\mathcal{L}^{\ast}_{X}Z_{t}(X)=\frac{\partial^{2}Z_{t}}{\partial X^{2}}(X) and the observation process is Yti=YtϕiY_{t}^{i}=Y_{t}^{\phi_{i}} with NY=NN_{Y}=N and λi(θ,z)=(λ¯/N)(1+αrϕi[z]𝐞(ϕi),𝐞(θ)\lambda^{i}(\theta,z)=(\bar{\lambda}/N)\cdot(1+\alpha r^{\phi_{i}}[z]\left\langle\bm{\mathrm{e}}(\phi_{i}),\bm{\mathrm{e}}(\theta)\right\rangle where rϕi[z]=1+σϵcos(2(ϕθ^[z]))r^{\phi_{i}}[z]=1+\sigma\epsilon\cos(2(\phi-\hat{\theta}[z])) and θ^[z]\hat{\theta}[z] is the posterior expectation of θt\theta_{t} when the posterior density is Zt=zZ_{t}=z.

We then consider the time evolution of the joint density function of XtX_{t} and ZtZ_{t}: t(x,z):=(Xt=x,Zt=z)\mathbb{P}_{t}(x,z):=\mathbb{P}(X_{t}=x,Z_{t}=z). Assuming that Zt()Z_{t}(\cdot) can be treated in the same way as a finite-dimensional stochastic variable, the dynamics of the joint density is formally represented by the following forward Kolmogorov equation:

t(x,z)t\displaystyle\frac{\partial\mathbb{P}_{t}(x,z)}{\partial t} =Xt(x,z)𝒟x~δδz(x~)[Xz(x~)z(x~)j=1NY(λi(x~,z)λ^i[z])]t(x,z)\displaystyle=\mathcal{L}_{X}^{\ast}\mathbb{P}_{t}(x,z)-\int\mathcal{D}\tilde{x}\frac{\delta}{\delta z(\tilde{x})}\left[\mathcal{L}_{X}^{\ast}z(\tilde{x})-z(\tilde{x})\sum_{j=1}^{N_{Y}}(\lambda^{i}(\tilde{x},z)-\hat{\lambda}^{i}[z])\right]\mathbb{P}_{t}(x,z) (48)
+𝒟wi{λZi(wi,x,zwi)t(x,zwi)λZi(wi,x,z)t(x,z)}\displaystyle+\int\mathcal{D}w^{i}\left\{\lambda_{Z}^{i}(w^{i},x,z-w^{i})\mathbb{P}_{t}(x,z-w^{i})-\lambda_{Z}^{i}(w^{i},x,z)\mathbb{P}_{t}(x,z)\right\} (49)

where λZi(w,x,z):=λi(x,z)δ(wi()z()λi(,z)λ^i[z]λ^i[z])\lambda_{Z}^{i}(w,x,z):=\lambda^{i}(x,z)\delta\left(w^{i}(\cdot)-z(\cdot)\frac{\lambda^{i}(\cdot,z)-\hat{\lambda}^{i}[z]}{\hat{\lambda}^{i}[z]}\right) describes the rate of jump of zz with size ww when Xt=xX_{t}=x and Zt=zZ_{t}=z, and δ/δz(x~)\delta/\delta z(\tilde{x}) represents a functional derivative, which describes how a functional changes with respect to variations in its input function zz at the point x~\tilde{x} (see Ref. [28], Sec. II.C.5, for the treatment of functional derivatives in the context of stochastic processes). We note that there is a distinction between xx and x~\tilde{x} which represent the realization of the hidden process XtX_{t} and the argument of the posterior density ZtZ_{t}, respectively.

To marginalize over hidden variable XtX_{t}, we integrate each term with respect to xx by noting that t(x,z)=t(xz)t(z)=z(x)t(z)\mathbb{P}_{t}(x,z)=\mathbb{P}_{t}(x\mid z)\mathbb{P}_{t}(z)=z(x)\mathbb{P}_{t}(z):

t(z)t\displaystyle\frac{\partial\mathbb{P}_{t}(z)}{\partial t} =𝒟x~δδz(x~)[Xz(x~)t(z)z(x~)j=1NY(λi(x~,z)λ^i[z])t(z)]\displaystyle=-\int\mathcal{D}\tilde{x}\frac{\delta}{\delta z(\tilde{x})}\left[\mathcal{L}_{X}^{\ast}z(\tilde{x})\mathbb{P}_{t}(z)-z(\tilde{x})\sum_{j=1}^{N_{Y}}(\lambda^{i}(\tilde{x},z)-\hat{\lambda}^{i}[z])\mathbb{P}_{t}(z)\right] (50)
+j=1NY{𝒟wi{λ^Zi[wi,zwi]t(zwi)λ^Zi[wi,z]t(z)}},\displaystyle+\sum_{j=1}^{N_{Y}}\left\{\int\mathcal{D}w^{i}\left\{\hat{\lambda}_{Z}^{i}[w^{i},z-w^{i}]\mathbb{P}_{t}(z-w^{i})-\hat{\lambda}_{Z}^{i}[w^{i},z]\mathbb{P}_{t}(z)\right\}\right\}, (51)

where λ^Zi[w,z]:=λ^i[z]δ(wi()z()λi(,z)λ^i[z]λ^i[z])\hat{\lambda}_{Z}^{i}[w,z]:=\hat{\lambda}^{i}[z]\delta\left(w^{i}(\cdot)-z(\cdot)\frac{\lambda^{i}(\cdot,z)-\hat{\lambda}^{i}[z]}{\hat{\lambda}^{i}[z]}\right) represents the posterior expectation of λZi(w,x,z)\lambda_{Z}^{i}(w,x,z).

Next, we approximate the jump term by a diffusion based on the Kramers-Moyal expansion. To do that, we add an assumption on the observation process. The jump rate λi\lambda^{i} in the observation process can be decomposed into a sum of XtX_{t}-dependent and XtX_{t}-independent terms: λi(x,z)=λ¯i(1+αi(x,z))\lambda^{i}(x,z)=\bar{\lambda}^{i}(1+\alpha^{i}(x,z)). For the XtX_{t} dependent term, we assume that αi(x,z)\alpha^{i}(x,z) is sufficiently small so that αi(x,z)=O(α)\alpha^{i}(x,z)=O(\alpha) holds for any jj and xx with a small parameter α\alpha and that terms of order o(α2)o(\alpha^{2}) are ignorable. Then, we expand the jump term in the Taylor series about the jump size wiw^{i}:

𝒟wi{λ^Zi[wi,zwi]t(zwi)λ^Zi[wi,z]t(z)}\displaystyle\int\mathcal{D}w^{i}\left\{\hat{\lambda}_{Z}^{i}[w^{i},z-w^{i}]\mathbb{P}_{t}(z-w^{i})-\hat{\lambda}_{Z}^{i}[w^{i},z]\mathbb{P}_{t}(z)\right\} =𝒟wi𝒟x~p|p|=1(wi(x~))p|p|!δpδz(x~)p[λ^Zi[wi,z]t(z)]\displaystyle=\int\mathcal{D}w^{i}\int\mathcal{D}\tilde{x}^{p}\sum_{|p|=1}^{\infty}\frac{(-w^{i}(\tilde{x}))^{p}}{|p|!}\frac{\delta^{p}}{\delta z(\tilde{x})^{p}}\left[\hat{\lambda}_{Z}^{i}[w^{i},z]\mathbb{P}_{t}(z)\right] (52)
=𝒟x~p|p|=1(1)|p||p|!δpδz(x~)p[𝒟wiwi(x~)pλ^Zi[wi,z]t(z)]\displaystyle=\int\mathcal{D}\tilde{x}^{p}\sum_{|p|=1}^{\infty}\frac{(-1)^{|p|}}{|p|!}\frac{\delta^{p}}{\delta z(\tilde{x})^{p}}\left[\int\mathcal{D}w^{i}w^{i}(\tilde{x})^{p}\hat{\lambda}_{Z}^{i}[w^{i},z]\mathbb{P}_{t}(z)\right] (53)
=𝒟x~p|p|=1(1)|p||p|!δpδz(x~)p[λ^i[z](z(x~)λi(x~,z)λ^i[z]λ^i[z])pt(z)]\displaystyle=\int\mathcal{D}\tilde{x}^{p}\sum_{|p|=1}^{\infty}\frac{(-1)^{|p|}}{|p|!}\frac{\delta^{p}}{\delta z(\tilde{x})^{p}}\left[\hat{\lambda}^{i}[z]\left(z(\tilde{x})\frac{\lambda^{i}(\tilde{x},z)-\hat{\lambda}^{i}[z]}{\hat{\lambda}^{i}[z]}\right)^{p}\mathbb{P}_{t}(z)\right] (54)

where pp is a multi-index. By noting that (z(x~)λi(x~,z)λ^i[z]λ^i[z])p\left(z(\tilde{x})\frac{\lambda^{i}(\tilde{x},z)-\hat{\lambda}^{i}[z]}{\hat{\lambda}^{i}[z]}\right)^{p} is of the order O(α|p|)O(\alpha^{|p|}), we ignore the third and higher order terms with respect to |p||p|. Because the first order term of Eq. (54) cancels with the term j=1NY(λi(x~,z)λ^i[z])t(z)-\sum_{j=1}^{N_{Y}}(\lambda^{i}(\tilde{x},z)-\hat{\lambda}^{i}[z])\mathbb{P}_{t}(z) in the first line of Eq. (51), we obtain the dynamics of t(z)\mathbb{P}_{t}(z) as the following Fokker-Planck equation:

t(z)t\displaystyle\frac{\partial\mathbb{P}_{t}(z)}{\partial t} =𝒟x~δδz(x~)[Xz(x~)t(z)]\displaystyle=-\int\mathcal{D}\tilde{x}\frac{\delta}{\delta z(\tilde{x})}\left[\mathcal{L}_{X}^{\ast}z(\tilde{x})\mathbb{P}_{t}(z)\right] (55)
+12j=1NY𝒟x~𝒟x~δ2δz(x~)δz(x~)[z(x~)(λi(x~,z)λ^i[z])z(x~)(λi(x~,z)λ^i[z])λ^i[z]t(z)]\displaystyle\quad+\frac{1}{2}\sum_{j=1}^{N_{Y}}\int\mathcal{D}\tilde{x}\mathcal{D}\tilde{x}^{\prime}\frac{\delta^{2}}{\delta z(\tilde{x})\delta z(\tilde{x}^{\prime})}\left[\frac{z(\tilde{x})(\lambda^{i}(\tilde{x},z)-\hat{\lambda}^{i}[z])\cdot z(\tilde{x}^{\prime})(\lambda^{i}(\tilde{x}^{\prime},z)-\hat{\lambda}^{i}[z])}{\hat{\lambda}^{i}[z]}\mathbb{P}_{t}(z)\right] (56)
=𝒟x~δδz(x~)[Xz(x~)t(z)]\displaystyle=-\int\mathcal{D}\tilde{x}\frac{\delta}{\delta z(\tilde{x})}\left[\mathcal{L}_{X}^{\ast}z(\tilde{x})\mathbb{P}_{t}(z)\right] (57)
+12j=1NY𝒟x~𝒟x~δ2δz(x~)δz(x~)λ¯i[z(x~)(αi(x~,z)α^i[z])z(x~)(αi(x~,z)α^i[z])t(z)]\displaystyle\quad+\frac{1}{2}\sum_{j=1}^{N_{Y}}\int\mathcal{D}\tilde{x}\mathcal{D}\tilde{x}^{\prime}\frac{\delta^{2}}{\delta z(\tilde{x})\delta z(\tilde{x}^{\prime})}\bar{\lambda}^{i}\left[z(\tilde{x})(\alpha^{i}(\tilde{x},z)-\hat{\alpha}^{i}[z])\cdot z(\tilde{x}^{\prime})(\alpha^{i}(\tilde{x},z^{\prime})-\hat{\alpha}^{i}[z])\mathbb{P}_{t}(z)\right] (58)

This Fokker-Planck equation describes the following SDE of ZtZ_{t}:

dZt(x)\displaystyle\mathrm{d}Z_{t}(x) =XZt(x)dt+Zt(x)j=1NYλ¯{αi(x,z)α^i[Zt]}dWti\displaystyle=\mathcal{L}_{X}^{\ast}Z_{t}(x)\mathrm{d}t+Z_{t}(x)\sum_{j=1}^{N_{Y}}\sqrt{\bar{\lambda}}\left\{\alpha^{i}(x,z)-\hat{\alpha}^{i}[Z_{t}]\right\}\mathrm{d}W_{t}^{i} (59)

where dWti\mathrm{d}W_{t}^{i} are the standard Wiener processes. By applying the result to our problem in the main text, the prediction term does not change from Eq. (13) and the observation process is transformed from Eq. (14) to Eq. (59).

Appendix C Relaxation of the estimation error to the stationary state

To show that the estimation error reaches a stationary state after time depending on the parameter Λ\Lambda, we plot the estimation error as a function of the scaled time Λ1/2t\Lambda^{1/2}t in Figure C.4.

Refer to caption
Figure C.4: Estimation error as a function of the scaled time Λ1/2t\Lambda^{1/2}t for different values of the signal-to-noise ratio Λ\Lambda and the shape parameter σϵ\sigma\epsilon. (A1, A2) Estimation error for an elliptic shape (A1) and a shape described by the second Fourier component (A2). Different curves represent different values of σϵ\sigma\epsilon and Λ\Lambda, obtained by varying the ligand binding rate λ¯\bar{\lambda} while keeping the gradient steepness α\alpha fixed. (B1, B2) Estimation error for an elliptic shape (B1) and a shape described by the second Fourier component (B2). Different curves represent different values of ϵ\epsilon and Λ\Lambda, obtained by varying the gradient steepness α\alpha while keeping the ligand binding rate λ¯\bar{\lambda} fixed. In all panels, solid, dotted, and dashed curves represent σϵ=0\sigma\epsilon=0, σϵ=0.7\sigma\epsilon=0.7, and σϵ=0.7\sigma\epsilon=-0.7, respectively. The shaded regions around curves indicate ±\pm the standard error estimated by bootstrap sampling. The vertical dashed lines indicate the time interval used for time averaging, from t=2Λ1/2t=2\Lambda^{-1/2} to t=3Λ1/2t=3\Lambda^{-1/2}.

Appendix D Sample trajectories of the filtering processes under low signal-to-noise ratio regime

To complement the parameter dependence of the estimation error, Figure D.5 shows the time trajectory of the true and estimated gradient direction with the signal-to-noise ratio being lower than the value used in Figure 1.

Refer to caption
Figure D.5: Gradient sensing processes under a signal-to-noise ratio lower than that in Fig. 1. (A1,A2) Time lapse of a sample process of the gradient sensing when the signal-to-noise ratio is decreased to Λ=1\Lambda=1 by reducing λ¯\bar{\lambda} from 21032\cdot 10^{3} to 21022\cdot 10^{2} (A1) and by reducing α\alpha from 10110^{-1} to 101.510^{-1.5} (A2). Blue and red arrows indicate the true and estimated gradient direction, respectively. The length of the red arrows indicate the confidence of estimate |𝐞^t||\hat{\bm{\mathrm{e}}}_{t}|. Orange dots show the locations of receptors at which ligand bound during time [t,t+Δt)[t,t+\Delta t) with Δt\Delta t being the discretization width. (B1,B2)Time evolution of variables in the sample process shown in (A1,A2) when the signal-to-noise ratio is decreased to Λ=1\Lambda=1 by reducing λ¯\bar{\lambda} from 21032\cdot 10^{3} to 21022\cdot 10^{2} (A1) and by reducing α\alpha from 10110^{-1} to 101.510^{-1.5} (A2). The first row shows the true gradient direction θt\theta_{t} (blue curve) and the angular locations of receptors at which ligand binds dYt\mathrm{d}Y_{t} (orange dots). The second row shows the posterior density of the gradient direction ZtZ_{t} and the estimated gradient direction θ^t\hat{\theta}_{t} (red line). The third row shows the true (blue curve) and estimated (red curve) gradient direction. The final row shows the error between the true and estimated direction measured by (θtθ^t)2/2(\theta_{t}-\hat{\theta}_{t})^{2}/2 (black curve) and the uncertainty of estimate represented by 1|𝐞^t|1-|\hat{\bm{\mathrm{e}}}_{t}| (red curve). The other parameter is set to σϵ=0.4\sigma\epsilon=-0.4 during the simulation.