Inferring potential landscapes from noisy trajectories of particles within an optical feedback trap

Summary While particle trajectories encode information on their governing potentials, potentials can be challenging to robustly extract from trajectories. Measurement errors may corrupt a particle’s position, and sparse sampling of the potential limits data in higher energy regions such as barriers. We develop a Bayesian method to infer potentials from trajectories corrupted by Markovian measurement noise without assuming prior functional form on the potentials. As an alternative to Gaussian process priors over potentials, we introduce structured kernel interpolation to the Natural Sciences which allows us to extend our analysis to large datasets. Structured-Kernel-Interpolation Priors for Potential Energy Reconstruction (SKIPPER) is validated on 1D and 2D experimental trajectories for particles in a feedback trap.

Inferring naturally occurring energy landscapes or verifying artificially created potentials demands a method free of a priori assumptions on the potential's shape. This requirement rules out many commonly used methods devised for harmonic systems (Neuman and Block, 2004;Berg-Sørensen and Flyvbjerg, 2004;Jones et al., 2015;Gieseler et al., 2021) or alternative, otherwise-limited, methods to deduce potentials from data (Reif, 2009;Tü rkcan et al., 2012;García et al., 2018;Wang et al., 2019;Frishman and Ronceray, 2020;Yang et al., 2021;Stilgoe et al., 2021). For example, some methods (Reif, 2009;Tü rkcan et al., 2012) necessarily rely on binned data, relating potential energies to Boltzmann weights or average apparent force, thereby limiting the frequency of data in each bin and requiring that equilibrium be reached before data acquisition. Other methods assume stitched locally harmonic forms (García et al., 2018). Still others use neural networks (Wang et al., 2019) to deduce potentials; the uncertainty originating from measurement error and data sparsity is then not easily propagated to local uncertainty estimates over the inferred potential.
In previous work (Bryan IV et al., 2020), we introduced a method starting from noiseless one-dimensional time series data to infer effective potential landscapes without binning, or assuming a potential form, even in the case of limited sampling, while admitting full posterior inference (and thus error bars or, equivalently, credible intervals) over any candidate potentials arising from sparse data. Our method was, however, fundamentally limited to one dimension (because of the poor scaling of the computation with respect to the dataset size). As such, the method could not distinguish the effects of inherent stochasticity in the dynamics due to temperature, say, and noise introduced from the measurement apparatus. It would therefore be helpful to develop a method capable of discriminating between both noise sources and learning features of both just as the perennial hidden Markov model achieves for discrete state spaces (Rabiner and Juang, 1986;Sgouralis and Pressé , 2017).
Here, we introduce a method to infer potentials from noisy, multidimensional data often sampled in a limited fashion. We take advantage of tools from Bayesian nonparametrics to place priors on potentials assuming no prior functional forms for these potential. To do so, we introduce structured-kernel-interpolation Gaussian processes (Wilson and Nickisch, 2015) (SKI-GP) to the Natural Sciences in order to circumvent the otherwise-prohibitive computational scaling of widely used Gaussian processes. In our application, we utilize the power of SKI-GP to create Structured-Kernel-Interpolation Priors for Potential Energy Reconstruction (SKIPPER). SKIPPER can be inferred from trajectories while meeting all the following criteria simultaneously: 1) no reliance on binning or pre-processing, 2) no assumed analytic potential form, 3) inferences drawn from posteriors, allowing for spatially nonuniform uncertainties to be informed by local density of available data in specific regions of the potential (e.g., fewer data points around barriers), 4) treatment of multidimensional trajectories, 5) rigorous incorporation of measurement noise through likelihoods, and 6) compatible with lightly sampled trajectories. No other existing method meets all six criteria simultaneously.
This work serves the dual role of a method for inferring potentials with SKIPPER, as well as a demonstration of how to incorporate SKI-GPs into an inference framework. In particular, we demonstrate how to construct an effective kernel matrix for an unknown variable (the potential landscape) and how to directly sample the variable from a conditional posterior.

RESULTS
We benchmark SKIPPER on experimental data on a double-well potential and show that we can accurately infer the shape of the potential. We then show that SKI-GP allows us to explore 2D time series data (previously infeasible due to large amounts of data using a naive GP). We finally apply SKIPPER to trajectories in a high-barrier landscape where traces are too short to reach equilibrium. A demonstration on data from a simple harmonic well, a complicated 2D potential, and robustness tests over parameters of interest can be found in the SI (Bryan IV et al., 2022).
For testing the accuracy and effectiveness of SKIPPER, we simultaneously collected two measurements of each trajectory, one using a detector with low measurement noise and one using a second detector with higher measurement noise as shown in section. We refer to the low-noise trajectory as the ''ground truth'' trajectory, although it itself is subject to a small amount of measurement noise. For each experiment, we impose a potential on the particle using our feedback trap. We refer to this applied potential as the ''ground truth'' potential, although it may differ from the actual potential the particle experiences due to errors in the feedback trap setup, as well as experimental limitations such as drift.

Demonstration on simulated data
We started by analyzing simulated trajectories from a simple harmonic well. Results are shown in Figure 1. Each column shows the inferred potential (top row) and inferred trajectory (bottom row) for each dataset analyzed. We provide uncertainties and ground truth estimates for both the potential and trajectory. Additionally, for sake of comparison, we also show the potential estimated using the Boltzmann method (Bryan IV et al., 2020;Reif, 2009) (discussed in section). We highlight that the Boltzmann method does not provide trajectory estimates. By contrast, SKIPPER infers those positions obscured by noise. Figure 1 shows that the ground truth potential and trajectory fall within our error bars (credible interval) for all datasets up to ND = 7. At ND = 7, the inferred potential develops bumps where it is unable to infer the trajectory accurately, resulting in a 2 nm shift of the potential well minimum. On the other hand, the Boltzmann method (see earlier discussion) does well in the low-noise case (Figure 1, left column) but fails as the noise introduced grows (Figure 1, right column). This is expected, as measurement noise broadens the iScience Article position histogram and thereby the potential. Using SKIPPER, the estimate of the potential drops at the edges of the spatial region sampled by the particle, since the high-potential edges are rarely visited by the particle. We thus have insufficient information through the likelihood to inform those regions, and the inferred potential reverts back to the prior (set at 0, as described earlier).
In addition to inferring estimates for the potential energy landscape and trajectory, we also estimate the magnitude of the measurement noise for each experiment. Figure 2 shows the inferred measurement noise magnitudes obtained from SKIPPER for each single-well experiment, with the mean of our PDFs within 5% of the ''ground truth'' measurement, for each dataset analyzed.
Next, to show the full utility of SKIPPER, we analyzed simulated data from a complicated 2D potential with wells in the shape of a winky face. The winky face was constructed from a grid of Gaussian wells with potential energy corresponding to the pixels of a winky face. The winky face is 13 3 13 pixels, where each pixel is 10 nm wide. The eyes and mouth are about 2 kT less than the rest of the face. The white space around the face has potential energy above 20 kT, to ensure that the particle stays inside the face. We simulated a 50,000 data point trace of a particle moving along the winky face. Figure 3 shows results on this potential. Note that potential energies greater than 5 kT and potential energies in areas where the particle is not encountered are displayed as white background. In the top row of Figure 3, when measurement noise is small compared to the size of the face (1 vs 130 nm), both SKIPPER and the Boltzmann method are able to reconstruct a potential energy landscape that is recognizably a winky face. However, with moderate noise (s = 5 nm), SKIPPER is able to pick up the fine detail in the face that is missed by the Boltzmann method (Figure 3 middle row). In the case of high noise (s = 10 nm), SKIPPER's advantage over the Boltzmann method is more dramatic, as SKIPPER's reconstruction is still recognizably a winky face, whereas the Boltzmann method's reconstruction not only fails to resemble a face but also overestimates the width of the face (Figure 3 bottom row).

Robustness tests on simulated data
In order to probe SKIPPER's robusteness with respect to parameters of interest, we demonstrate SKIPPER on data simulated with a single-well potential analyzed under different circumstances.
First, we tested SKIPPER's robustness with respect to the number of data points on simulated data. Supplementary figure SI-1 shows the results on simulated data. Trivially, when the trajectory is so short that the particle does not travel across the well (N = 50 as in Fig. SI-1 left panel), we cannot infer the shape , we infer the general shape of the potential, but the inference is highly impacted by small stochastic anomalies. See for example the far right of the potential, where a few higher-and lower-than-expected thermal kicks at the right side of the well caused SKIPPER to infer a second well. For reasonable inference, a few thousand points suffice (N = 5000 in Fig We then tested SKIPPER's robustness with respect to measurement noise in order to probe the regime in which SKIPPER can be used. Supplementary figure SI-2 shows the results. SKIPPER was robust to measurement noise for the four values of measurement noise variance chosen, but could not reproduce the potential energy landscape when the magnitude of the noise was of the same order as the maximum range of the particle. Even so, when the magnitude of the noise was up to 10% the maximum range of the particle (s = 10 nm), we nonetheless inferred the potential energy landscape accurately.

Double well
We analyzed data from a real particle in a double-well potential. Results are shown in Figure 4. Each column shows the inferred potential (top row) and inferred trajectory (bottom row) for each dataset analyzed. We provide uncertainties and ground truth estimates for both the potential and trajectory. Additionally, for sake of comparison, we also show the potential estimated using the Boltzmann method (Bryan IV et al., 2020;Reif, 2009). Briefly, the Boltzmann method, as outlined in the SI (Bryan IV et al., 2022), assumes that the bead localizations are sampled from the Boltzmann distribution and derives an estimate of the potential in discretized bins of space from the log of the relative frequency the bead is seen in each bin. We highlight that the Boltzmann method does not provide trajectory estimates. By contrast, SKIPPER infers those positions obscured by noise. Figure 4 shows that the ground truth potential and trajectory fall within the estimated range even when the measurement noise is so large that the particle is occasionally seen in the wrong well ( Figure 4, top right panel). Both SKIPPER and the Boltzmann method slightly overestimate the potential of the left well at the lowest noise level, because the (short) trajectory spends too much time in the right well, leaving the left well undersampled.

2D single well
Next, we analyzed data from a real particle in a 2D harmonic potential. Results are shown in Figure 5. For clarity, we do not show uncertainties, trajectories, or Boltzmann-method estimates for the 2D plot, but we do show them for a 1D potential slice. Despite the added complexity in inferring the potential in full 2D at once, our estimates fall within uncertainty in regions where data are appreciably sampled, even at high measurement noise.

Trajectories with limited sampling
One advantage of SKIPPER is that it does not rely on equilibrium assumptions. As such, we can analyze trajectories initiating from far from equillibrium conditions and with limited sampling such that the particle does not reach equillibrium in the duration of the sampling. To demonstrate this, we created real datasets where the iScience Article particle starts at the top of the potential well and ''rolls off'' to either side. The trajectories are short (5 ms), so that the particle does not reach equilibrium during the time trace. By including the likelihoods from 100 such trajectories into our posterior, we gain information on either side of the well and can recreate the full potential, even though each individual trajectory is initiated from the top and samples only one well.
In the first four panels of Figure 6, we illustrate 4 of 100 trajectories used to reconstruct the potential. We note that all trajectories start from the top of the barrier (defined as x = 0), and none of the trajectories fully sample both wells.
As all trajectories start from the top of the well, we collect a disproportionally large number of data at the top of the well. Indeed, for our 10 kT barrier with roughly 50,000 samples, we would not expect Here, we show results from inference on a 2D potential with wells in the shape of a winky face. The left column shows the ground truth potential energy landscape. The middle column shows inference using SKIPPER. The right column shows inference using the Boltzmann histogram method. The top row analyzes data simulated with 1 nm measurement noise. The middle row analyzes data simulated with 5 nm measurement noise. The bottom row analyzes data simulated with 10 nm measurement noise. Note that in our display pixels with potential energy greater than 5 kT are white. iScience Article to see any samples at the top of barrier in an experiment taken after the bead had reached equilibrium. Despite this extreme oversampling, SKIPPER is able to infer the height of the barrier to within 15% accuracy (SKIPPER predicts an 8.6 kT barrier; the barrier of the design potential is 10 kT).
Our error bar at the top of the barrier in Figure 6 is artificially low because every trajectory initiates from the top.

DISCUSSION
In conclusion, inferring potential landscapes is a key step toward providing a reduced dimensional description of complex systems (Wang et al., 2019;Chmiela et al., 2017;Espanol and Zuniga, 2011;Izvekov and Voth, 2005;Manzhos et al., 2015). Here, we go beyond existing methods by providing a means of obtaining potentials, among multiple other quantities, from time series data corrupted by measurement noise. We do so by efficiently learning the potential from the rawest form of data, point-by-point. That is, we achieve this without data pre-processing (e.g., binning), assuming an analytical potential form, nor requiring equilibrium conditions. As SKIPPER is Bayesian, it allows for direct error propagation to the final estimate of the inferred potential shape. In other words, SKIPPER differs from others assuming analytic potential forms (Pé rez-García et al., 2021) or projection onto basis functions (Frishman and Ronceray, 2020), as well as methods relying on neural nets (Manzhos et al., 2015) that cannot currently propagate experimental uncertainty or provide error bars reflecting the amount of data informing the potential at a particular location. Importantly, unlike the Boltzmann method (Reif, 2009), SKIPPER does not invoke any equilibrium assumption and can consider trajectories initiated from positions not sampled from an equilibrium distribution. This feature is especially relevant in studying landscapes with rarely sampled regions of space and, in particular, far from equilibrium.  In this work, we assume a Gaussian noise model, but we can, in principle, utilize different noise models by substituting Equation 4 for the desired model. As SKI-GP is general and the measurement noise model can be tuned, moving forward we could apply modified SKIPPER algorithms to map potential landscapes from force spectroscopy (Gupta et al., 2011) or even single molecule fluorescence energy transfer (Kilic et al., 2021;Sgouralis et al., 2018), with applications to inferring protein conformational dynamics or binding kinetics (Schuler and Eaton, 2008;Chung and Eaton, 2018;Sturzenegger et al., 2018;Pressé et al., 2013Pressé et al., , 2014. In inferring smooth potentials, we would move beyond the need to require discrete states inherent to traditional analyses paradigms such as hidden Markov models (Rabiner and Juang, 1986;Sgouralis and Pressé , 2017).
Beyond inferring potentials, we believe that SKI-GPs may become a powerful tool within the Natural Sciences especially in cases, such as this one, where we deal with noisy and abundant data while learning continuous functions from data. The SKI-GP prior shifts the computational burden from inferring the potential at every data point, to inferring the potential at select inducing points allowing the computation time to scale cubically with the number of inducing points rather than cubically with the number of data points. As the number of inducing points required to create a detailed map of the potential is often far less than the number of data points, this allows for significant computational cost reduction. Furthermore, as the inducing points are static, they allow for pre-computation of the kernel matrix despite the fact that locations of the data points change at each iteration of the Gibbs sampler. As demonstrated in this work, SKI-GPs can be used to accurately and tractably map fields of variables. Such fields are ubiquitous in nature, including temperature maps (Aigouy et al., 2005), optical absorption coefficient maps (Yuan and Jiang, 2006), and diffusion coefficient maps (Taylor and Bu shell, 1985;Weistuch and Pressé , 2017). We believe that this work can serve as a demonstration for how to incorporate SKI-GPs into inference frameworks that can be extended beyond potential learning.

Limitations of the study
The model formulation requires that potentials are time independent and measurement noise is Gaussian distributed. The time scales and drag coefficient are assumed to be large enough to approximate motion with overdamped Langevin dynamics.

Figure 6. Demonstration of data from experimental trajectories with limited sampling
We reconstruct a potential by analyzing many short (500 data points) trajectories with limited sampling. The left four panels show four of the 100 small data segments used to reconstruct the potential. Each trajectory starts at the top of the potential and rolls off to either side. The far right shows the inferred potential plotted with uncertainty overlaid on the ground truth potential and the potential inferred using the Boltzmann method. For comparison, the inferred and ground truth energy landscapes were shifted so that the lowest point is set to 0 kT. Geman, S., and Geman, D. (1984).

Materials availability
This study did not generate new unique materials.

Data and code availability
Data and code are available via Zenodo.
Any additional information required to reanalyze the data reported in this work paper is available from the Lead contact upon request.

Experimental apparatus
The experiment is done on a modified version of the feedback optical tweezer used in various stochastic thermodynamics experiments involving virtual potentials (Kumar and Bechhoefer, 2018a, 2020). The schematics of the setup are provided in Fig. SI-3. We used a continuous-wave diode-pumped laser (HÜ BNER Photonics, Cobolt Samba, 1.5 W, 532 nm) and a custom-built microscope to construct the optical trap on a vibration-isolation table (Melles Griot). We split the laser into a trapping-beam and a detectionbeam using a 90:10 beam splitter. The trapping beam then passes through a pair of acousto-optic deflectors (DTSXY-250-532, AA Opto Electronic), which allows us to deflect the beam in the orthogonal XY plane. We can steer the angle of the beam and also its intensity using analog voltage-controlled oscillators (DFRA10Y-B-0-60.90, AA Opto Electronic). We increase the beam-diameter using a two-lens system in telescopic configuration to overfill the back aperture of the trapping microscope objective (MO1 in Figure S3). Then a 4f relay system images the steering point of the AOD on the back aperture of MO2. A water-immersion, high-numerical-aperture objective (MO2, Olympus 60X, UPlanSApo, NA = 1.2) is used for trapping 1.5 m diameter spherical silica beads (Bangs Laboratories) in aqueous solution (SC).
The detection beam is passed through a half-wave plate (/2), so that its polarization is orthogonal to the trapping beam, to avoid unwanted interference. It is focused using a low-numerical-aperture microscope objective (MO1, 40X, NA = 0.4) antiparallel to the trapping objective MO2. The loosely focused detection beam (compared to the trapping beam) has a larger focal spot and offers a high linear range for the position detection. We can also adjust the detection plane using a 4f relay lens system. The forward scattered detection beam from the trapped bead is collected by the trapping objective MO2 and transmitted through the polarizing beam splitter (PBS). The PBS reflects the trapping beam and thus separates the detection beam. We also place another linear polarizer (P) after the PBS to minimize the amount of leaked trapping laser. The detection beam is then separated using a beam splitter on a pair of quadrant iScience Article photodiodes (QPD, First Sensor, QP50-6-18u-SD2). We control the intensity of the detection beam on QPD1 by placing a neutral density filter (ND) of desired strength and thus control the signal-to-noise ratio (SNR) on QPD1.
A red LED (660nm, Thorlabs, M660L4) illuminates the trapped particle for imaging onto a camera. The light from the LED enters the detection object through a long-pass filter (LPF, cutoff wavelength 585 nm, Edmond Optics) which reflects the detection beam. The collected illumination light is separated by a short-pass filter (SPF, cutoff wavelength 600 nm, Edmond Optics) and reflected onto a camera (FLIR BFS-U3-04S2MÀCS).
Using a LabView program, we digitize the analog signal from the pair of QPDs by a field programmable gate array (FPGA, National Instruments, NI PCIe-7857). The voltage signals are then calibrated as position signals using a QPD-AOD-camera method (Kumar and Bechhoefer, 2018b). The FPGA runs the control loop and uses a programmed feedback rule to send the appropriate control signal to the AODs.

Measurement noise
Here we discuss the estimation of measurement noise for our experimental setup. Briefly, we obtain an estimate of the noise by recording two simultaneous trajectories of a particle and finding the difference between the trajectories. As both measured trajectories should be equal in the limit of zero noise, we can interpret any deviation between the measured trajectories as coming from noise.
We start by applying a triangle wave voltage of frequency 1.6 Hz to the AOD, to move a bead linearly along the x axis. We record the beads motion at 100 kHz using two simultaneous detectors giving output x 1 ðtÞ and x 2 ðtÞ, respectively. We demonstrate this idea in Figure S4A. Figure S4B shows that the difference x 1 À x 2 follows a Gaussian distribution. The mean of the Gaussian distribution shown in Fig. SI-4B, has a nonzero value of 0:28G0:02 nm, which arises from nonlinear calibration errors in the two detectors. In Figure S4C, we verify that the autocorrelation is flat, indicating that each position measurement has independent measurement noise.
To calculate the variance of the measurement noise, we denote the SD of the detectors s 1 and s 2 , respectively. For Gaussian distributed measurement noise x 1 À x 2 is also Gaussian, with variance equal to s 2 x1 + s 2 x2 .
We can fit the power spectral density of x 1 À x 2 with the aliased Lorentzian expression for discretely sampled times series (Berg-Sørensen and Flyvbjerg, 2004). Because the measurement noise is Gaussian, we add a noise term to the fitting function. Thus we can estimate the noise of each trajectory by integrating the noise term over the frequency domain. From this, we estimate standard deviations of noise in x 1 and x 2 as 3:1G0:3 nm and 3:8G0:4 nm, respectively.

Data acquisition
We performed experiments using a feedback optical tweezer, whose details are given in section and have been described in previous work (Kumar and Bechhoefer, 2018a). Briefly, we trap a silica bead of 1.5 m diameter using an optical tweezer, which creates a harmonic well without feedback. By applying feedback, we change the shape of the potential to a double well along one of the axes. We measure the position of a bead with two different quadrant photodiodes (QPD) simultaneously to give us two trajectories (x 1 and x 2 ) with two different values of signal-to-noise (SNR) as explained in the section). One detector has high SNR and is used for feedback to create the desired virtual potential (Jun and Bechhoefer, 2012); the other has an adjustable SNR and is used to explore inferences from measured signals with lower SNR. We reduce the SNR in the other detector by placing neutral density (ND) filters of increasing optical density (OD) in front of it. Thus, we can use SKIPPER on the same trajectory over two different experimental SNRs and compare performance. We estimate the measurement noise and SNR in each detector from the noise floor of the power spectrum (see section).

QUANTIFICATION AND STATISTICAL ANALYSIS
Concretely, our goal is to use noisy positional measurements, y 1:N , to infer all unknowns: 1) the potential at each point in space, Uð ,Þ (with UðxÞ denoting the potential evaluated at x); 2) the friction coefficient, z; ll OPEN ACCESS iScience 25, 104731, September 16, 2022 iScience Article 3) the magnitude of the measurement noise, s 2 (under a Gaussian noise model); and 4) the actual position at each time, x 1:N . Toward achieving our goal, we construct a joint posterior probability distribution over all unknowns. As our posterior does not admit an analytic form, we devise an efficient Monte Carlo strategy to sample from it.

Dynamics
We describe the dynamics of the particle with an overdamped Langevin equation (Zwanzig, 2001), where xðtÞ is the possibly multidimensional position coordinate at time t; _ xðtÞ is the velocity; f ðxÞ is the force at position xðtÞ; and z is the friction coefficient. The forces acting on the particle include positional forces f ðxÞ expressed as the gradient of a conservative potential, f ðxÞ = À VUðxÞ. The stochastic (thermal) force, rðtÞ, is defined as follows: where C ,D denotes an ensemble average over realizations, T is the temperature of the bath and k is Boltzmann's constant. Under a forward Euler scheme (LeVeque, 2007) for Equation (1a) with time points given by t n = nDt, each position, given its past realization, is sampled from a normal distribution x n + 1 jx n ; f ð , Þ; z $ N x n + Dt z f ðx n Þ; g 2 I :

(Equation 3)
In words, ''the position x n + 1 given quantities x n ; f ð ,Þ; and z is sampled from a Normal distribution with mean x n + Dt z f ðx n Þ and variance g 2 = 2DtkT z I.'' As is typical for experimental setups, we use a Gaussian noise model and write y n jx n ; s 2 $ N À x n ; s 2 I Á :

(Equation 4)
In words, the above reads ''y n given quantities x n ; s 2 is drawn from a normal.'' Here s 2 is the measurement noise variance. In Equation (4), the measurement process is instantaneous, i.e., assumed to be faster than the dynamical time scales. Our choice of Gaussian measurement model here can be modified at minimal computational cost (e.g., (Hirsch et al., 2013)) if warranted by the data, provided the final measurement noise model is stationary and each measurement depends only on the position at that time level.
Important considerations dictate the prior on the potential. First, the potential may assume any shape (and, as such, is modeled nonparametrically) although it should be smooth (i.e., spatially correlated). A Gaussian process (GP) prior (Williams and Rasmussen, 2006) allows us to sample continuous curves with covariance provided by a pre-specified kernel. However, naive GP prior implementations are computationally prohibitive, with time and memory requirements scaling as the number of data points cubed (Bryan IV et al., 2020;Williams and Rasmussen, 2006 where U i is the potential of region i and p i is the fraction of the time that the particle spent in region i. A limitation of the Boltzmann method is that the presence of measurement noise implies that the fraction of time that a particle was seen in a region, p i , does not equal the fraction of time that the particle was in the region. As a consequence, measurement noise will smear the shape of the inferred potential over the range of the measurement noise (see Figure 5).