# Optimal estimates of the diffusion coefficient of a single Brownian trajectory

###### Abstract

Modern developments in microscopy and image processing are revolutionizing areas of physics, chemistry and biology as nanoscale objects can be tracked with unprecedented accuracy. The goal of single particle tracking is to determine the interaction between the particle and its environment. The price paid for having a direct visualization of a single particle is a consequent lack of statistics. Here we address the optimal way of extracting diffusion constants from single trajectories for pure Brownian motion. It is shown that the maximum likelihood estimator is much more efficient than the commonly used least squares estimate. Furthermore we investigate the effect of disorder on the distribution of estimated diffusion constants and show that it increases the probability of observing estimates much smaller than the true (average) value.

###### pacs:

05.40.Jc, 31.15.xk, 87.16.dp, 61.43.Er## I Introduction

Single particle tracking dates back to the classic study of Perrin on Brownian motion (BM) perrin (). It generates the position time series of an individual particle trajectory in a medium (see, e.g., Refs. bra (); saxton ()) and when properly interpreted, the information drawn from a single, or a finite number of trajectories, can provide insight into the mechanisms and forces that drive or constrain the motion of the particle. The method is thus potentially a powerful tool to probe physical and biological processes at the level of a single molecule moerner (). At the present time, single particle tracking is widely used to characterize the microscopic rheological properties of complex media mason (), and to probe the active motion of biomolecular motors greenleaf (). In biological cells and complex fluids, single particle trajectory (SPT) methods have, in particular, become instrumental in demonstrating deviations from normal BM of passively moving particles (see, e.g., Refs. golding (); weber (); bronstein (); seisenberger ()).

The reliability of the information drawn from SPT analysis, obtained at high temporal and spatial resolution but at expense of statistical sample size is not always clear. Time averaged quantities associated with a given trajectory may be subject to large fluctuations among trajectories. For a wide class of anomalous diffusions described by continuous-time random walks, time-averages of certain particle’s observables are, by their very nature, themselves random variables distinct from their ensemble averages rebenshtok (). An example is the square displacement time-averaged along a given trajectory, which differs from the ensemble averaged mean squared displacement he (). By analyzing time-averaged displacements of a particular trajectory realization, subdiffusive motion can actually look normal, although with strongly differing diffusion coefficients from one trajectory to another lubelski ().

Standard BM is a much simpler and exceedingly well-studied random process yuval () than anomalous diffusion, but still it is far of being as straightforward as one might be tempted to think. Even in bounded systems, despite the fact that the first passage time distribution has all moments, first passages to a given target of two independent identical BMs, starting at the same point in space, will most likely occur at two distinctly different time moments carlos (), revealing a substantial manifestation of sample-to-sample fluctuations. Ergodicity, that is, equivalence of time- and ensemble-averages of square displacement holds only in the infinite sample size limit. In practice, it means that standard fitting procedures applied to finite (albeit very long) trajectories of a given particle will unavoidably lead to fluctuating estimates of the diffusion coefficient . Indeed, variations by orders of magnitude have been observed in SPT measurements of the diffusion coefficient of the LacI repressor protein along elongated DNA austin () (see also Section VI.1). Significant sample-to-sample fluctuations resulting in broad histograms for the value of the diffusion coefficient have been observed experimentally for two-dimensional (2D) diffusion in the plasma membrane saxton (), as well as for diffusion of a single protein in the cytoplasm and nucleoplasm of mammalian cells goulian ().

Such a broad dispersion of the value of the diffusion coefficient extracted from SPT measurements, raises important questions about the correct or optimal methodology that should be used to estimate . Indeed, these measurements are performed in rather complex environments and each SPT has its own history of encounters with other species, defects, impurities, etc., which inevitably results in rather broad histograms for observed . On the other hand, it is highly desirable to have a reliable estimator of the diffusion coefficient even for the hypothetical ”pure” cases, such as, e.g., unconstrained standard BM. A reliable estimator should produce a distribution of as narrow as possible and with the most probable value as close as possible to the ensemble average one. A knowledge of the distribution of such an estimator could provide a useful gauge to identify effects of the medium complexity as opposed to variations in the underlying thermal noise driving microscopic diffusion. Commonly used methods of extraction of from the SPT data are based on a least square (LS) estimate of the time-averaged square displacement and some of its derivatives (see, e.g., saxton (); goulian (); saxton2 () and the next section). A recent study, Ref. boyer (), focussed on estimators for for 1D BM, the statistics of which is amenable to analytical analysis. Several methods for estimating from the SPT data were studied and it was shown that a completely different approach - consisting of maximizing the unconditional probability of observing the whole trajectory - is superior to those based on the LS minimization. As a matter of fact, at least in 1D systems the distribution of the maximum likelihood (ML) estimator of the diffusion coefficient not only appears narrower than the LS ones, resulting in a smaller dispersion, but also the most probable value of the diffusion coefficient appears closer to the ensemble average boyer ().

In this paper we focus first on the case of pure standard BM and calculate exactly, for arbitrary spatial dimension , the distribution of the maximum likelihood estimator

(1) |

of the diffusion coefficient of a single BM trajectory . The parameter here is the lag time (at which the measurement is started) which can be set equal to zero for standard BM without any lack of generality. However for anomalous diffusion, or BM in presence of disorder will play a significant role. The symbol denotes ensemble average, so that

(2) |

being the ensemble-average diffusion coefficient. Consequently, the random variable is defined as the ratio of the realization-dependent diffusion coefficient, calculated as the weighted time-average along a single trajectory, and the ensemble average diffusion coefficient. Clearly, .

Further on, we analyze here a useful measure of sample-to-sample fluctuations - the distribution function of the random variable

(3) |

where and are two identical independent random variables with the distribution . Hence, the distribution probes the likelihood of the event that the diffusion coefficients drawn from two different trajectories are equal to each other.

Finally, we discuss the effect of disorder on the distributions and for 1D BM in random media. We consider two different models of diffusion in 1D random environments - diffusion in presence of a random quenched potential with a finite correlation length, as exemplified here by the Slutsky-Kardar-Mirny model slutsky (), and diffusion in a random forcing landscape - the so-called Sinai model Sinai1982 (). The former is appropriate for diffusion of proteins on DNA, which is affected by the base-pair reading interaction and thus is sequence dependent, while the latter describes, for example, the dynamics of the helix-coil boundary in a melting heteropolymer pgg (). Note that in the former case, at sufficiently large times, one observes a diffusive-type motion with , while in the latter case dynamics is strongly anomalous so that is logarithmically confined, .

The paper is outlined as follows: In Section II we recall some common fitting procedures used to calculate the diffusion coefficient from single particle tracking data. In Section III we focus on the maximum likelihood estimator and, generalizing the approach developed in Ref. boyer () for 1D systems, obtain new results for the moment generating function and the probability density function of the ML estimator for arbitrary spatial dimension . In that Section we also obtain the asymptotical behavior of the probability distribution function , as well as its kurtosis and skewness. Next, in Section IV we focus on the probability distribution function of the random variable - a novel statistical diagnostics of the broadness of the parental distribution which probes the likelihood of the event that two estimates of the diffusion coefficient drawn from two different trajectories are the same. Further on, Section V presents a comparison of the commonly used least squares estimator and the maximum likelihood estimator. We show that the latter outperforms the former in any spatial dimension producing a lower variance and the most probable value being closer to the ensemble average value. Next, in Section VI we focus on Brownian motion in presence of disorder. As exemplified by two models of dynamics in systems with quenched disorder - Sinai diffusion (random force) and Slutsky-Kardar-Mirny model (random potential), disorder substantially enhances the importance of sample-to-sample fluctuations. We show that the observation of values of the diffusion coefficient significantly lower than the ensemble average becomes more probable. We show, as well, that as the strength of disorder is increased, the distribution undergoes a surprising shape-reversal transition from a bell-shaped unimodal to a bimodal form with a local minimum at . Finally, we conclude in Section VII with a brief recapitulation of our results and some outline of our further research.

## Ii Fits for the diffusion coefficient of a single trajectory

To set up the scene, we first briefly recall several fitting procedures commonly used to calculate the diffusion coefficient from the SPT data. More detailed discussion can be found in Refs. saxton (); goulian (); saxton2 (); boyer (). We focus here on estimators which yield a first power of . Non-linear estimators, e.g. a mean maximal excursion method ralf () which has been used to study anomalous diffusion and produces , will be analyzed elsewhere.

One of the simplest methods consists in calculating a least squares estimate based on the minimization of the integral

(4) |

where the diffusion law is taken either as a linear, , or an affine function, . In particular, for the linear case the least squares minimization yields the following linear-least-squares estimator:

(5) |

where is the normalization factor, , conveniently chosen so that .

A second, more sophisticated, approach is based on

(6) |

which is the temporal moving average over a sufficiently long trajectory produced by the underlying process of duration . The diffusion coefficient is then extracted from fits of , or from a related least squares estimator, which is given by the following functional of the trajectory

(7) |

where is the same normalization constant as in Eq. (5). Note that the random variable is again conveniently normalized so that , which enables a direct comparison of the respective distributions of different estimators. As shown in Ref. boyer (), only provides a slightly better estimate of than . A conceptually different fitting procedure has been discussed in Ref. boyer () which amounts to maximizing the unconditional probability of observing the whole trajectory , assuming that it is drawn from a Brownian process with mean-square displacement (see Eq. (2)). This is the maximum likelihood estimate which takes the value of that maximizes the likelihood of , defined as:

(8) |

where the trajectory is appropriately discretized. Differentiating the logarithm of with respect to and setting , one finds the maximum likelihood estimate of , which upon a proper normalization is defined by Eq. (1).

Below we will derive the distribution function of the ML estimator and compare it against numerical results for the distribution function of the LS estimator for .

## Iii Distribution of the ML estimator

### iii.1 The moment generating function

Let denote the moment generating function of the random variable defined in Eq. (1),

(9) |

The squared distance from the origin, of -dimensional BM at time for a given realization, decomposes into the sum

(10) |

being realizations of trajectories of independent 1D BMs (for each spatial direction). Thus, factorizes

(11) |

where

(12) |

Here, in order to calculate , we follow the strategy of Ref. boyer () and introduce an auxiliary functional:

(13) |

where the expectation is for a BM starting at at time . We derive a Feynman-Kac type formula for considering how the functional in Eq. (13) evolves in the time interval . During this interval the BM moves from to , where is an infinitesimal Brownian increment such that and , where denotes now averaging with respect to the increment . For such an evolution we have to order :

(14) | |||||

Expanding the right-hand-side of the latter equation to second order in , linear order in and performing averaging, we find eventually the following Schrödinger equation:

(15) |

The solution of this equation has been obtained in Ref. boyer () and gives

(16) |

where is the modified Bessel function abramowitz (). Consequently, we find the following general result

(17) |

Note that is independent of and , as it should in virtue of the scaling properties of the BM.

### iii.2 The distribution function

We turn next to the analysis of the distribution of the ML estimator defined in Eq. (1). First of all, we calculate several first moments of by merely differentiating the result in Eq. (17):

(18) |

Consequently, one may expect that all moments tend to as , so that . For fixed , the variance , the coefficient of asymmetry and the kurtosis . All these characteristics vanish when .

Note next that since , the poles of are located at , where is the th zero of the Bessel function abramowitz (). Consequently, for even , we can straightforwardly find in form of an infinite series in the zeros of the Bessel function . For , has only simple poles so that the expansion theorem gives

(19) |

For and the standard residue calculus yields

(20) |

and

(21) |

where

(22) |

Similar results can be readily obtained for greater even .

For arbitrary , including odd values, the distribution is defined by inverting the Laplace transform and is given by the following integral:

(23) |

where the phase is given by

(24) |

and being the th order Kelvin functions abramowitz ().

Finally, we consider the small- and large- asymptotic behavior of the probability density function . To extract the small- asymptotic behavior of we consider the large- form of . From Eq. (17) we get

(25) |

as . Consequently, we find the following strongly non-analytical behavior:

(26) |

The large- behavior of the distribution is defined by the behavior of the moment generating function in the vicinity of ,

(27) |

Consequently, we find that for , decays as

(28) |

This behavior is, of course, consistent with the series expansions in Eqs. (19), (20) and (21). Our results on the distribution are summarized in Fig. 1.

## Iv The distribution of the random variable

Suppose next that we have two different independent realizations of BM trajectories, and which we use to generate to independent random variables and . A natural question arising about their suitability as estimators is how likely is it that they will have the same value? Of course the distributions and thus moments of these two random variables are the same, however a measure of their relative dispersion can be deduced by studying the distribution function of the random variable carlos (); carlos1 (), defined in Eq. (3). This distribution is given explicitly by gleb ()

(29) |

and hence, it suffices to know in order to determine .

For , (and, in fact, for any other even ), can be evaluated exactly. Plugging Eq. (19) into (29), we get:

(30) |

Performing the sum over , we arrive at the following result for the distribution in 2D systems

(31) |

Numerically obtained distributions for and dimensional systems are presented in Fig. 2. Notice that in all dimensions is the most probable value of the distribution so that most probably . Nevertheless, the distributions are rather broad which signifies that sample-to-sample fluctuations are rather important.

## V Comparison of the LS and ML estimators.

We will now show that the ML estimator defined in Eq. (1) substantially outperforms the MS estimator as defined in Eq. (7) in any spatial dimension . This is a very surprising result as one would intuitively expect, and it is often stated in the literature, that using the process has the effect of a reducing the fluctuations of the estimate of because the process is partially averaged in time.

To demonstrate this, we present in Fig. 3 a comparison of the analytical results for of the ML estimator with the corresponding distributions of the LS estimator for obtained numerically.

Indeed, we find that the variance of the distribution equals , and for and , respectively. The distribution of the ML estimator appears to be substantially narrower so that the variance is significantly lower, , and . Moreover, the most probable values of are closer to the ensemble average value than the most probable values of to : we observe that the distribution attains its maximal values at , and for and , respectively, while the corresponding maxima of the distribution are located at and . Last but not least, the distribution appears to be significantly broader than , as revealed by Fig. 2. The worst performance of the LS estimator is in 1D systems in which the distribution has a bimodal -shape with a local minimum at , and maxima (most likely values) around and . This means that the values and drawn from two different trajectories will most probably be different by an order of magnitude!

## Vi 1D Brownian motion in presence of disorder.

In this final section we address the question of how the distribution of the ML estimator of a single trajectory diffusion coefficient will change in presence of quenched disorder. We will consider two different models of BM in random 1D environments - diffusion in presence of a random correlated potential and diffusion in presence of a random force.

### vi.1 Diffusion in presence of a random potential

First we consider a BM in a 1D inhomogeneous energy landscape, where disorder is correlated over a finite length . This model gives a simple description of diffusion of a protein along a DNA sequence, for instance, where the particle interacts with several neighboring base pairs at a time slutsky (). The total binding energy of the protein is assumed to be a random variable. When the particle hops one neighboring base further to the right or to the left, its new energy is highly correlated to the value it had before the jump. Slutsky et al. slutsky () modeled this process as a point-like particle diffusing on a 1D lattice of unit spacing with random site energies , whose distribution is Gaussian with zero mean, variance and is correlated in space as . At each time step, the particle located at some site jumps to the left or to the right with probabilities and , respectively, where . Diffusion is asymptotically normal for any disorder strength . Nevertheless, the particle can be trapped in local energy minima for long periods of time. During an extended intermediate time regime, it is observed that first passage properties fluctuate widely from one sample to another slutsky ().

Our numerical simulations reveal that disorder has a dramatic effect on the distributions and . As shown by Fig. 4, the distribution broadens significantly in the small regime: very small values of the time-average diffusion constant (compared to the thermal and disorder average) become increasingly more probable as the disorder strength increases. However, the right tail of is much less affected. Similarly, two independent measurements are likely to differ significantly, even in moderately disordered media (see inset of Fig. 4). When , the distribution undergoes a continuous shape reversal transition - from a unimodal bell-shaped form to a bimodal -shape one with the minimum at and two maxima approaching and at larger disorder strengths. Unfortunately, it does not seem possible to obtain this critical value analytically. Even for the case of a pure Brownian motion considered in Ref. carlos (), such an analysis appears to be extremely difficult.

Therefore, for sample-to-sample fluctuations becomes essential and it is most likely that the diffusion coefficients drawn from two different trajectories will be different.

### vi.2 Diffusion in presence of a random force

We discuss now the effect of disorder on the distributions and for 1D BM in presence of a quenched uncorrelated random force - the so-called Sinai diffusion Sinai1982 (). In this model, one considers a random walk on a 1D infinite lattice and site dependent hopping probabilities: for hopping from to the site and for hopping to the site . Here, are independent, uncorrelated, identically distributed random variables with distribution and the strength of disorder is bounded away from and . It is well-known that in the large- limit the model produces an anomalously slow sub-diffusion , where the angle brackets denote averaging with respect to different realizations of disorder. At shorter times, however, one observes an extended stage with a transient behavior which is substantially different from the asymptotic one. As a consequence, the statistics of , defined in equation (1), and will depend not only on the integration time but also on the lag time from which a single trajectory is analyzed.

We have numerically computed the distributions and . In Fig. 5 we present the dependence of the and statistics on the strength of the disorder . As in the previous disordered potential case, we find that the maximum of shifts toward zero as the disorder gets stronger. For comparison, the solid black line in Fig. 5 represents observed for standard BM, . Moreover, the stronger the disorder is, the broader the distribution becomes, yielding more peaked maxima in (see the inset of Fig. 5). We also note that has a bimodal -shaped form even for the weakest disorder that we have considered, suggesting that the zero-disorder limit is non-analytic compared to the continuous transition observed for diffusion in the random energy landscape of the previous section.

In Fig. 5 and Fig. 5, we show the statistics of for different values of and . Increasing (or ) we observe that changes from an almost uniform form, for which any relation between and is equally probable, to a bimodal -shaped distribution, which signifies that in this regime two ML estimates and will most likely have different values. Bimodality is a property of the Sinai regime carlos (), as also noticed earlier in Ref. redner (), and thus it shows up only at sufficiently long times when the trajectories follow the asymptotic ultra-slow diffusion. The distribution is remarkably sensitive to the characteristic aging of Sinai diffusion.

As a final observation, one may also study the statistics of for trajectories evolving in the same random force field. In this case (not shown here) one gets a narrow distribution for of and a unimodal that converges to a delta singularity at when the disorder becomes infinite. This is due to the fact, known as Golosov phenomenon, that two trajectories in the same disorder will move together Golosov1983 ().

## Vii Discussion

We have analyzed the reliability of the ML estimator for the diffusion constant of standard Brownian motion and shown its superiority over the more commonly used LS estimator in a number of important aspects, notably the variance of the estimator, the proximity of the most probable value to the true mean value and the distribution of the random variable which is a measure of the extent to which two estimations of vary. Going beyond the important test case of pure Brownian motion we have also analyzed the effect of quenched disorder, modeling fluctuations of the local energy landscape and forces. As one may have intuitively expected, the presence of short range disorder tends to broaden the distribution of the so measured value of , as it presents an additional source of fluctuation. However in the Sinai model, in the same realization of the force field, trajectories are disorder dominated and are almost independent of the thermal noise, leading to highly peaked distributions of . Analytically understanding the distribution of in the presence of disorder presents an interesting mathematical challenge which will involve analysis of the corresponding Schrödinger equation with a random drift term. Further interesting questions remain to be addressed, in particular can one use the two time correlation function of measured trajectories to obtain better estimators ? Single particle tracking technology will undoubtedly further improve in the coming years and many interesting mathematical, statistical and physical challenges will arise in the ultimate goal of getting the most out of the trajectories so obtained.

###### Acknowledgements.

DSD, CMM and GO gratefully acknowledge support from the ESF and hospitality of NORDITA where this work has been initiated during their stay within the framework of the ”Non-equilibrium Statistical Mechanics” program.## References

- (1) J. Perrin, C. R. Acad. Sci. 146, 967 (1908); Ann. Chim. Phys. 18, 5 (1909).
- (2) C. Bräuchle, D. C. Lamb, and J. Michaelis, Eds., Single particle tracking and single molecule energy transfer (Wiley-VCH, Weinheim, 2010).
- (3) M. J. Saxton and K. Jacobson, Ann. Rev. Biophys. Biomol. Struct. 26, 373 (1997).
- (4) W. E. Moerner, Proc. Natl. Acad. Sci. USA 104, 12596 (2007).
- (5) T. G. Mason and D. A. Weitz, Phys. Rev. Lett. 74, 1250 (1995).
- (6) W. J. Greenleaf, M. T. Woodside, and S. M. Block, Annu. Rev. Biophys. Biomol. Struct. 36, 171 (2007).
- (7) I. Golding and E. C. Cox, Phys. Rev. Lett. 96, 098102 (2006).
- (8) S. C. Weber, A. J. Spakowitz, and J. A. Theriot, Phys. Rev. Lett. 104, 238102 (2010).
- (9) I. Bronstein et al., Phys. Rev. Lett. 103, 018102 (2009).
- (10) G. Seisenberger et al., Science 294, 1929 (2001).
- (11) A. Rebenshtok and E. Barkai, Phys. Rev. Lett. 99, 210601 (2007).
- (12) Y. He, S. Burov, R. Metzler and E. Barkai, Phys. Rev. Lett. 101, 058101 (2008).
- (13) A. Lubelski, I. M. Sokolov, and J. Klafter, Phys. Rev. Lett. 100, 250602 (2008).
- (14) P. Mörters and Y. Peres, Brownian motion (Cambridge University Press, Cambridge, 2010).
- (15) C. Mejía-Monasterio, G. Oshanin, and G. Schehr, J. Stat. Mech. P06022 (2011)
- (16) Y. M. Wang, R. H. Austin and E. C. Cox, Phys. Rev. Lett. 97, 048302 (2006).
- (17) M. Goulian and S. M. Simon, Biophys. J. 79, 2188 (2000).
- (18) M. J. Saxton, Biophys. J. 72, 1744 (1997).
- (19) D. Boyer and D. S. Dean, J. Phys. A: Math. Gen. 44, 335003 (2011).
- (20) M. Slutsky, M. Kardar, and L. A. Mirny, Phys. Rev. E 69, 061903 (2004).
- (21) Ya. G. Sinai, Theory Probab. Appl. 27, 256 (1983).
- (22) P. G. de Gennes, J. Stat. Phys. 12, 463 (1975).
- (23) V. Tejedor et al., Biophys. J. 98, 1364 (2010).
- (24) M. Abramowitz and I. R. Stegun, Handbook of mathematical functions (Dover, New York, 1972).
- (25) C. Mejía-Monasterio, G. Oshanin, and G. Schehr, Phys. Rev. E 84, 035203(R) (2011).
- (26) G. Oshanin, Yu. Holovatch, and G. Schehr, Physica A 390, 4340 (2011).
- (27) G. Oshanin and S. Redner, Europhys. Lett. 85, 10008 (2009).
- (28) A. O. Golosov, Soviet Math. Dokl. 28, 18 (1983).