Learning Force Fields from Stochastic Trajectories

When monitoring the dynamics of stochastic systems, such as interacting particles agitated by thermal noise, disentangling deterministic forces from Brownian motion is challenging. Indeed, we show that there is an information-theoretic bound, the capacity of the system when viewed as a communication channel, that limits the rate at which information about the force field can be extracted from a Brownian trajectory. This capacity provides an upper bound to the system’s entropy production rate and quantifies the rate at which the trajectory becomes distinguishable from pure Brownian motion. We propose a practical and principled method, stochastic force inference, that uses this information to approximate force fields and spatially variable diffusion coefficients. It is data efficient, including in high dimensions, robust to experimental noise, and provides a self-consistent estimate of the inference error. In addition to forces, this technique readily permits the evaluation of out-of-equilibrium currents and the corresponding entropy production with a limited amount of data.

Popular Summary

Stochastic differential equations are commonly used to model the dynamics of systems subject to quickly varying, apparently random forces. The most popular of these equations is Brownian dynamics, where the observable variables are subject to deterministic forces and random white noise. It represents well the dynamics of submicrometer particles, which abound in biological and soft-matter systems. Brownian dynamics is also used to model the time evolution of more complex systems such as climate, cell motion, and financial markets. Here, we address the inverse problem of Brownian dynamics: how to reconstruct the force field and the noise from an analysis of experimental observations.

We first identify a key bottleneck in this inference problem: By introducing a new mapping of Brownian dynamics onto communication theory, we find that information about the force field can be obtained only at a finite rate, quantified by a channel capacity. We relate this capacity to the entropy production rate, providing a novel link between stochastic thermodynamics and information theory.

We then propose a practical method, stochastic force inference, that uses this information to reconstruct the force and diffusion fields. Briefly, it consists of a smooth regression of these fields with known functions, where both underfitting and overfitting are controlled. This method can readily infer entropy production, currents, interactions, and state-dependent diffusion fields in out-of-equilibrium many-body systems in a way that is robust to experimental noise.

Our approach provides a powerful and data-efficient solution to reconstructing the force field underlying Brownian dynamics, which could be used on a broad class of stochastic systems where inferring a dynamical model from limited noisy data is of interest.

Figure 1

Typical trajectories of example Brownian systems studied in this article. (a) Pure Brownian motion in 2D without forces. (b) A drifted Brownian motion trajectory. (c) The stochastic Lorenz process (see Fig. 5). (d) Time series of a 6D out-of-equilibrium Ornstein-Uhlenbeck process (see Fig. 4). (e) The same trajectories as in (d), with additional time-uncorrelated measurement noise. (f) Self-propelled active Brownian particles with soft repulsion and harmonic confinement (see Fig. 7). (g) Simulated single-molecule trajectories in a complex environment with space-dependent diffusion (see Fig. 9).

Figure 2

The dynamics of an overdamped system can be seen as a noisy data-transmission channel encoding information about the force field, with a rate bounded by the channel capacity $C$ as defined in Eq. (2). Note that this definition does not include the information loss stemming from the measurement device. This analogy is further discussed in the Appendix pp1.

Figure 3

Stochastic force inference for a 2D Ornstein-Uhlenbeck process with force field $F_{\mu }(\mathbf{x})=-{\mathrm{\Omega }}_{\mu \nu }{x}_{\nu }$ and isotropic diffusion. (a) An example trajectory. The inferred force field for this trajectory using SFI with functions $b=\{1,x_{\mu }\}$ (blue arrows) is compared to the exact force field (black arrows). Inset: The inferred force components along the trajectory versus the exact force components with normalized mean-squared error (MSE). (b) The average of the relative error $[({\widehat{F}}_{\mu \alpha }-F_{\mu \alpha }^{\tau })D_{\mu \nu }^{-1}({\widehat{F}}_{\nu \alpha }-F_{\nu \alpha }^{\tau })]/[{\widehat{F}}_{\mu \alpha }{D}_{\mu \nu }^{-1}{\widehat{F}}_{\nu \alpha }]$ on the inferred projection coefficients ${\widehat{F}}_{\mu \alpha }$ and its self-consistent estimate $N_b/2{\widehat{I}}_b$ both converge to $N_b/2I_b$, as expected from theory (see Appendix pp3). Here, $F_{\mu \alpha }^{\tau }=\int F_{\mu }\mathbf{(}\mathbf{x}(t)\mathbf{)}{\widehat{c}}_{\alpha }\mathbf{(}\mathbf{x}(t)\mathbf{)}(dt/\tau )$ is the projection of the exact force on the empirical projectors.

Figure 4

(a) Time series of a 6D out-of-equilibrium Ornstein-Uhlenbeck process, with anisotropic harmonic confinement and diffusion tensor, and circulation. The force field is $F_{\mu }(\mathbf{x})=-{\mathrm{\Omega }}_{\mu \nu }{x}_{\nu }$. The matrix $\mathrm{\Omega }$ and the diffusion matrix are chosen from a random ensemble. The antisymmetric part of $D^{-1}\mathrm{\Omega }$ has rank two, thus inducing circulation in a randomly chosen plane. (b) The same trajectories as in (d) with additional time-uncorrelated measurement noise. (c) SFI for the trajectory in (a) allows precise identification of the plane of circulation and reconstruction of the force along the trajectory. (d) SFI applied to the trajectory in (b) with measurement noise. It can still detect forces accurately. (e) Convergence of the angular error for cycle detection with increasing trajectory length for the process shown in (d) and (e). (f) Inferred entropy production rate for this process with and without measurement noise (we subtract here the systematic bias $2N_b/\tau$). The shadowed area indicates the self-consistent confidence interval for the inferred entropy production. The dotted line shows the exact value of the entropy produced; for the noisy process, SFI underestimates this value due to blurring of the currents. (g) Entropy production captured when observing a $d$-dimensional projection of the trajectory averaged over direction of observation for long trajectories. In (e)–(g), the error bars indicate the standard deviation over an ensemble of 32 trajectories. Parameters of the simulations are presented in Appendix pp8.

Figure 5

Stochastic force inference with nonlinear force fields. (a) Trajectory of an out-of-equilibrium process with harmonic trapping and circulation, and a Gaussian repulsive obstacle in the center. The force field is given by $F_{\mu }(\mathbf{x})=-{\mathrm{\Omega }}_{\mu \nu }{x}_{\nu }+\alpha e^{-x^2/2{\sigma }^2}{x}_{\mu }$ where $\mathrm{\Omega }$ has both a symmetric and antisymmetric part. (b) Trajectory of the stochastic Lorenz process, a 3D process with a chaotic attractor. The force field is $F_x=s(y-x)$, $F_y=rx-y-zx$, $F_z=xy-bz$, where we choose $r=10$, $s=3$, and $b=1$. (c)–(h) SFI for these two trajectories, respectively, with polynomials of order $n=1$, 3, 5 and $n=1$, 2, 3: inferred force versus exact force (left) and bootstrapped trajectory using the inferred force field (right). (i),(j) Capacity (top) and entropy production (bottom) of each process projected on different bases for an asymptotically long trajectory as a function of the number of degrees of freedom $N_b$ in the basis. These bases are polynomial and Fourier functions with order $n=0,\dots ,7$, and a coarse-grained approximation with a variable number of grid cells $n=2,\dots ,7$ in each dimension. Parameters and details of the simulations are presented in Appendix pp8.

Figure 6

Influence of the size of the basis on the precision of SFI. (a),(b) SFI error as a function of the number of fit parameters, respectively, for the models presented in Figs. 5 and 5 with a Fourier basis, and for different numbers of time steps in the trajectory. Specifically, the $y$ axis is the mean-squared relative error on the inferred force along the trajectory, $⟨({\widehat{F}}_{\mu }-F_{\mu })D_{\mu \nu }^{-1}({\widehat{F}}_{\nu }-F_{\nu })⟩/⟨{\widehat{F}}_{\mu }{D}_{\mu \nu }^{-1}{\widehat{F}}_{\nu }⟩$. The crossover from under- to overfitting is apparent and takes place at larger $N_b$ and lower error with longer trajectories. The star symbols indicate the optimal basis size predicted by our self-consistent criterion of maximizing ${\widehat{I}}_b-\delta {\widehat{I}}_b$. (c),(d) The squared error as a function of the amount of information $C\tau$ in a trajectory of duration $\tau$ for the optimal basis, averaged over $n=3$ trajectories. For the Lorenz process with a polynomial basis [(d) orange squares], the convergence is fast as the basis is adapted to the exact force field, and the saturation of the error to a lower plateau is due to the finite time step (see Appendix pp6).

Figure 7

Stochastic force inference for harmonically trapped active Brownian particles with soft repulsive interactions $F(r)=1/(1+r^2)$ between particles at distance $r$. (a) Snapshot of a configuration for 25 active particles. The black dots indicate the direction of self-propulsion. We perform SFI on a trajectory of only 25 frames blurred to mimic measurement noise. Background shows the trajectory of one particle and force on each particle, inferred (blue arrows) and exact (black arrows). The fitting basis for SFI consists of a combination of harmonic trapping, constant-velocity self-propulsion, and radial interactions between particles with the form $r^k{e}^{-r/r_0}$ with $k=0,\dots ,5$ and $r_0$ a typical nearest-neighbor distance between particles. (b) Inferred versus exact components of the force on all particles along the trajectory. (c) Inferred radial force between interacting particles, compared to the exact force.

Figure 8

Stochastic inference of inhomogeneous diffusion and forces. (a) A trajectory of a 1D ratchet model with $F(x)=F_0\cos (2\pi x)$ and $D(x)=1+a\cos (2\pi x)$, with periodic boundary conditions. (b),(c) For the trajectory presented in (a), inferred and exact diffusion coefficient [using Eq. (12)] and force field [using Eq. (15)] as a function of the position. We use a first-order Fourier basis to infer both force and diffusion. (d) Analysis of the convergence of the diffusion (blue) and force (orange) estimators as a function of the trajectory duration for the process presented in (a). The dotted and dashed black lines are the self-consistent estimates for the squared error, respectively, for the diffusion and the force. The plateau for the diffusion inference is due to the finite time step. (e) A trajectory of a minimal 2D model, an isotropic harmonic trap at equilibrium $F_{\mu }(\mathbf{x})=-D_{\mu \nu }(\mathbf{x})x_{\nu }$ in a constant gradient of isotropic diffusion $D_{\mu \nu }(\mathbf{x})=(1+a_{\rho }{x}_{\rho }){\delta }_{\mu \nu }$. (f),(g) Inferred versus exact diffusion coefficient [using Eq. (12)] and force components [using Eq. (15)] along trajectory (a). A linear polynomial basis is used to fit the diffusion coefficient and a quadratic basis to fit $F_{\mu }$. (d) Convergence of the diffusion projection estimator (normalized by the average diffusion tensor) to its exact value for the process shown in (a). Circles, using Eq. (12); diamonds, using Eq. (12) in the presence of time-uncorrelated measurement noise; triangles, using the bias-corrected local estimator. Error bars represent the standard deviation over 64 samples. Details and parameters in Appendix pp8.

Figure 9

Quantitative comparison of SFI with other methods on a simulated system mimicking 2D single-molecule trajectories in a complex cellular environment with multiple potential wells, out-of-equilibrium circulation, and space-dependent isotropic diffusion. (a) The diffusion field (blue gradient) and drift field (white arrows, scaled as $|\mathrm{\Phi }{|}^{1/2}$ for better legibility). (b) The steady-state probability distribution function (PDF) of the process. The blue traces show two representative trajectories with $n=100$ time steps. The red traces show trajectories blurred by moderate Gaussian measurement error (with amplitude shown as a red kernel). (c)–(f) Comparison of the performance of SFI with adaptive Fourier basis (green circles) and two widely used inference methods: InferenceMAP [25], a Bayesian method for single-molecule inference (blue triangles), and grid-based binning with maximum-likelihood estimation [23, 33] [Eq. (7)] and an adaptive mesh size (orange squares). We evaluate the performance of these methods on the approximation of the drift field [(c),(e)] and diffusion field [(d),(f)] as a function of the number $N$ of single-molecule trajectories [similar to those in (b)] used, with ideal data [(c),(d)] and in the presence of measurement noise [(e),(f)]. The performance is evaluated as the average mean-squared error on the reconstructed field along trajectories. SFI outperforms both other methods in all cases; for noisy data, SFI is the only one that provides an unbiased estimation of the drift. Details and parameters in Appendix pp8-s9.

Figure 10

Comparison of the performance of the three methods on a set of $N=128$ noisy trajectories. (a) The trajectories in the raw dataset exploited by each method. (b)–(d) The inferred drift field (thick light blue arrows) and the exact one used to generate the data (thin red arrows). The lower left insets show scatter plots of the inferred versus exact drift components, with an indication of the normalized mean-squared error (MSE). Top right inset of (b): Screen capture of the InferenceMAP software used for this drift field. While methods (b) and (c) capture the qualitative shape of the drift field, they appear biased, and consistently overestimate the drift for noisy data. In contrast, our method shows quantitative agreement with the exact drift field.