Statistical Tests for Force Inference in Heterogeneous Environments

TL;DR

This paper introduces a Bayesian inference method to detect and estimate forces in heterogeneous environments from stochastic trajectories, accounting for spurious forces due to spatially varying diffusivity, and validates it analytically, numerically, and experimentally.

Contribution

It presents a novel Bayesian approach for force inference that marginalizes over spurious forces and incorporates a statistical test, implemented in open-source software.

Findings

Method reliably detects forces in heterogeneous environments.

Analytical and numerical validation shows high accuracy.

Tested successfully on experimental biomolecule trajectory data.

Abstract

We devise a method to detect and estimate forces in a heterogeneous environment based on experimentally recorded stochastic trajectories. In particular, we focus on systems modeled by the heterogeneous overdamped Langevin equation. Here, the observed drift includes a "spurious" force term when the diffusivity varies in space. We show how Bayesian inference can be leveraged to reliably infer forces by taking into account such spurious forces of unknown amplitude as well as experimental sources of error. The method is based on marginalizing the force posterior over all possible spurious force contributions. The approach is combined with a Bayes factor statistical test for the presence of forces. The performance of our method is investigated analytically, numerically and tested on experimental data sets. The main results are obtained in a closed form allowing for direct exploration of their properties and fast computation. The method is incorporated into TRamWAy, an open-source software platform for automated analysis of biomolecule trajectories.

Figures

Equations

$$d{\mathbf{X}}_t=\frac{\mathbf{f}\left({\mathbf{X}}_t\right)}{\gamma \left({\mathbf{X}}_t\right)}dt+\sqrt{2D\left({\mathbf{X}}_t\right)}\circ d{\mathbf{W}}_t,$$

$$d{\mathbf{X}}_t=\mathbf{a}({\mathbf{X}}_t)dt+b({\mathbf{X}}_t)\circ d{\mathbf{W}}_t,$$

$$\begin{gathered} {\mathbf{X}}_T={\mathbf{X}}_0+\underset{N\to \infty }{\lim }[\sum _{i=0}^{N-1}\mathbf{a}({\mathbf{X}}_{{\xi }_i})(t_{i+1}-t_i) \\ +\sum _{i=0}^{N-1}b({\mathbf{X}}_{{\xi }_i})({\mathbf{W}}_{t_{i+1}}-{\mathbf{W}}_{t_i})], \end{gathered}$$

$$d{\mathbf{X}}_t=\mathbf{\alpha }({\mathbf{X}}_t,\lambda )dt+b({\mathbf{X}}_t)d{\mathbf{W}}_t,$$

$$\mathbf{\alpha }(\mathbf{X},\lambda )\equiv \mathbf{a}(\mathbf{X})+\lambda b(\mathbf{X})\nabla b(\mathbf{X}).$$

$$\mathrm{Pr}(H_i\mid T)=\frac{p(T\mid H_i)\pi (H_i)}{p(T)}.$$

$$\frac{\mathrm{Pr}(H_1\mid T)}{\mathrm{Pr}(H_0\mid T)}=\frac{p(T\mid H_1)}{p(T\mid H_0)}\frac{\pi (H_1)}{\pi (H_0)}.$$

$$K\equiv \frac{p(T\mid H_1)}{p(T\mid H_0)}.$$

$$p(T\mid H_i)=\int d{\mathbf{\theta }}_{i}p(T\mid {\mathbf{\theta }}_i,H_i)\pi ({\mathbf{\theta }}_i\mid H_i).$$

$$\begin{gathered} p(\{\Delta \mathbf{r}\}\mid \mathbf{\alpha },b^2) \\ =(2\pi b^2\Delta t{)}^{-nd/2}\exp \left(-\frac{n(\overline{\Delta \mathbf{r}}-\mathbf{\alpha }\Delta t{)}^2+nV}{2b^2\Delta t}\right). \end{gathered}$$

$$\pi (\mathbf{a},\lambda ,\mathbf{g},b^2\mid H_1)\approx \pi (a,b^2\mid \lambda ,H_1)\pi (\lambda )\pi (\mathbf{g}).$$

$$\begin{gathered} p(\{\Delta \mathbf{r}\}\mid H_0) \\ =\frac{A_d{C}_d}{n^{d/2}}{\int }_0^{1}d\lambda {\left(nV+n_{\pi }{V}_{\pi }+n{\left(\overline{\Delta \mathbf{r}}-\lambda \mathbf{g}\Delta t\right)}^2\right)}^{-\kappa (d)}. \end{gathered}$$

$$K^{\mathrm{M}}$$

$$K(\lambda )$$

$$\begin{gathered} p({\zeta }_a\mid T) \\ \propto {\int }_0^{1}d\lambda {\left[v+\frac{{\eta }^2}{1-{\eta }^2}{\zeta }_a^2+({\zeta }_{\mathrm{t}}-{\zeta }_a-\lambda {\zeta }_{\mathrm{sp}}{)}^2\right]}^{-\kappa (d)-d/2}, \end{gathered}$$

$$p(b^2\mid T)=\frac{{\left(\frac{nV{G}^{\prime \prime \prime }}{2b^2\Delta t}\right)}^{\kappa (d)}\exp \left(-\frac{nV{G}^{\prime \prime \prime }}{2b^2\Delta t}\right)}{b^2\gamma \left(\kappa (d),\frac{nV{G}^{\prime \prime \prime }}{4{\sigma }_L^2}\right)}.$$

$$K^{\mathrm{M}}\equiv \frac{{\int }_{\eta (R-1)}^{\eta R}d{q}^{\prime }(v+{\zeta }_{sp}^2{q}^{\prime 2}{)}^{-\kappa (d)}}{{\int }_{R-1}^{R}dq(v+{\zeta }_{sp}^2{q}^2{)}^{-\kappa (d)}}.$$

Full text

Keywords: overdamped Langevin equation, Bayesian inference, inverse problems, biomolecule dynamics, Itô-Stratonovich dilemma, random walks

1 Introduction

Random walks are encountered throughout biology and other domains of science, and so is the associated inverse problem of inferring their properties from experimental data. Random walkers can be considered probes of their environment, and their recorded trajectories thus contain information on the properties of both the walker and its environment. In the context of biophysics, the random walkers are typically colloidal particles or biomolecules, but, in a general context, they may, for example, represent the motion along an abstract coordinate of a chemical reaction or the fluctuating price of a stock asset. Transport of biomolecules within cells [79], conformational dynamics of proteins and RNA molecules [6], diffusion of proteins on the DNA [76], dynamics of nanosized objects in the cytosol [20], dynamics of receptors in neurons  [17, 13, 66], complex random walks in mixed biological environments [71, 30, 35, 58], bacteria performing chemotaxis [53, 80], immune-cell dynamics [65, 64, 26], and directionally persistent cell movement [50] are all examples of cases where biologically relevant information can be extracted from recorded stochastic trajectories.

Empirical systems featuring biomolecule random walks are typically characterized by high heterogeneity, so the inverse problem often translates into inferring the properties of the heterogeneous environment from the trajectories of tracer molecules [14, 5, 19, 52, 34, 11, 56, 45, 48, 70, 61, 29, 12, 55]. A paradigmatic model for random walks in such systems is the heterogeneous overdamped Langevin equation (OLE):

[TABLE]

which describes the continuous-time dynamics underlying a discrete-time recorded random walk [34]. Here ${\mathbf{X}}_{t}$ is the tracer’s position at time $t$, ${\mathbf{f}}({\mathbf{X}})$ is the force acting on it in the point ${\mathbf{X}}$, $D({\mathbf{X}})$ is its diffusivity, $\gamma({\mathbf{X}})$ is the viscous friction coefficient, and ${\mathbf{W}}_{t}$ is a Gaussian zero-mean continuous-time white noise process with uncorrelated increments and unit variance [42]. Owing to its simplicity, the OLE (1) is a popular model for biological random walks, providing an effective mesoscopic description of the dynamics [34]. As is often the case for models of biological systems, the OLE is empirically postulated rather than derived from the first principles, since this derivation is complex and requires taking into account many factors, such as the heterogeneity of the environment composition, presence of boundaries [44] and hydrodynamic properties [16], as well as possible noise correlations [24, 4]. We refer the interested reader to (i) [33, 63] for an in-depth discussion of the derivation of the OLE and for a microscopic model of crowded environments, to (ii) [74, 36, 46] for some approaches to the derivation of the equation of motion in media featuring diffusivity or temperature gradients, and to (iii) [43, 4, 24] for a discussion and experimental measurements of the extent to which Brownian noise is truly uncorrelated.

When the diffusivity $D({\mathbf{x}})$ varies in space, Eq. (1) is only well-defined after a convention for calculating the (stochastic) integral of the noise term has been defined [73]. Two well-known examples are the Itô and Stratonovich conventions. Each convention leads to a different extra “spurious” drift term, which is proportional to the diffusivity gradient [73, 46, 81]. This feature of the OLE is known as the Itô-Stratonovich dilemma [75, 46, 69, 21, 22, 60]. For an experimental illustration of the presence of two physically different components see [78, 8].

Any stochastic convention can be used in the OLE to statistically describe a given experimental random walk. However, they are not equivalent to an external observer attempting to interpret the parameters of the random walk. Indeed, the two components of the drift in the OLE — the spurious force and the non-diffusive force — have a different physical nature. The spurious force is proportional to the diffusivity gradient, hence its value changes when the diffusivity, viscosity, or the temperature of the system change in space (time) or across systems. On the contrary, the non-diffusive component of the drift does not depend on the diffusivity and represents specific or non-specific interactions. The spurious force does not contribute to the equilibrium (Boltzmann) distribution.

Spurious forces are due to the interactions of the tracked particle with the surrounding thermal bath, while the non-diffusive forces represent its interactions with other objects or fields [75, see also the discussion of the “internal noise” in]. The separation of interactions between these two groups naturally depends (i) on the scale on which the system is analyzed and (ii) on which parts of the environment are included into the thermal bath. In Sect. 3.7 below, we show how in the same simulated system, the drift due to non-diffusive forces on the microscopic scale is perceived as a spurious force on the mesoscopic scale, when the contribution of individual interacting partners can no longer be identified.

For practical applications, it is thus important to develop a method allowing to distinguish between diffusive and non-diffusive forces or at least to develop a test allowing to confirm the presence of the non-diffusive forces on a given scale. The need for such approaches is further emphasized by the inaccessibility of the equilibrium distributions and of the exact boundary conditions at the nanometer scale in numerous biological setups.

Since the seminal work of [1], the inverse problem of drift and diffusivity inference from random walks has been attracting attention [27, 38], especially in financial applications [15, 59]. In this article, we address a more specific problem of distinguishing between diffusive and non-diffusive forces, since the value of the latter is generally unknown. Although the spurious force is proportional to the diffusivity gradient $\lambda\nabla D({\mathbf{X}})$, it includes an unknown proportionality factor $\lambda$. It is known that for physical systems in equilibrium, described by Boltzmann distribution, $\lambda=1$, but its value is not known in general for out-of-equilibrium systems [31, 39, 43, 46, 78, 81, 21, 22, 60]. Each value of $\lambda$ represents specific symmetries of transition probabilities in these systems [39, 40, 69].

Our goal is hence two-fold: (i) to develop a statistical test for the presence of non-diffusive forces, and (ii) to infer the posterior distribution for the intensity of the non-diffusive forces while taking into account all possible contributions of the “spurious” forces as well as experimental localization errors and motion blur. The method we introduce here is statistically robust to changes in the spurious force contribution in the OLE due, for example, to changes of the diffusivity or viscosity. We validate our approach on numerical trajectories and demonstrate its efficiency on experimental data.

2 The Itô-Stratonovich Dilemma for the Inverse Problem

In this section we give a brief review of the Itô-Stratonovich dilemma [75]. Numerous discussions of the dilemma have been focused on choosing the appropriate integral convention for the forward problem of integrating the OLE in a particular system. In contrast, we here focus on how the dilemma affects the inverse problem of inferring the underlying physical parameters of a model from recorded data. To underline the generality of the problem, we rewrite the OLE (1) in the form of a general stochastic differential equation (SDE):

[TABLE]

where ${\mathbf{a}}$ and $b$ are differentiable functions of ${\mathbf{X}}_{t}$. We will refer to ${\mathbf{a}}$ and $b$ as the drift and diffusivity respectively.

The integral of Eq. (2) is defined as the limit of Riemann sums

[TABLE]

where each point $\xi_{i}$ is chosen in the interval $[t_{i};t_{i+1}]$. The standard conventions — Itô, Stratonovich-Fisk and Hänggi-Klimontovich — correspond to $\xi_{i}=t_{i}$, $\xi_{i}=(t_{i}+t_{i+1})/2$ and $\xi_{i}=t_{i+1}$ respectively [69, 46]. More generally, $\xi_{i}$ can be set to any point $\xi_{i}=t_{i}+\lambda(t_{i+1}-t_{i})$ within the $[t_{i};t_{i+1}]$ interval. This allows one to rewrite Eq. (2) with any convention $\lambda$ in the Itô form [77, 41]:

[TABLE]

where the total drift $\alpha$ is the sum of ${\mathbf{a}}$ and the spurious drift $\lambda b({\mathbf{X}})\mathbf{\nabla}b({\mathbf{X}})$:

[TABLE]

From the perspective of the forward problem, Eq. (5) shows that the often arbitrary choice of the value of $\lambda$ influences the value of the drift $\bm{\alpha}$ when ${\mathbf{a}}$ and $b$ are fixed, — this is the essence of the Itô-Stratonovich dilemma [75]. In the context of the inverse problem, one is given fixed values of $\bm{\alpha}$ and $b$ estimated from the recorded trajectories, so different choices of $\lambda$ result in different estimates of the non-diffusive drift ${\mathbf{a}}$. If the chosen $\lambda$ does not agree with its true value in the empirical system, the resulting estimate of ${\mathbf{a}}$ becomes biased.

We emphasize that we do not address here the forward problem, i.e. the question of finding the correct $\lambda$ for a given system [75, 46, 69, 21, 22, 60, discussed at length elsewhere, see]. The correct $\lambda$ values are often inaccessible in real biological systems. Instead, we aim to solve the inverse problem of whether non-diffusive forces are observed in the system and to infer their values if the appropriate value of $\lambda$ cannot be determined. It is an inverse problem with an uncertainty in the underlying physical model. This ambiguity in $\lambda$ may stem, for example, from the lack of a priori knowledge about the out-of-equilibrium fluxes in the system, noise correlations or the particle density distribution. In all cases, the method developed below allows one to obtain estimates of the non-diffusive forces and to circumvent the Itô-Stratonovich dilemma by marginalizing over all possible $\lambda$ values. The estimates are robust to changes in the spurious force contribution in the OLE.

Above, we have formulated the main question of this paper from a physical point of view as that of inferring non-diffusive forces, when the correct $\lambda$ is unknown. It is interesting to note that the same question can also be asked from a purely statistical point of view: Given the OLE, does there exist a value of $0\leqslant\lambda\leqslant 1$ that would allow to describe the given system with zero non-diffusive forces $(\mathbf{a}=0)$? This would allow to describe the same system with fewer parameters ($D$ and $\lambda$ instead of $D$, $\lambda$, $\mathbf{a}$), thus minimizing the description length among all the descriptions proposed by the OLE family [2, 62].

From this point of view, the Bayes factor developed below is a Bayesian analog of the difference in the description lengths between the models with $\mathbf{a}\neq 0$ and $\mathbf{a}=0$ for the given data. It evaluates how much more efficient the non-diffusive-force description is, as compared to the spurious-force-only description of the same data. If the spurious-force description is preferred, as a byproduct, one can calculate the value of $\lambda$ that provides the most efficient description of the data.

3 The Bayesian Approach

Our goal is to discriminate between the following two nested hypotheses:

• •

$H_{0}$: the only forces present are spurious forces due to heterogeneous diffusivity (the null hypothesis).

• •

$H_{1}$: there are other, non-diffusive, forces acting on the random walker in addition to the spurious forces.

We use the Bayes factor to decide between these hypotheses [84].

3.1 The Bayes factor

According to Bayes’ rule [28], the posterior probability $\operatorname{\mathrm{Pr}}(H_{i}\mid T)$ of a hypothesis $H_{i}$ given data $T$ is

[TABLE]

Here, $T$ is a trajectory, $T\equiv\{{{\mathbf{r}}}_{i}\}_{i=1}^{n}$, or a set of trajectories; $\operatorname{\mathit{p}}(T\mid H_{i})$ is the marginal likelihood for the data $T$ to be observed under the hypothesis $H_{i}$; $\pi(H_{i})$ is the prior probability of $H_{i}$; and $\operatorname{\mathit{p}}(T)$ is the probability to observe $T$ under either hypothesis. For the two competing hypotheses $H_{1}$ and $H_{0}$, the ratio of their posterior probabilities reads

[TABLE]

The first fraction on the right-hand side is called the Bayes factor for $H_{1}$ over $H_{0}$ [84]:

[TABLE]

Each marginal likelihood $\operatorname{\mathit{p}}(T\mid H_{i})$ is calculated by marginalizing the corresponding conditional likelihood $\operatorname{\mathit{p}}(T\mid{\bm{\theta}}_{i},H_{i})$ over all model parameters:

[TABLE]

For $H_{0}$, the likelihood $\operatorname{\mathit{p}}(T\mid{\bm{\theta}}_{0},H_{0})$ thus depends on 3 parameters: ${\bm{\theta}}_{0}=\{b^{2},{\mathbf{g}},\lambda\}$, where $b^{2}$ is the diffusivity and ${\mathbf{g}}\equiv\nabla b$ is the diffusivity gradient. For $H_{1}$, the likelihood $\operatorname{\mathit{p}}(T\mid{\bm{\theta}}_{1},H_{1})$ additionally includes the drift ${\mathbf{a}}$, so ${\bm{\theta}}_{1}=\{b^{2},{\mathbf{g}},\lambda,{\mathbf{a}}\}$. Note that we treat ${\mathbf{g}}$ as independent from $b^{2}$, which allows us to obtain the results in the analytical form. This assumption is further discussed in Appendix A1.

3.2 Likelihood

The likelihood $\operatorname{\mathit{p}}(T\mid{\bm{\theta}}_{i},H_{i})$ is obtained as the fundamental solution of the Fokker-Planck equation corresponding to the OLE (2). However, it cannot in general be obtained analytically. Instead, one can approximate it locally by assuming that $\bm{\alpha}$ and $b$ are constant within small spatial domains [52, 19]. In this case, the likelihood of observing a set of displacements $\{\Delta\mathbf{r}\}$ inside a given domain is [52]:

[TABLE]

Here the mean displacement $\overline{{\Delta{\mathbf{r}}}}\equiv\sum_{i=1}^{n}{\Delta{\mathbf{r}}}_{i}/n$ and the biased sample variance $V\equiv\sum_{i=1}^{n}|{\Delta{\mathbf{r}}}_{i}-\overline{{\Delta{\mathbf{r}}}}|^{2}/n$ are the sufficient statistics of the model [28], and $d$ is the number of dimensions. The equations below are valid for $d=1$ and $d=2$, but the framework can also be extended to $d=3$.

Note that calculations would be similar if one relaxed the approximation of the locally constant values of $\bm{\alpha}$ and $b$. Computations would be performed numerically but the analytical explorations such as those of Appendix A3 would not be possible. Meanwhile, the assumption of bin independence is paramount to the presented method.

3.3 Priors

The likelihood (6) belongs to the exponential family [28]. Therefore, a natural choice for the prior is a conjugate prior for the parameters ${\mathbf{a}}$ and $b^{2}$. Among other advantages, conjugate priors provide a closed form of the posterior distribution. We furthermore assume a factorized form for the prior distributions for $\lambda$ and the diffusivity gradient ${\mathbf{g}}$:

[TABLE]

We have no a priori information available about the true value of $\lambda$ other than that $0\leqslant\lambda\leqslant 1$, so we use the flat prior $\pi(\lambda)\equiv 1$. The diffusivity gradient prior is approximated by a delta function $\pi({\mathbf{g}})\equiv\delta({\mathbf{g}}-\hat{\mathbf{g}})$ centered around its maximum a posteriori (MAP) value $\hat{{\mathbf{g}}}$. Details of $\hat{g}$ estimation are given in Appendix A1.

Under $H_{1}$, the full conjugate prior is then (cf. (6)):

where $A_{d}\equiv(n_{\pi}/(2\pi))^{d/2}(n_{\pi}V_{\pi}/2)^{m(d)}\Delta t^{d+1}/\Gamma(m(d))$; $m(d)\equiv d(n_{\pi}-1)/2-1$; $\bm{\mu}_{\pi}$, $V_{\pi}$ and $n$ are the parameters of the prior (called hyper-parameters). The models $H_{0}$ and $H_{1}$ are nested models. In such case, it is common practice to obtain the $H_{0}$ prior (Eq. (LABEL:fm:H0-prior)) by integrating the $H_{1}$ prior over ${\mathbf{a}}$ [84].

We further set the hyper-parameters to maximally favor the null model. More specifically, $n_{\pi}$ acts as an effective number of prior observations. The least constraining prior is obtained by setting $n_{\pi}=4$ for 1D data and $n_{\pi}=3$ for 2D, which are the minimal number of observations, for which the prior is proper (normalized). Furthermore, we center the prior on zero force by setting $\bm{\mu}_{\pi}=\lambda{\mathbf{g}}{\Delta t}$. The remaining hyper-parameter $V_{\pi}$ defines the prior distribution for the diffusivity. Sensitivity of the results to $u\equiv V_{\pi}/V$ is explored in Appendix A2.

3.4 Model evidence and the Bayes factor

The evidence for the $H_{1}$ and $H_{0}$ models is the central ingredient in the Bayes factor. Given the likelihood ${\operatorname{\mathit{p}}(\{\Delta\mathbf{r}\}\mid\bm{\alpha}({\mathbf{a}},\lambda),b^{2})}$ (Eq. (6)) and prior ${\pi({\mathbf{a}},\lambda,{\mathbf{g}},b^{2}\mid H_{1})}$ (Eq. (LABEL:fm:H1-prior)), the evidence for $H_{1}$ is calculated by marginalizing ${\operatorname{\mathit{p}}(\{\Delta\mathbf{r}\}\mid\bm{\alpha}({\mathbf{a}},\lambda),b^{2})}{\pi({\mathbf{a}},\lambda,{\mathbf{g}},b^{2}\mid H_{1})}$ over all the parameters ${\bm{\theta}}_{1}=\{a,b^{2},g,\lambda\}$. This gives

where $C_{d}=2^{\kappa(d)}\Gamma(\kappa(d))(2\pi)^{\frac{d(1-n)}{2}}{\Delta t}^{-d-1}$ and $\kappa(d)\equiv d(n+n_{\pi}-1)/2-1$.

For $H_{0}$, the likelihood $\operatorname{\mathit{p}}(\{\Delta\mathbf{r}\}\mid\bm{\alpha}(\lambda),b^{2})$ is given by Eq. (6) with $\alpha(\lambda)=\lambda b\nabla b$, and the prior ${\pi(\lambda,{\mathbf{g}},b^{2}\mid H_{0})}$ by Eq. (LABEL:fm:H0-prior). Marginalization of ${\operatorname{\mathit{p}}(\{\Delta\mathbf{r}\}\mid\bm{\alpha}(\lambda),b^{2})}{\pi(\lambda,{\mathbf{g}},b^{2}\mid H_{0})}$ over ${\bm{\theta}}_{0}=\{b^{2},{\mathbf{g}},\lambda\}$ gives

[TABLE]

Expressions (LABEL:fm:posterior-H1) and (10) let us finally calculate the marginalized Bayes factor $K^{\mathrm{M}}$, which takes into account all possible values for the unknown parameter $\lambda$. For comparison, we also provide the Bayes factor $K(\lambda)$ for fixed-$\lambda$ inference procedures (Itô, Stratonovich or Hänggi), which is calculated in the same manner:

[TABLE]

All the integrals appearing in Eqs (LABEL:fm:posterior-H1–12) are 1D integrals that are numerically evaluated using the trapezoid rule.

The natural parameter combinations appearing in Eq. (12) are: (i) $\mathbf{\zeta}_{\mathrm{t}}\equiv\overline{{\Delta{\mathbf{r}}}}/\sqrt{V}$, the signal-to-noise ratio for the total force in a single displacement; (ii) $\mathbf{\zeta}_{\mathrm{sp}}\equiv\hat{\mathbf{g}}{\Delta t}/\sqrt{V}$, the signal-to-noise ratio for the spurious force in a single displacement; (iii) $\eta\equiv\sqrt{n_{\pi}/(n+n_{\pi})}$, the relative strength of the prior compared to the observed data; (iv) $v\equiv 1+n_{\pi}V_{\pi}/(nV)$, a weighted ratio of jump variances in the prior and in the data.

Figure 1A plots the marginalized Bayes factor $K^{\mathrm{M}}$ (12) as a function of $\zeta_{sp}$, and of the component of the total force $\zeta_{t\parallel}$ parallel to $\mathbf{\zeta}_{\mathrm{sp}}$. The lowest values of $K^{\mathrm{M}}$ are achieved in the region $0\leqslant\zeta_{t\parallel}/\zeta_{sp}\leqslant 1$. The value of $K^{\mathrm{M}}$ changes relatively little within this region but grows rapidly at its boundary. The absolute minimum of $K^{\mathrm{M}}$ is achieved for $\zeta_{sp}=0$ and $\zeta_{t\parallel}=0$ with $\min K^{\mathrm{M}}=\min K(\lambda)=\eta^{d}[(v+\eta^{2}\zeta_{t\perp}^{2})/(v+\zeta_{t\perp}^{2})]^{-\kappa}$. A mathematical analysis of Eq. (12) is provided in Appendix A3, where it is shown that non-diffusive forces cannot in principle be detected in certain intervals of $\mathbf{\zeta}_{\mathrm{sp}}$, $\mathbf{\zeta}_{\mathrm{t}}$ regardless of the number of collected data points. Appendix A4 extends the Bayes factors (12) to the experimentally relevant case with localization errors and motion blur.

3.5 Force posterior

When $H_{1}$ is met, we can infer the value of the non-diffusive force by marginalizing the force posterior over all possible values of $\lambda$:

[TABLE]

where a signal-to-noise ratio for the force $\mathbf{\zeta}_{a}\equiv\mathbf{a}\Delta t/\sqrt{V}$ was introduced. Figure 1B plots an example force posteriors obtained with the marginalized method and with fixed-$\lambda$ inference schemes. The wider marginalized method posterior takes into account all possible $\lambda$ values. Appendix A5 demonstrates that in contrast to the fixed-$\lambda$ posteriors, the marginalized posterior is in general non-symmetric.

3.6 Numerical results

The performance of the marginalized method was investigated on simulated trajectories. Random trajectories were simulated in a 2D box with periodic boundary conditions, a uniform total force, and a triangular diffusivity profile along the $x$ axis (Fig. 2A,B). Other simulation parameters are given in Appendix A6. For each trial, the simulated trajectories were then analyzed using the TRamWAy software platform [72] and following a procedure similar to the one used in [19, 61] and consisting of (i) individual spatial tessellation in each trial; (ii) assignment of recorded displacements to spatial domains; (iii) inference of $\mathbf{\zeta}_{\mathrm{sp}}$ and $\mathbf{\zeta}_{\mathrm{t}}$ in each domain; (iv) calculation of the Bayes factor in each domain.

The marginalized Bayes factor $\hat{K}^{\mathrm{M}}$, inferred in each domain, was then plotted against its expected value $K^{\mathrm{M}}$ to test the accuracy of the method (Fig. 2C,D). The figure shows good correspondence between the inferred Bayes factor and the expected Bayes factor. 95 % confidence intervals (CIs) show the extent of the deviation of the results from the true values due to the stochastic nature of the simulated trajectories.

3.7 Microscopic model of heterogeneous diffusivity

The next simulation was performed with two goals: (i) to illustrate how spurious forces may originate from crowding at the molecular scale, and (ii) to illustrate a case, wherein our developed statistical test successfully indicates the absence of non-diffusive forces. For this purpose, we simulated free diffusion of particles with no microscopic drift within a square region with periodic boundaries and with impenetrable immobile beads evenly spaced on a square lattice (Fig. 3A), similar to schemes suggested in[49, 32, 10]. The microscopic diffusivity of the particle was the same throughout the system. A spatial variation in the radii of the immobile beads created a spatial variation in the effective diffusivity on a much larger “mesoscopic” scale, where each analysis bin included $\sim 100$ small beads (Fig. 3B). As a result, recordings at the mesoscopic scale exhibit a diffusivity gradient (Fig. 3C), which contributes to the drift observed on the same scale (Fig. 3D). Note that at long time scales, the system is in physical equilibrium, although particles experience a stationary non-zero drift. Simulation details and parameters are provided in Appendix A7.

The diffusivity gradient contribution to the drift is the spurious force, its exact value depends on $\lambda$. Assuming the value of $\lambda$ is unknown, one can use the Bayes factor test developed above to estimate the a posteriori likelihood of that the observed drift is due to non-spurious forces (Fig. 3E). In our simulation, the inferred Bayes factors were small ($\log_{10}K<-1$) in most parts of the region, supporting the claim that only spurious forces were present (Fig. 3F). Statistical noise in several bins resulted in weaker evidence, which did not let us draw statistically significant conclusions in those zones. These results confirm the capacity of the method to detect spurious forces. Its capacity to detect non-spurious forces will be illustrated in the next section.

This simulation captures one possible microscopic mechanism behind the observation of a diffusivity gradient on the mesoscopic scale in biological systems. However, note that the homogeneous composition required for a uniform microscopic diffusivity is probably achievable in the biological systems only on the molecular scale ($10^{-9}\,\mathrm{m}$ and smaller). On this scale though, it is not clear whether the diffusivity itself is well-defined, since by definition it is the result of the action of millions of individual molecules and Fick’s law describes an intrinsically mesoscopic phenomenon.

Other microscopic mechanisms for the diffusivity gradient include (i) confinement, wherein it was shown that the diffusivity in a homogeneous system changes with the distance to a wall [7, 44, 46, 78], (ii) corralled motion [57], (iii) hydrodynamic coupling to other objects in the medium [16, 3], (iv) temperature gradients [18, 9, 81], or (v) intermittent trapping [54].

4 Applications

The developed method was tested on two experimental systems. The first one was a well-controlled setup of a bead in the optical tweezers. The second one was a complex biological process of HIV virion assembly in a T cell [25], where the OLE is potentially only an approximation to the true biomolecule dynamics (ignoring inertial effects, colored noise or memory of the previous states).

4.1 Optical tweezers

Optical tweezers combine physical trapping of the bead with local laser heating of the medium, leading to a heterogeneous diffusivity field. Therefore, the heating effect and the ensuing spurious forces may interfere with the inferred trapping potential. Figure 4A–C compares the results of Bayes factor calculations for the same system subjected to three different laser powers. The tessellation procedure was designed to assign the same number of jumps to each domain. In all 3 cases, the particle is confined and the Bayes factor favors the presence of forces ($\log_{10}K^{\mathrm{M}}>1$) in a large number of domains, which form a connected region. With the decrease of the laser power, the confinement at the center of the trap becomes more shallow, so that the statistical test only detects confining forces on the trap border.

4.2 Assembly of HIV-1 Virus-Like Particles

The HIV virus-like-particle (VLP) assembly experiments that provided the data are described in [23]. The VLPs derive from the human immunodeficiency virus type 1 (HIV-1), but are immature and deprived of envelope proteins. One of their main components is the group-specific antigen (Gag) protein. It is a viral structural protein produced by the virus that anchors and oligomerizes at the plasma membrane of the host T cells, eventually assembling into a VLP [25]. In the experiments, the HIV-1 Gag precursor was genetically modified to contain a photoactivable fluorescent tag mEOS2 protein. It allowed to record VLP assembly in human CD4 T cells using single-particle tracking photoactivated localization microscopy (sptPALM) [51, 23]. Several VLPs can assemble in parallel in the same observation region.

The TRamWAy software platform was used to tessellate the observation region and infer maps of diffusivities (Fig. 4G) and drift [19, 72, Fig. 4I, see]. The Bayes factor map was then computed. The localization uncertainty was $\sigma_{L}=30\,\mathrm{nm}$, requiring corresponding corrections to the Bayes factor (Appendix A4). Inference results and Bayes factors for a $2\,\mathrm{\mu m}\times 2\,\mathrm{\mu m}$ zone of one T-cell membrane are shown in Fig. 4D-F.

Plots of the trajectories, the density of the recorded points and the diffusivity (Fig. 4D,E,G) indicate that there are two regions of interest (ROI) in the data set. However, the plots of the diffusivity gradient and the drift (Fig. 4H,I) suggest that the two parameters are of the same scale, hence it is not a priori clear, whether the localization of the particles is due to non-diffusive or spurious forces. Only the calculation of the Bayes factors for these regions allowed us to confirm that it is not solely due to heterogeneities in the diffusivity but rather to non-spurious forces ($\log_{10}(K)\gg 1$, Fig. 4F).

Some other individual domains in Fig. 4F bear evidence of a force with rather high Bayes factors ($\log_{10}(K)\geq 1$). In such a complex system, the high $K$ values in these individual domains may stem from local membrane activity, failed capsid assembly [23] or be false detection. In the rest of the region, the Bayes factor is $|\log_{10}(K)|<1$ meaning neither of the models is favored at the chosen level of statistical significance.

As demonstrated in the simulation of Sect. 3.7, the results of any inference procedure depend on how the spatial scale on which the analysis is performed, corresponds to the internal scale of the observed system. An illustration of this fact for the VLP data can be seen in Fig. 5. Here, our statistical test was applied to the same VLP data set on three different spatial scales. When the bins are much larger than the typical structures present in the biological system (Fig. 5A, $0.5\,\mathrm{\mu m}$), the interactions are averaged out and our statistical test confirms the absence of interactions or is inconclusive. On the scale of the structures (Fig. 5B, $0.25\,\mathrm{\mu m}$), one may identify the potential regions of interest, but is unable to resolve their internal structure. When the bins are smaller than the regions of interest (Fig. 5C, $0.05\,\mathrm{\mu m}$), given enough data, the internal structure of the regions can be resolved. At even smaller scales (not shown), when few points are available per bin, one starts losing the connectivity of the regions of interest, and the statistical tests becomes inconclusive or (by design) favor the model with only spurious forces ($H_{0}$). We suggest that one should aim for a scale smaller than the scale of the analyzed structure, but maintain enough points per bin to reach statistically significant conclusions.

5 Discussion

In this paper, we introduced a method to address the inverse problem for the spatially heterogeneous OLE that is robust in regards to changes in the spurious-force contribution. We leveraged Bayesian inference and Bayesian model comparison to account for the uncertainty in the values of the spurious force caused by a heterogeneous diffusivity field. The method provides a test for the presence of non-diffusive forces and returns the values of the non-diffusive forces and diffusivity.

The marginalized posterior takes into account the error in the inferred forces due both to stochastic errors and to possible spurious forces when the true value of $\lambda$ is unknown. The expression for the Bayes factor was derived in a closed form, allowing for identification of natural parameters associated with the dynamics, namely, the signal-to-noise ratios for the total force (drift) and spurious forces, $\mathbf{\zeta}_{\mathrm{t}}$ and $\mathbf{\zeta}_{\mathrm{sp}}$, and the relative strength of the localization uncertainty $4\sigma_{L}^{2}/(nV)$. Interestingly, we showed that under some configurations, the discrimination between active and spurious forces is impossible without introducing additional assumptions.

As for any statistical method, a prerequisite for our method is that one observe the trajectories on the “right” spatial and temporal scales, which depend on the individual system. In particular, the spatial tessellation employed here should be constructed on the appropriate spatial scale, i.e. finer than the spatial heterogeneities of interest and coarse enough to provide sufficient measurements in each mesh domain (as discussed above). Another condition required by our method is that the number of points per bin be $>4$ in 1D and $>3$ in 2D, which are the equivalent numbers of points contained in the prior. Otherwise, due to the choice of $\mu_{\pi}$, the model with only spurious forces is likely to always be favored. Our experience with the method indicates that for the biological data we tested, $n\geqslant 20$ typically provides a reasonable compromise between the spatial and statistical resolution.

The VLP example demonstrated successful utilization of the method for the detection of biological activity. The test was applied in an unsupervised way, which makes it useful for automatic analysis of single-molecule dynamics. In general, however, the results may depend on the spatial meshing. For the VLP data set, we had the advantage of a priori knowing the characteristic spatial scale of the underlying biological processes [23]. In a general case, one may need to sample multiple spatial scales in an attempt to optimize the detection. An optimal mesh scale in this case can be seen as a trade-off between increasing statistical significance (by getting more data per domain) and increasing resolution (by reducing the domain size).

Potential ways to circumvent this fundamental trade-off of spatial versus statistical resolution could be to regularize the inference of the diffusivity and drift fields [19, 47] or to cluster the regions with similar Bayes factor values based on a certain rule. However, the former approach induces correlations between the results inferred in different domains making analytical calculations intractable and hindering interpretation of the results. The main difficulty with the latter approach consists in defining the appropriate clustering criterion and in accounting for how the uncertainty in the individual Bayes factors propagates to the Bayes factors of the clusters.

One should also keep in mind that the validity of the main result (12) relies on the assumption that the diffusivity $b$ is smooth enough, so that the gradient $\nabla b$ exists on the spatial scale on which the system is experimentally probed. Additionally, we stress that $\bm{\alpha}$ and $b$ are mesoscopic quantities and their values may change depending on the analyzed scale, [82, see]. In practice, the choice of the spatial and temporal resolutions for the analysis is limited by the particular experimental setup and the properties of the biological system.

The Bayesian approach that we proposed here is general and not limited to the OLE equation. The ambiguity of the stochastic integration is encountered in numerous other scientific fields involving stochastic equations with multiplicative heterogeneous noise. The effect is usually ignored and an arbitrary standard convention is used. The marginalized method allows us to avoid arbitrarily choosing an integral convention in the absence of system-specific information, therefore providing more robust results.

The marginalized method code is available as a module of the open source project TRamWAy [67] and the microscopic crowding simulation code is available at [88]. Two Jupyter notebooks are provided as illustration of the module interface [72].

Acknowledgments. We thank Aleksandra Walczak, Vincent Hakim, Bassam Hajj, Mathieu Coppey and Maxime Dahan (deceased) for helpful discussions. This study was funded by the Institut Pasteur, L’Agence Nationale de la Recherche (TRamWAy, ANR-17-CE23-0016), the INCEPTION project (PIA/ANR-16-CONV-0005, OG), and the programme d’investissement d’avenir supported by L’Agence Nationale de la Recherche ANR-19-P3IA-0001. The funding sources had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Conflicts of interest. The authors declare to have no financial or non-financial conflicts of interest.