Time-averaged MSD of Brownian motion

Time-averaged MSD of Brownian motion

Alexei Andreanov and Denis S. Grebenkov Abdus Salam ICTP - Strada Costiera 11, 34151, Trieste, Italy aandrean@ictp.it Laboratoire de Physique de la Matière Condensée (UMR 7643), CNRS – Ecole Polytechnique, F-91128 Palaiseau, France denis.grebenkov@polytechnique.edu
July 24, 2019

We study the statistical properties of the time-averaged mean-square displacements (TAMSD). This is a standard non-local quadratic functional for inferring the diffusion coefficient from an individual random trajectory of a diffusing tracer in single-particle tracking experiments. For Brownian motion, we derive an exact formula for the Laplace transform of the probability density of the TAMSD by mapping the original problem onto chains of coupled harmonic oscillators. From this formula, we deduce the first four cumulant moments of the TAMSD, the asymptotic behavior of the probability density and its accurate approximation by a generalized Gamma distribution.

02.50.-r, 05.60.-k, 05.10.-a, 02.70.Rr

Keywords: diffusion, Brownian motion, MSD, time average, single-particle tracking

1 Introduction

The statistical properties of Brownian motion play a crucial role in disciplines as various as mathematics, physics, chemistry, biology, engineering, economics and ecology [1, 2, 3, 4, 5]. In many cases, an observation or an experiment can be repeated a large number of times to get a representative statistics of outcomes. In other situations, one measures an observable which reflects a cumulative effect of a large number of molecules (e.g., the macroscopic nuclear magnetic resonance signal formed by an extremely large number of nuclei [5]). In both cases, a system can be characterized by ensemble-averaged quantities. For instance, the ensemble average (or expectation) of the mean-square displacement (MSD) of a freely diffusing particle, , yields a measure of the diffusion coefficient through the Einstein’s relation: . However, there are other situations when few or even single stochastic trajectory is available that makes the ensemble averaging inaccurate or impossible. The common examples are stock prices in finance [6] or trajectories of macromolecules inside living cells acquired by single particle tracking (SPT) techniques [3, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]. In the former case, the observation (i.e., a price time series) is evidently unique for each asset. In the latter case, although an experiment can be repeated, the ensemble average may still be problematic, especially for tracers that move in spatially heterogeneous and time-evaluating media such as living cells. In order to infer the dynamical information from such individual trajectories, one needs therefore to replace ensemble averages by time averages. For instance, the diffusion coefficient is often inferred by considering the time-averaged MSD (TAMSD) at a lag time over a sample trajectory of duration :


Although this is still a random variable, the time average reduces its fluctuations around the mean which, for Brownian motion, is .

The statistical properties of quadratic functionals of a Gaussian process have been intensively studied. In mathematical statistics, several series representations of the probability density of a general quadratic form of a Gaussian process have been proposed, e.g. a mixture of chi-squared distributions [21, 22, 23, 24], a power series expansion [25, 26], and a series expansion over Laguerre polynomials [27, 28, 29] (see also [30, 31] for a review). For biophysical applications, Qian and co-workers analyzed the TAMSD of a discrete random walk which is known to have a Gamma (or chi-squared) distribution only at the smallest lag time [7]. The distribution of diffusion coefficients was also studied via Monte Carlo simulations by Saxton [8, 9] (see also [10, 32] for experimental data). More recently, Boyer and Dean investigated the distribution of several estimators of diffusion coefficients for Brownian motion [33]. They reduced the analysis of least-squares and maximum likelihood estimates to the distribution of a local quadratic functional of Brownian motion. Mapping this problem onto an imaginary time Schrödinger equation, they derived explicit analytical results for the distribution of the diffusion coefficient estimator. In a separate paper, Boyer and co-workers discussed the optimal way to extract diffusion coefficients from single trajectories of Brownian motion and showed the superior efficiency of the maximum likelihood estimator over least-squares estimates [34]. An optimal Bayesian method for quantifying biomolecule diffusivity was presented by Voisinne et al. [35].

A different approach to the analysis of inferring schemes with quadratic functionals was developed by Grebenkov [36, 37]. First, the explicit formulas for the mean and variance of the TAMSD (and another quadratic functional) were derived for a large class of Gaussian processes governed by a generalized Langevin dynamics [36]. The formula for the variance allows one to estimate the spread of measurements and to choose an appropriate sample duration . Second, a matrix representation for the characteristic function of was used in order to analyze the probability density of the TAMSD of a general discrete Gaussian process [37]. Moreover, an empirical approximation for by a generalized Gamma distribution (GGD) was proposed:


with three parameters , and , and is the modified Bessel function of the second kind (in the limit , one retrieves the Gamma distribution ). Although the high accuracy of the GGD approximation (2) was confirmed numerically for a discrete Brownian motion, the theoretical status of this approximation remained unclear. In addition, the dependence of the parameters , and on the lag time was not analyzed.

In this paper, we rigorously characterize the TAMSD of Brownian motion. In Sect. 2, we introduce a functional integral representation for a Laplace transform of the probability density which is a cornerstone for all consecutive derivations. For , it leads to a simple closed formula for . For , an exact representation of involves a determinant of an explicit matrix of size , being an integer part of . In Sect. 3, we discuss some practical consequences of these theoretical results. In particular, we derive the explicit formulas for the first four cumulant moments of the TAMSD, as well as the asymptotic behavior of the probability density at small and large . For instance, we show that decays as at small that justifies the use of a GGD as a convenient approximation for . It also explains a remarkable accuracy of this approximation which correctly reproduces the asymptotic behavior of at both small and large . Finally, we derive the asymptotic relations for the parameters , and of the GGD as functions of the lag time by matching the first moments of the true and approximate distributions. These relations allow one to get an accurate approximation of the probability density for a given lag time that opens a number of practical applications for single-particle tracking experiments. These findings and future perspectives are summarized in Conclusion.

2 Theoretical results

Throughout this section, we consider to be Brownian motion with mean zero and covariance , with diffusion coefficient set to . The sample duration is fixed to , while the lag time is dimensionless. We are interested in the probability density of the random variable which denotes the TAMSD of :


We aim at calculating the Laplace transform of the probability density


Once is obtained, one can retrieve the probability density either through the inverse Laplace transform of , or through the inverse Fourier transform of the characteristic function (although both ways are formally equivalent, Fourier transforms are more convenient for numerical implementation). From a theoretical point of view, the knowledge of is therefore fully equivalent to that of (see Sect. 3 for practical consequences).

We start by expressing the Laplace transform in Eq. (4) as a functional integral over trajectories of a Brownian motion [4, 17]


where we assume that Brownian motion starts from the origin: (anyway, the TAMSD is translation invariant and thus the above functional integral is independent of the starting point). The next step would be to use the Feynman-Kac formula [4, 18] and map the functional integral to a propagator of a certain quantum mechanical problem, to compute the propagator and to evaluate . This method was successfully used by Boyer and Dean for local quadratic functionals of Brownian motion [33]. However, in our problem different parts of the trajectory, at times and , are coupled, and the action is non-local, forbiding a direct use of the Feynman-Kac formalism. To overcome this difficulty, we partition the trajectory into pieces and relabel them as new Brownian walkers in order to reduce it to an effectively many-body but local problem. The latter is not surprising since a system with memory is often equivalent to some many-body system, for example dynamics of a classical spin-glass or quantum mean-field spin-glass are typically mapped onto single particle problems with memory [19, 20]. The transformation appearing in our problem is similar in spirit, but uses the inverse mapping to turn a non-local single particle problem into a local many-body problem. We illustrate this in detail for the case where the mapping and the computation are relatively simple. The case is a straightforward generalization which is treated in the next Subsection.

Figure 1: (Color online) Splitting of a trajectory into blocks. (a) For , there are only three blocks: two coupled parts and and one free part ; (b) When (), a trajectory is split into equal blocks yielding a chain of successively coupled oscillators , …, ; (c) When (for some ), the time interval is covered by subintervals of duration , plus a smaller interval of duration . In order to split the interval regularly, each subinterval is divided into two parts, of durations and , respectively. As a consequence, the trajectory is mapped onto two chains of coupled oscillators and . The two chains and are coupled through boundary conditions for successive parts of the trajectory: , which express continuity of the original trajectory .

2.1 Case

The partitioning of the trajectory in the case is illustrated on Fig. 1a. We split the whole trajectory into three parts:

The original trajectory is then replaced by a collection of two harmonically coupled particles and a free walker. Boundary conditions express continuity of the original trajectory :

Note that the subscripts of , and are just notations for these integration variables. We have defined for convenience

The original non-locality is translated into harmonic coupling of particles , , the overall action is now local. Since there are just two particles involved, this case is easy to resolve: , are decoupled in a standard way by defining their normal modes ,

where represents the center of mass motion and describes harmonic oscillations around the center of mass. Boundary conditions for and are:

Therefore we get a system of two free particles , and one harmonic oscillator :

The Feynman-Kac formula can now be applied to each particle separately expressing as a product of propagators of respective quantum problems:

i.e. two free propagators and one propagator for a harmonic oscillator:


where and . In the limit , one retrieves the diffusive propagator of a free walker:


Pluging expressions for propagators, one gets

where , . Substituting the expressions for , , , and evaluating directly the Gaussian integrals, we get the following closed formula


This formula is one of the main results of the paper.

2.2 Case

Although the analysis for the case relies on the same ideas and it is a straightforward generalization of the above calculus, the computation is much more involved because the number of particles is growing linearly to infinity as . Technically one has to distinguish subcases of for all integer . This is a reflection of the fact that the smaller is, the more images with integer can be placed inside the interval , yielding coupled to , which is coupled to and so forth. As previously, the original problem with self-interaction maps onto an exactly solvable model of coupled harmonic oscillators.

For discrete lag times with , the whole time interval is covered by equal subintervals so that the trajectory can be mapped onto a single chain of coupled harmonic oscillators as shown on Fig. 1b. Using this mapping, we derive in A.1 the exact formula for the Laplace transform :


where , with . The matrix of size , which represents the coupling of harmonic oscillators, is given by Eqs. (48) or (49). Although Eq. (9) provides an exact solution for the original problem, it has a rather formal nature, as evaluation of for large matrix sizes (i.e. small times ) is quite complicated since is a generic matrix with no particular simplifying structure (see A.3). This is the main difficulty encountered in the computation of for small values of .

In the generic case (for some ), the interval is covered by subintervals of duration , plus a smaller interval of duration . In order to split the interval regularly, each subinterval is divided into two parts, of durations and . This covering maps the whole trajectory onto two interacting chains of and harmonic oscillators as illustrated on Fig. 1c. In A.2, we derive the following exact formula


where , with . The matrix of size , which represents the coupling of harmonic oscillators, is given by Eq. (53). Once again, this solution is exact though still quite formal, as the dependence of on and is partly “hidden” in the determinant of the matrix .

3 Discussion

Quite surprisingly, a theoretical analysis of the statistical properties of the TAMSD turned out to be a challenging problem, even in the simplest case of Brownian motion. This reflects the non-locality of the involved quadratic functional and the many-body character of the underlying quantum system. It is also worth stressing that a “resolution” of this problem would mean different things for theoreticians and experimentalists. In the previous Section, we have fully resolved the problem from a theoretical point of view, by providing the exact formulas (8), (9), (10) for the Laplace transform of the probability density. However, this solution remains somewhat formal since an explicit “shape” of the probability density is still unknown, even for the simplest case . In this Section, we discuss some practical consequences of these results.

3.1 Moments

3.1.1 Case .

Using the exact solution (8), one can compute the moments of the TAMSD:


from which the cumulant moments are


Here, is the mean value, the variance, while the skewness and kurtosis are given as and , respectively. We note that the expression for the variance reproduces the result by Qian et al. which was obtained for a discrete random walk [7]. In turn, the formulas for the third and fourth moments are derived for the first time.

Figure 2: (Color online) (Left) The reduced variance , skewness and kurtosis of the TAMSD as functions of the lag time . Symbols present the values obtained from an exact matrix formula for a discrete random walk (see Eq. (6) from [37], the sample length ), while lines refer to the analytical formulas (13), (14), (15) for and (17), (18), (3.1.2) (or (20)) for . (Right) The relative error between the analytical formulas and their counterparts for a discrete random walk. The vertical dotted lines stand at and . All the formulas are accurate for the whole region (an increase of the relative error at small is related to intrinsic differences between discrete and continuous processes).

3.1.2 Case .

For the case (), the explicit formula for is provided in A.3. Taking the derivatives with respect to , we obtain the first cumulant moments:


As the explicit formula (10) for is too cumbersome for , another representation (61) for the cumulant moments was deduced in B. In particular, we checked that Eqs. (17) and (18) for and are also valid for . In turn, the formula for for becomes


Figure 2 shows the reduced variance , skewness and kurtosis of the TAMSD as functions of the lag time . The values obtained from an exact matrix formula for a discrete random walk (see Eq. (6) from [37], the sample length ), are compared to the analytical formulas (13), (14), (15) for and (17), (18), (3.1.2) (or (20)) for . Notice an excellent agreement between the curves.

3.2 Small asymptotic behavior

3.2.1 Case

In the limit of large , Eq. (8) can be approximated as

Using the asymptotic relation for the inverse Laplace transform ,


with and , one derives the asymptotic behavior of the probability density at small :


In other words, the asymptotic decay of as at large implies the sharp decay of as at small . This is the first rigorous justification for using a GGD in Eq. (2) to approximate the probability density .

3.2.2 Case

In this case, we have to analyze the asymptotic behavior of Eqs. (9) or (10). For large , one can replace by (with exponentially small corrections) to get in both cases


where the functions and are explicitly given in A. The asymptotic behavior of at large is mainly determined by the exponential term, while turns out to be an algebraic correction. In C, we derive a representation of in the form , where and become independent of for large (or ). As a consequence, the asymptotic behavior of reads as

Neglecting and using again Eq. (21), one deduces the asymptotic behavior of :


in which is given by Eq. (65). Similar to the case , the probability density sharply decays (as ) when .

3.3 Large asymptotic behavior

Since the inverse Laplace transform can be written as the Bromwich integral in the complex plane, the singularities of for play an important role. In particular, the singularity with the largest (negative) real part, which is called the abscissa of convergence of the Laplace transform, determines the asymptotic behavior of the probability density at large [38].

For , one can easily check that all singularities of are located on the negative real axis. In fact, taking , one can rewrite Eq. (8) as


where is a positive parameter between and (given that ). This function has two sets of singularities at points and which are solutions of the equations

The smallest solution of the second equation can be determined perturbatively in powers of :

The related singularity is thus located at

Note that varies from at to at so that varies from at to at . At the same time, the smallest solution is independent of and is always larger than .

In close vicinity of , one has

from which the inverse Laplace transform yields an asymptotic approximation


with .

The analysis can be extended to the case . For instance, when , one can substitute into Eq. (9) in order to locate the singularities of either as solutions of (with ), or as solutions of . For a general situation for some , the singularties of can be determined in a similar way. In both cases, the large asymptotic behavior is determined by the largest negative root (i.e., the root with the smallest absolute value ) of the equation . Given the square root type of singularity in Eq. (10), the exponential asymptotic decay of at large in the form similar to Eq. (26) is expected to be valid for any lag time .

The accuracy of the asymptotic Eq. (22) for small and Eq. (26) for large is illustrated on Fig. 3. For all lag times, Eqs. (22),(24) correctly describe the asymptotic behavior of at small . Moreover, for larger (e.g., and ), Eq. (22) turns out to reproduce the whole shape of , except for large . The deviations at large are expected because the exponential decay of is not captured in the small asymptotic relation (22).

Figure 3: (Color online) The probability density (solid line), its small asymptotic relation (22) (shown by circles), and large asymptotic relation (26) (shown by crosses), at lag times , and , in semilog (left) and loglog (right) scales.

3.4 Limiting case

In the limit , one gets , from which


We retrieved therefore the limiting behavior of the discrete case for the lag time [7, 37]. Since the density describes the distribution of a squared Gaussian variable , and then Eq. (27) could be directly obtained by integrating the Gaussian probability density for with .

The function diverges at and monotonously decreases on the positive semi-axis , in sharp contract to (with ), which must have a maximum, as being equal to at and . In the limit , the most probable value of (i.e., the position of the maximum) approaches , while the mean value approaches . This also illustrates the fact that has the largest skewness among , as seen on Fig. 2.

3.5 Limiting case

In the opposite limit of it is harder to get analytic behavior from theoretical results, as the limit describes the situation where the number of particles grows to infinity, which, combined with complicated boundary conditions for normal modes, makes a direct evaluation unfeasible. We provide here a few qualitative arguments, while a more rigorous analysis is reported in D. Since the mean and variance of vanish as , it is convenient to rescale by , in order to fix the mean value to . For small , the rescaled TAMSD can be approximated as

where and are independent Gaussian displacements (standing for , with ) with mean zero and variance . The above sum is known to have a chi-squared (or Gamma) distribution with degrees of freedom [39], from which


In the limit , we get and as expected since concentrates around in this limit. Although the limiting results for and are correct, an approximation (28) of for small by a Gamma distribution is not accurate.

3.6 Dimensional units

All the above results were obtained for and , with being a dimensionless lag time between and . Dimensional units can be easily incorporated by replacing by , by and by in the formulas for and . For instance, Eq. (8) reads for a given diffusion coefficient and sample duration as

In the formulas for the cumulant moments, one has to replace by and multiply the right-hand side by , with being the order of the moment.

3.7 Discrete versus continuous cases

As we discussed in Introduction, the TAMSD is commonly employed for inferring the diffusion coefficient from an individual random trajectory of a diffusing tracer. For this purpose, one records positions of a tracer at successive times, with a fixed time step . An individual trajectory of the tracer is therefore discretized, and the TAMSD becomes


where is the position of the tracer at time , is the discrete lag time, and is the number of points in the sample. In the limit , the discrete TAMSD is expected to converge to the continuous one defined by Eq. (1). The probability distribution of the TAMSD for a discrete Gaussian process was investigated in [37]. For a given covariance matrix (and zero mean), one can use the classical matrix representation


where the matrix represents the discrete quadratic form in Eq. (29). When the sample length is smaller than few thousands, the Schur decomposition of the matrix allows one to rapidly compute the Laplace transform or the characteristic function , from which the probability density can be reconstructed through the inverse Fourier transform [37].

Although Eq. (30) provides an exact solution for and an efficient numerical way for computing for moderate , the dependence on the parameters (e.g., the lag time) remains elusive. In contrast, Eq. (8) for Brownian motion is explicit and allows one to derive the asymptotic behavior and the dependence on the lag time. Since approaches in the limit , Eq. (8) is expected to be an accurate approximation to for large enough . This point is illustrated on Fig. 4, in which for is compared to , for different sample lengths . The curves for exhibit deviations at large , while the curves for larger are barely distinguishable from at the present range of values. As a consequence, one can use the explicit formula (8) for discrete samples with . At the same time, the functions and will always be different for large enough . In fact, for a fixed , exhibits a power law asymptotic decay (as is a polynomial in ), while decays exponentially fast. However, the distinction between these behaviors may appear for so large , for which both functions become already extremely small so that a distinction would be irrelevant for practical applications. In summary, Eq. (8) and its extensions for can in practice be used for studying the TAMSD of discretely sampled positions of a freely diffusing particle.

Figure 4: (Color online) (Left) Comparison between for Brownian motion (with ) from Eq. (8) (thick black line) and for discrete random walks with different sample length , with . The curves for are barely distinguishable from the theoretical one. (Right) The relative error between and . As expected, the relative error increases with and it is larger for smaller .

3.8 Relation to generalized Gamma distribution

As shown numerically in [37], the probability density for a discrete random walk can be accurately approximated by a generalized Gamma distribution (GGD) from Eq. (2). Its Laplace transform reads as


On one hand, its comparison to Eqs. (8), (9), (10) immediately shows that a GGD is not an exact distribution for (except for the trivial case ). On the other hand, we have seen earlier that the GGD correctly reproduces the asymptotic behavior at small and large , and it turns out to be a remarkably accurate approximation of (see Fig. 6). Moreover, a simple form of Eq. (2) is particularly attractive for practical purposes (e.g., inferring statistics of diffusion coefficients by maximum likelihood methods, estimating the confidence intervals, etc.). For practical uses of this approximation, it is important to relate the parameters , and to the lag time .

In E, we derive the asymptotic relations for the parameters , and in the limit by matching the first three moments: