Maximum likelihood estimations of force and mobility from single short Brownian trajectories

Abstract

We describe a method to extract force and diffusion parameters from single trajectories of Brownian particles. The analysis, based on the principle of maximum likelihood, is well-suited for out-of-equilibrium trajectories, even when a limited amount of data is available and the dynamical parameters vary spatially. We substantiate this method with experimental and simulated data, and discuss its practical implementation, strengths, and limitations.

---

1 Introduction

Brownian particles are ubiquitous in soft matter and biological sciences, from colloidal particles to fluorescently-tagged molecules. The trajectories of these particles contain precious information about their structure and interactions. For example, particle sizes are routinely quantified by measuring diffusion coefficients of a dilute suspension in a well-characterized fluid.^1,2 Alternatively, when the particles are well-characterized, the trajectory of a Brownian particle can probe the solvent's rheological properties.³ Analysis of Brownian trajectories can also reveal the conservative and dissipative forces acting on particles due to external fields⁴ or interactions with other particles.⁵

Since Brownian trajectories are stochastic, their analysis is necessarily statistical. The physical theory describing the statistics of Brownian particles is well-established.⁶ When particles fluctuate near an equilibrium position, conservative forces acting on them are readily extracted using Boltzmann statistics.⁷ In the absence of external forces, the dissipative forces acting on particles can be calculated from their diffusion coefficient using the Stokes–Einstein relation. More generally, the Smoluchowski equation describes the trajectory of Brownian particles when conservative forces and diffusion coefficients vary over space.⁸ In this case, conservative and dissipative forces can simultaneously be determined from the distribution of displacements at each configuration.

Importantly, this requires the measurement of many trajectories through the same system configuration. This approach has been successfully implemented for micron-sized colloidal particles, which can be repeatedly trapped in the desired configuration and released using optical traps.^5,9–11 However, in some cases Brownian particles cannot be manipulated with optical traps, and in others the laser will perturb the rest of the system. In such situations, one requires a method to determine the forces with less data, typically a single trajectory.

Maximum likelihood estimators have recently been shown to be powerful tools for measuring the diffusion coefficient of Brownian particles.^12–14 They are more efficient and more accurate than the conventional approach, where one determines a diffusion coefficient by calculating a mean-squared displacement and fitting. Given a discretely sampled trajectory, maximum likelihood estimators provide an explicit analytical expression for the best estimate of the diffusion coefficient. This approach explicitly accounts for common experimental artifacts such as localization error from exposure time-induced blurring,¹² finite trajectory length,¹³ or limited photon statistics.¹⁴

We recently introduced a method, based on the same principle, to estimate force and diffusion parameters from single Brownian trajectories.¹⁵ Our numerical approach allows one to efficiently estimate spatially varying dynamical parameters. It can provide useful results for even short trajectories, because it assumes a specific functional form for the spatial dependence of the dynamical parameters. One can readily characterize the statistical error and bias in these estimates through the analysis of ensembles of simulated data. Here, we provide a detailed description of our numerical maximum likelihood analysis (MLA) of Brownian trajectories. First, we provide a simple review of the principles of Brownian motion and maximum likelihood. Second, we demonstrate the method based on simulations and experimental data, including the analysis of strong interactions between pairs of paramagnetic particles. Finally, we discuss the practical implementation of the method, and provide the corresponding MATLAB code in the ESI.†

2 Theoretical background

In this section, we briefly review the solution to the Smoluchowski equation and the principle of maximum likelihood. We illustrate how to combine these two concepts with a simple example.

Let us consider the one-dimensional position x of a particle over time, and call F the applied force and D the diffusion coefficient. In general, F and D may depend on x. However, on sufficiently short time intervals Δt the particle samples a region where F and D are uniform. In this case, the solution of the Smoluchowski equation states that the displacements δx = x(t + Δt) − x(t) over Δt follow a Gaussian distribution

of mean μ and variance σ² given by

μ = v̄Δt, (2)

σ² = 2DΔt, (3)

where v̄ is the drift velocity:

v̄ = bF + ∇D. (4)

Here, b is the mobility and is related to D by the Stokes–Einstein relation:

D = bk_BT, (5)

with k_B the Boltzmann constant and T the absolute temperature.

The time evolution of the probability distribution for a constant force and diffusion coefficient is illustrated in Fig. 1a. From an experimental trajectory X = {x₁,…,x_N+1} with time interval Δt, F and D can be directly calculated from the mean and variance of the displacements, δx_i = x_i+1 − x_i, according to eqn (2) and (3).

Fig. 1 Brownian basics. (a) Probability density (normalized at each time point for visibility) of observing a room temperature Brownian particle with D = 10⁻¹³ m² s⁻¹ and F = 50 fN at position x at time t, starting from x(0) = 0, and simulated random trajectory (black line). Darker colors indicate higher probability. (b) Log-likelihood landscape in the (D,F) plane associated with the trajectory in (a). For better visibility, the log-likelihood is renormalized by the maximum value on the grid. Darker colors indicate higher log-likelihood values.

In contrast, the principle of MLA is to identify the set of dynamical parameters, α, that are most likely to describe an observed set of displacements, δ = {δx_i}_i=1…N. The likelihood that α describes the observed displacements δ can be expanded as

Here, p(α|δx_i) is the probability that α underlies the particle dynamics given a single observed displacement δx_i. Eqn (6) assumes the displacements are independent. This is true in the non-inertial regime, where the displacements are measured over a time Δt ≫ m/γ, with m the particle mass and γ its drag coefficient. The inertial relaxation time scales with the square of the particle diameter. For micrometric particles in water, m/γ ∼ 10⁻⁶ s.

Smoluchowski's theory provides the probability of observing a specific displacement given a set of parameters, p(δx_i|α). Using Bayes' theorem, we can determine p(α|δx_i):

Assuming the priors p(α) and p(δx_i) to be uniform on some reasonable intervals, the likelihood that a set of parameters α describes an observed set of displacements δ readswhere ω is a constant.

For numerical stability, it is more convenient to work with log-likelihood functions. For an observed trajectory X of corresponding set of displacements δ, we define the log-likelihood function of argument α as

It is important to note that L_δ(α) depends on δ, that is, on the trajectory considered.

The determination of α = (F,D) from a single Brownian trajectory is illustrated in Fig. 1. A simulated trajectory for a room-temperature Brownian particle with D = 10⁻¹³ m² s⁻¹ and F = 50 fN is super-imposed over the probability distribution as the black curve in Fig. 1a. The log-likelihood landscape associated to this trajectory in the α = (F,D) phase space is plotted in Fig. 1b. The log-likelihood function is maximized at the appropriate values of F and D, as visible in Fig. 1b and mathematically supported in the ESI.†

The utility and limitations of this approach are made apparent in the following sections.

3 Illustrative examples

3.1 Experimental trajectory of a single particle with constant F and D

We compare measurements of the force and diffusion coefficient of a single Brownian particle using two different methods: (1) direct calculation of the time-dependent mean and variance of the trajectory,⁹ and (2) MLA. Paramagnetic spheres (2.8 μm diameter) are suspended in water and sediment against a glass coverslip. A permanent magnet is positioned next to the sample to create a locally uniform gradient of the magnetic field B along the x-direction. The field gradient drives the paramagnetic particle with constant external force. Further experimental details are provided in the ESI.† A representative trajectory, x(t), of the particle is shown in Fig. 2a. Here, the frame-to-frame time interval is Δt = 2 ms, and the exposure time is τ_ex = 0.5 ms.

Fig. 2 Comparison of MLA and conventional analysis with an experimental trajectory (a) of a paramagnetic particle in a constant magnetic field gradient (inset). Mean (b) and variance (c) of the displacements δx at different lag times nΔt (blue dots), and corresponding linear fits (red lines). The slopes are equal to v̄ and 2D, respectively. (d) Frame-to-frame displacements δx_i as a function of position x_i. Red line indicates MLA fit to the mean (1.78 nm), and red dashed lines MLA fit to the standard deviation (±15.8 nm). (e) Log-likelihood landscape (renormalized) of the trajectory in (a). Red circle indicates results from statistical analysis.

From this trajectory, we calculate v̄ (Fig. 2b) and D (Fig. 2c) from (respectively) the mean and variance of the displacements at different lag times δx = x_i+n − x_i = x(t_i + nΔt) − x(t_i). We obtain the following estimations: F = 58 ± 3 fN and D = 6.16 ± 0.14 × 10⁻¹⁴ m² s⁻¹.

Alternatively, we can apply MLA to the distribution of the frame-to-frame displacements {δx_i = x_i+1 − x_i}, using parametrization α = (F,D). The results for the MLA fits of the mean and standard deviation are presented in Fig. 2d. The log-likelihood landscape for this trajectory in the (F,D) plane is represented in Fig. 2e. It shows a maximum at (F,D) = (59 ± 2 fN, 6.28 ± 0.15 × 10⁻¹⁴ m² s⁻¹), which corresponds to the values obtained from the conventional statistical analysis. This demonstrates that MLA is reliable in this simple case of constant F and D.

3.2 Simulated trajectories of pair interactions

We focus now on the more complex case of pair interactions. We assume that both the interparticle force, F, and the relative diffusion coefficient, D, depend on the center-to-center distance r (Fig. 3a-inset). We perform simulations corresponding to pairs of micron-sized particles interacting through capillary interaction with a rough contact-line. The force profile is given by

F(r) = −ϕ(x/a)^−β, (10)

with β = 5, following well-known theories,¹⁶ and a = 1 m a scaling factor (necessary for dimensional consistency). We choose a diffusion coefficient dependence corresponding to two particles in an isotropic fluid¹⁷which works well for a capillary interaction when the two fluids meeting at the interface have the same viscosity. Here, D₀ is the one-particle diffusion coefficient at infinite separation, and R₀ is the hydrodynamic radius of the particle.

Fig. 3 Testing MLA for spatially varying force and diffusion coefficients with simulated trajectories. (a) Sample trajectory (green line), and corresponding MLA fit (black line) from a simulation of colloidal spheres interacting via (inset) force F(r) and diffusion D(r) as described in eqn (10) and (11). (b) All simulated trajectories. Darker shades signify a higher density of curves. (c–f) Probability density function (PDF) of the fit parameters obtained from MLA. Red curves correspond to trajectories with Δt = 0.1 ms and τ_ex = 0 ms. Blue curves correspond to trajectories with Δt = 4 ms and τ_ex = 1 ms. (g and h) Fitted displacement (g) and force (h) profiles obtained from MLA for all trajectories. Black curves correspond to input profiles.

We perform N_s = 1000 simulations with timestep of Δt_s = 0.1 ms. These trajectories depend on 4 parameters α = {ϕ, β, D₀, R₀}. The input values are {1.2 × 10⁻⁴⁰ N, 5, 5.5 × 10⁻¹⁴ m² s⁻¹, 1 × 10⁻⁶ m}. A single simulated trajectory is shown in Fig. 3a, along with its corresponding MLA fit. It is important to note that due to the stochastic nature of Brownian motion, each trajectory is different. We plot a superposition of all trajectories in Fig. 3b to show this variability.

We estimate the underlying physical parameters for each trajectory using MLA. Details on the numerical optimization are given in the ESI.† As shown in Fig. 3c–f, MLA estimates (red histograms) agree very well with the input parameters (vertical black line). The spatial dependence of the particle displacements and the force are also accurately captured by MLA, as shown in Fig. 3g and h. The red bands are a super-position of all the fitted profiles, and the black lines show the input to the simulation. The estimated profiles match very well the input profiles, with most of the deviations not exceeding about 10% of the actual profiles. We discuss the effect of trajectory blurring, present in video microscopy experiments, on the accuracy of MLA estimates in Section 4.1.

These results show that MLA gives a very good estimation of the actual force profile and dynamic parameters for simulated trajectories.

3.3 Experimental trajectory of magnetic dipole–dipole pair interactions

We now investigate the reliability of MLA on experimental trajectories of isolated pairs of paramagnetic spheres in a magnetic field B. The magnetic field induces a magnetic dipole in both spheres, causing them to attract each other. Neglecting the mutual induced dipoles effect occurring at short distances, the interaction is well described by a power-law force (eqn (10)) with exponent¹⁸β = 4. Because the beads are heavy and sedimented, their hydrodynamic coupling should also include a contribution from the bottom surface of the observation chamber. For the diffusion coefficient of the separation, we use a three-parameter functional form:

D(r) = (D_r⁻¹(r, D₀, R₀) + (D_h⁻¹ − D₀⁻¹)/2)⁻¹ (12)

where D_r(r, D₀, R₀) is given by eqn (11), D₀ and R₀ are defined as before, and D_h is the diffusion coefficient of an isolated, sedimented bead at a small distance h from the bottom surface.¹⁹ As we will see later, this simple three-parameter form matches the exact numerical solution of the Stokes equations²⁰ well in the range of moderate interparticle distances considered in this paper.

We apply MLA to a pair trajectory using functional forms described in eqn (10) and (12), hence α = {ϕ, β, D₀, R₀, D_h}. The raw trajectory and displacement profile are shown in Fig. 4a and b. The fitted trajectory, diffusion profile and force profile are shown in Fig. 4a, c and d.

Fig. 4 Experimental test of MLA for spatially varying force and diffusion coefficient. (a) Trajectory of the separation between two paramagnetic beads in a magnetic field (inset). (b) Corresponding displacements as a function of separation for the part of the trajectory before contact. (c) Grey dots: diffusion coefficients obtained from the variance of the displacements (3 × 10⁵ time points). Red curve: diffusion profile obtained from MLA of the trajectory shown in (b), (1400 time points). Green curve: computational solution for two particles near a wall. (d) Force profiles obtained from MLA of the trajectory shown in (b) with some of the colored points removed. The color and spatial extent of the curve indicate the range of separations considered. For example, the orange curve is the force fit to the trajectory including the black, light yellow and yellow (r_i, δr_i) points, but excluding the dark orange, red and dark red points in (b).

We investigate the reliability of the fits in two ways. First, we compare the diffusion profile obtained from this single trajectory to an independent measurement. Second, we look at how the MLA fit is modified when a few points of the trajectory are artificially removed: a robust fit ought to be resilient to small changes in the trajectory.

3.3.1 Diffusion profile. To independently measure the diffusion profile, we tracked the displacements of the same particles in the absence of a magnetic field. Here, there is no long-range conservative interaction between the particles. We record the positions of the particles at over 3 × 10⁵ timepoints, which allows us to thoroughly sample the variance of the displacements at randomly-sampled separations from near contact to a few radii¹⁷. The diffusion profile obtained from the binning of this large dataset is presented in Fig. 4c as the grey dots. In addition, we calculated the height- and separation-dependent mutual diffusion coefficient with an exact numerical approach.^20,21 By fitting this calculation to the data, we determined the height above the wall to be 65 nm, and obtained the profile shown as the green curve. The experimental and theoretical diffusion profiles are in good agreement with the MLA estimate of the diffusion profile measured for the strongly-interacting particles in a magnetic field, shown in red. Notably, the MLA estimate required only 1400 timepoints, about 200 times less data than the direct measurement of the displacement variance. In the presence of the magnetic field, the trajectory length was limited by the strong attraction between the beads and there was not enough data to accurately sample the variance in each spatial bin.

3.3.2 Fit resilience to trajectory sampling. We assess the robustness of MLA by investigating how the estimated force profile is modified when a few points are removed from the trajectory.

The displacements δr_i as a function of separation r for the part of the trajectory before contact (Fig. 4a, t ≤ 2.8 s) are presented in Fig. 4b. The large displacements at small separations (colored points in Fig. 4b) may be expected to dominate the force profile. To test this hypothesis, we manually removed some of these (r_i, δr_i) points, and performed MLA of the truncated trajectories. The resulting force profiles are compared in Fig. 4d, where the full trajectory is dark red and the line colors approach yellow as the trajectory is more strongly truncated (the colors correspond to the last points included in Fig. 4b). Truncation of the data at small separations does change the estimated force profile, but the effects are very small (of the order of a few percents), even for the shortest trajectory (light yellow), where most of data points close to contact have been removed. This shows that MLA provides a robust fit, in the sense that small changes in a trajectory have only small effects on the results from the fit.

4 Practical considerations

In this section, we present some general guidelines regarding the practical implementation of MLA.

4.1 Data acquisition

Experimental trajectories of Brownian particles are typically acquired using video microscopy, which involves two key temporal parameters: the frame-to-frame time interval Δt, and the exposure time τ_ex. Finite exposure times are well known to cause systematic errors in the observed variance of displacements at short time intervals.^22,23 To mimic realistic experimental conditions, we undersample and blur the simulated trajectories of Section 3.2. We use a readily accessible time interval Δt = 4 ms and an exposure time τ_ex = 1 ms. (To emulate finite exposure times, we simply average the simulated positions over a 1 ms interval.) The results of MLA estimates of the parameters from these trajectories are shown as blue histograms in Fig. 3c–f. They agree reasonably well with the input parameters, with small but significant systematic errors. Similarly, blurring causes small but significant systematic errors in the displacement and force profiles, shown by the blue bands in Fig. 3g and h. The displacements tend to be more dispersed near contact, and the forces are systematically overestimated. To minimize these systematic errors due to finite camera exposure time, it is important that τ_ex be significantly smaller than Δt. More generally, typical particle displacements over the considered time interval, Δt, should be large compared to uncertainties in measuring the particle displacement.

4.2 Statistical error estimates

Statistical errors for a force profile obtained by MLA can be evaluated by simulating trajectories. Consider an experimental trajectory fitted using functional forms F_e(r) and D_e(r), where MLA has returned α₀ as the most likely estimate of the parameters. One can then perform a Brownian dynamics simulation to simulate a large number of trajectories with the same time interval, length and spatial domains as the real experiments, using F_e(r), D_e(r), and α₀. Each simulated trajectory can be analyzed to provide estimates of the parameters as well as the force and diffusion profiles. The spreads in these values and fits capture the statistical uncertainty on the MLA estimates. To illustrate, we plot in Fig. 5a the local density of force profile curves corresponding to the 1000 simulations described in Section 3.2. From this two-dimensional histogram, we can calculate the 5th and 95th percentiles at each separation r (black dashed lines), and hence recover a 90% confidence interval. The MATLAB code at the core of such analyses is included in the ESI.†

Fig. 5 Statistical and systematic errors. Density of force curves obtained from MLA of simulated trajectories (Δt = 0.1 ms, τ_ex = 0). Darker colors indicate higher density. Solid black curve indicates the force or diffusion profile. The two black dashed curves indicate the 90% confidence interval, inferred using correct functional forms, from the density of curves at each separation. (a) Histogram of the density of force profiles obtained from MLA of simulated trajectories using the correct functional forms for force and diffusion. (b) Force profiles from MLA of simulated trajectories using the approximate force functional form F*(r) described in Section 4.3. (c) Force profiles from MLA of simulated trajectories using the approximate diffusion functional form D*(r) described in Section 4.3.

4.3 Functional form choice and resilience

When the force and diffusion profiles vary with position, MLA requires functional forms for each. In some cases, these functional forms may be known exactly from the literature, but in other cases, they are not. In the latter case, one must use approximate functional forms F*(r) and D*(r). Reasonable ones can be chosen based on some basic knowledge of the underlying physics, such as continuity, monotony, and limiting behavior. To characterize the impact of using an incorrect but physically reasonable functional form, we analyzed trajectories from Section 3.2 using functional forms different from the ones used for the simulations, eqn (10) and (11). MLA estimates of the force profile using an incorrect force profile F*(r) = −ϕ exp(−r/β) are shown in Fig. 5b. Although they agree well with the input profile in the far-field, they significantly underestimate the force in the near-field. Similarly, estimates of the force profile using an incorrect functional form for the diffusion profile D*(r) = 2D₀(1 − (3/2)r₀/r), shown in Fig. 5c, lead to small but significant errors in the force estimate. However, in both cases we see that even with these wrong functional forms, the fitted force profiles remain essentially confined to the 90% confidence interval, except maybe at very small separations with the exponential force profile, where the fitted profiles are on average ∼20% below what they should be (Fig. 5b).

5 Conclusions

Maximum likelihood analysis is a powerful method to extract forces and diffusion coefficients from individual Brownian trajectories. It is fast and easy to implement. Our numerical approach to MLA is particularly useful when the amount of data is limited and the force and diffusion profiles vary over space. It can be very efficient, achieving small statistical uncertainties with a modest amount of data, by exploiting the smooth spatial dependence of the force and diffusion profiles. It can be very accurate, when the correct forms for the spatial dependence are employed and the typical particle displacements are large compared to any static or dynamic particle locating errors.

In this paper and our previous implementation of the method,¹⁵ we considered experimental trajectories of micron-sized colloidal particles. Fundamentally, this technique could also be employed with much smaller objects, such as fluorescently labeled proteins. Advanced microscopic techniques now allow for measurements of molecule positions with spacial resolution around 10 nm and time resolution better than 1 ms.²⁴ However, due to photobleaching and other experimental difficulties, the recorded trajectories are typically very short, hence making the study of molecular interactions very difficult. We believe that extension of MLA in the domain of single molecule tracking could not only be useful for measuring diffusion coefficients, as already established, but also help to understand long-range interactions between biological components, hence clarifying the molecular cell processes which work by an interplay of diffusion and interaction.

This work was supported by the National Science Foundation (CBET 12-36086). We would like to thank Jason Merrill, Robert Style, and Larry Wilen for helpful discussions.

References

1. N. A. Clark, J. A. Lunacek and G. B. Benedek, Am. J. Phys., 1970, 38, 575.
2. B. J. Berne and R. Pecora, Dynamic Light Scattering, Dover Publications, Inc., Mineola, New York, 2000.
3. T. G. Mason, K. Ganesan, J. H. van Zanten, D. Wirtz and S. C. Kuo, Phys. Rev. Lett., 1997, 79, 3282.
4. Y. Roichman, B. Sun, Y. Roichman, J. Amato-Grill and D. G. Grier, Phys. Rev. Lett., 2008, 100, 013602.
5. J. C. Crocker and D. G. Grier, Phys. Rev. Lett., 1994, 73, 352.
6. J. Perrin, Atoms, D. Van Nostrand company, New York, 1916.
7. D. C. Prieve and N. A. Frej, Langmuir, 1990, 6, 396.
8. S. Chandrasekhar, Rev. Mod. Phys., 1943, 15, 1.
9. S. K. Sainis, V. Germain and E. R. Dufresne, Phys. Rev. Lett., 2007, 99, 018303.
10. J. W. Merrill, S. K. Sainis, J. Bławzdziewicz and E. R. Dufresne, Soft Matter, 2010, 6, 2187.
11. D. J. Evans, A. D. Hollingsworth and D. G. Grier, Phys. Rev. E, 2016, 93, 042612.
12. A. J. Berglund, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 82, 011917.
13. X. Michalet and A. J. Berglund, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 061916.
14. C. L. Vestergaard, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2016, 94, 022401.
15. R. Sarfati and E. R. Dufresne, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2016, 94, 012604.
16. D. Stamou, C. Duschl and D. Johannsmann, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2000, 62, 5263.
17. P. L. Biancaniello and J. C. Crocker, Rev. Sci. Instrum., 2006, 77, 113702.
18. D. Du, F. Toffoletto and S. L. Biswal, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 89, 043306.
19. E. R. Dufresne, T. M. Squires, M. P. Brenner and D. G. Grier, Phys. Rev. Lett., 2000, 85, 3317.
20. S. Bhattacharya, J. Bławzdziewicz and E. Wajnryb, Physica A, 2005, 356, 294.
21. S. Bhattacharya, J. Bławzdziewicz and E. Wajnryb, J. Fluid Mech., 2005, 541, 263.
22. W. P. Wong and K. Halvorsen, Opt. Express, 2006, 14, 12517.
23. T. Savin and P. S. Doyle, Biophys. J., 2005, 88, 623.
24. A. Kusumi, T. A. Tsunoyama, K. M. Hirosawa, R. S. Kasai and T. K. Fujiwara, Nat. Chem. Biol., 2014, 10, 524–532.

---

Footnote

† Electronic supplementary information (ESI) available: Materials & methods, and mathematical derivations. See DOI: 10.1039/c7sm00174f

---