High-Performance Reconstruction of Microscopic Force Fields from Brownian Trajectories
Laura Pérez García, Jaime Donlucas Pérez, Giorgio Volpe, Alejandro V. Arzola, Giovanni Volpe
aa r X i v : . [ phy s i c s . d a t a - a n ] A ug High-Performance Reconstruction of Microscopic Force Fields from BrownianTrajectories
Laura P´erez Garc´ıa, Jaime Donlucas P´erez, Giorgio Volpe, Alejandro V. Arzola, and Giovanni Volpe Instituto de F´ısica, Universidad Nacional Aut´onoma de M´exico,Apdo. Postal 20-364, 01000 Cd. M´exico, Mexico Department of Chemistry, University College London,20 Gordon Street, London WC1H 0AJ, United Kingdom, EU Department of Physics, University of Gothenburg, 41296 Gothenburg, Sweden, EU (Dated: August 17, 2018)The accurate measurement of microscopic force fields is crucial in many branches of science andtechnology, from biophotonics and mechanobiology to microscopy and optomechanics [1–4]. Theseforces are often probed by analysing their influence on the motion of Brownian particles [4–7]. Here,we introduce a powerful algorithm for microscopic Force Reconstruction via Maximum-likelihood-estimator (MLE) Analysis (FORMA) to retrieve the force field acting on a Brownian particle fromthe analysis of its displacements. FORMA yields accurate simultaneous estimations of both theconservative and non-conservative components of the force field with important advantages overestablished techniques, being parameter-free, requiring ten-fold less data and executing orders-of-magnitude faster. We first demonstrate FORMA performance using optical tweezers [4]. We thenshow how, outperforming any other available technique, FORMA can identify and characterise stableand unstable equilibrium points in generic extended force fields. Thanks to its high performance,this new algorithm can accelerate the development of microscopic and nanoscopic force transducerscapable of operating with high reliability, speed, accuracy and precision for applications in physics,biology and engineering.
In many experiments in biology, physics, and materials science, a microscopic colloidal particle is used to probelocal forces [1–4]; this is the case, for example, in the measurement of the forces produced by biomolecules, cells,and colloidal interactions. Often particles are held by optical, acoustic, or magnetic tweezers in a harmonic trappingpotential with stiffness k so that a homogeneous force acting on the particle results in a displacement ∆ x fromthe equilibrium position and can therefore be measured as k ∆ x . To perform such measurement, it is necessary todetermine the value of k , which is often done by measuring the Brownian fluctuations of the particle around its stableequilibrium position. This is achieved by measuring the particle position as a function of time, x ( t ), and then usingsome calibrations algorithms; the most commonly employed techniques are the potential [5], the power spectral density(PSD) [6], and the auto-correlation function (ACF)[7] analyses (see Methods for details). The first method samplesthe particle position distribution ρ ( x ), calculates the potential using the Boltzmann factor, and then fits the value of k ; this method requires a series of independent particle positions acquired over a time much longer than the systemequilibration time to sample the probability distribution, and depends on the choice of some analysis parameters suchas the size of the bins. The latter two methods respectively calculate the PSD and ACF of the particle trajectory inthe trap and fit them to their theoretical form to find the value of k ; both methods require a time series of correlatedparticle positions at regular time intervals with a sufficiently short timestep ∆ t , and depend on the choice of someanalysis parameters that determine how the fits are made.FORMA exploits the fact that in the proximity of an equilibrium position the force field can be approximated bya linear form [4, 8] and, therefore, optimally estimated using a linear MLE [9, 10]. This presents several advantagesover the methods mentioned above. First, it executes much faster because the algorithm is based on linear algebra forwhich highly optimised libraries are readily available. Second, it requires less data and therefore it converges fasterand with smaller error bars. Third, it has less stringent requirements on the input data, as it does not require aseries of particle positions sampled at regular time intervals or for a time long enough to reconstruct the equilibriumdistribution. Fourth, it is simpler to execute and automatise because it does not have any analysis parameter to bechosen. Fifth, it probes simultaneously the conservative and non-conservative components of the force field. Finally,since it does not need to use the trajectory of a particle held in a potential, it can identify and characterise bothstable and unstable equilibrium points in extended force fields, and therefore it is compatible with a broader range ofpossible scenarios where a freely-diffusing particle is used as a tracer, e.g., in microscopy and rheology. t n x n x n f n − kx n a b Figure 1 . FORMA: Force Reconstruction using MLE Analysis.
Schematic of FORMA in 1D for a particle held in an opticaltweezers. ( a ) While a particle is held within a harmonic optical trap generated by an optical tweezers (the green background illustrates thedepth of the potential), samples x n of its trajectory (solid line) are acquired at times t n . ( b ) FORMA exploits the fact that the stiffness k of the optical tweezers is related to the correlation between x n and the friction force f n acting on the particle (each dot represents adifferent ( x n , f n ) pair). FORMA uses an MLE to quickly, precisely, and accurately estimate this correlation and, thus, k . The spread ofthe dots around the linear regression line − kx n provides information about the diffusion coefficient D of the particle. To introduce the algorithm in the simplest 1D situation, we start by considering a spherical microparticle of radius R immersed in a liquid with viscosity η , at temperature T , and held in a harmonic confining potential of stiffness k .Experimentally, we have used a standard optical tweezers using a single focused laser beam to create a harmonic trap;in this configuration, the motion along each dimension is independent and can be treated separately as effectively1D [4]. The details of the experimental setup are described in the Methods and supplementary Fig. S1. Briefly, wehave employed a focused laser beam (Gaussian profile, linear polarisation, wavelength 532 nm, power at the sample0 . R = 0 . ± . µ m) in an aqueous solution. We have tracked theparticle position using video microscopy [11] with a frame frequency f s = 4504 . − , corresponding to a samplingtimestep ∆ t = 0 .
222 ms. The corresponding overdamped Langevin equation is [12]:˙ x = − kγ x + √ D w, (1)where x ( t ) is the 1D particle position, γ = 6 πηR , D = k B Tγ , and w ( t ) is a 1D white noise. We assume to have N measurements of the particle displacement ∆ x n at position x n during a time ∆ t n , where n = 1 , ..., N (note that thetime intervals do not need to be equal). Discretizing Eq. 1, we obtain that the average viscous friction force in the n -th time interval is f n = γ ∆ x n ∆ t n = − k x n + σ w n , (2)where σ n = q Dγ ∆ t n and w n is a Gaussian random number with zero mean and variance ∆ t − n [12]. The centralobservation is that Eq. 2 is a linear regression model [9, 10], whose parameters, k and σ , can therefore be optimallyestimated with a maximum likelihood estimator from a series of observations of the dependent ( f n ) and independent( x n ) variables, as schematically illustrated in Fig. 1. The detailed derivation is provided in the Methods for thegeneral case. The MLE estimation of the trap stiffness is then k ∗ = P n x n f n P n x n . (3)Eq. 3 is indeed a very simple expression that can be executed extremely fast using standard highly-optimised linearalgebra libraries, such as LAPACK [13], which is incorporated in most high-level programming languages, includingMatLab and Python. Using the fact that f n + k ∗ x n = σ w n (Eq. 2), we can estimate the diffusion coefficient as D ∗ = 1 N X n ∆ t n γ [ f n + k ∗ x n ] (4)and compare it to the expected value D , which provides an intrinsic quantitative consistency check for the quality ofthe estimation.In Fig. 2a, we show the estimation of the trap stiffness as a function of the number of samples (orange line). Alreadywith as little as 10 samples (corresponding to a total acquisition time of ∼ . < < samples ( ∼ D (Fig. 2b), which indeed converges to the expected value (dashed line) already for 10 samples with a 3%error. Even for the highest number of samples, the algorithm execution time is in the order of a few ms on a laptopcomputer (Fig. 2c). We have further verified these results on simulated data with physical parameters equal to theexperimental ones (Figs. 2d-f). These simulations are in very good agreement with the results of the experiments and,given that the value of k is fixed and known a priori (dashed line in Fig. 2d), they demonstrate the high accuracyof FORMA estimation: FORMA converges to the ground truth value of k for about 10 samples with 10% relativeerror, which, as in the experiments, reduces to 2% for 10 samples.In Fig. 2, we also compare the performance of FORMA with other established methods typically used in thecalibration of optical tweezers, i.e. the potential, PSD, and ACF analyses [4–7] (the details and parameters used forthese analysis are provided in the Methods). Overall, these results show that FORMA is more precise than othermethods when estimating k for a given number of samples, as the other methods typically need 10 to 100 times moredata points to obtain comparable relative errors (Figs. 2a, 2b, 2d, and 2e). FORMA also executes faster by one totwo orders of magnitude than the other methods (Figs. 2c and 2f). The relative errors of the FORMA estimation arealso typically smaller, being 2% for 10 samples, while the potential, PSD and ACF analyses achieve values largerthan 6%, 7% and 20%, respectively. FORMA is also more accurate in estimating the value of k , as can been seencomparing experiments and simulations (Fig. 2a and 2d): while the potential and ACF analyses also converge tothe correct k value, the PSD analysis introduces a significant bias in the estimation (around 8%). Although for thespecific case of the potential analysis (green lines), FORMA’s performance can be considered a marginal improvementin terms of accuracy, it actually provides access to additional information that the potential analysis does not provide,namely the estimation of the diffusion coefficient D . Nonetheless, FORMA is significantly more accurate than thePSD analysis (blue lines), whose estimated k and D present a systematic error; with experience this error can bereduced by tweaking the fitting PSD range, although this process can still be tricky without knowing a priori thevalue of k . Finally, FORMA is also significantly more precise and about two orders of magnitude faster than the ACFanalysis (pink lines), which in fact is the least precise and the slowest method with 20% relative error and 100 msexecution time for 10 samples, while FORMA has 2% relative error and 2 ms execution time for the same numberof samples.We now generalise FORMA to the 2D case. Beyond a conservative component that is also present in 1D forcefields, 2D force fields can also feature a non-conservative component [8]; as we will see, FORMA is able to estimate k [ p N µ m − ] FORMApotpsdacf δ k / k [ % ] D [ µ m s − ] δ D / D [ % ] N s -1 t c [ m s ] N s experiments simulations abc def Figure 2 . High accuracy and performance of FORMA compared to alternative techniques.
Experimentally determined values of ( a )the trap stiffness k and its relative error δk/k , ( b ) the diffusion coefficient D and its relative error δD/D , and ( c ) the computational executiontime t c for FORMA (orange lines); and ( d - f ) corresponding results from numerical simulations. The comparisons with potential (green lines), PSD(blue lines), and ACF (pink lines) analyses show that FORMA converges faster (i.e. for smaller N ), is more precise (i.e. smaller relative errors),is more accurate (i.e. it converges to the expected value represented by the black dashed line), and executes faster than the other methods. Inall cases, we have acquired/simulated 24 trajectories of the motion of a spherical microparticle with radius R = 0 . µ m in an aqueous medium ofviscosity η = 0 . f s = 4504 . − . The relative errors are obtained as the standard deviations of the estimationsover the 24 trajectories. The execution times are measured using a MatLab implementation of the algorithms on a laptop computer (MacBookAir, 2,2 GHz Intel Core i7, 8 GB 1600 MHz DDR3). both simultaneously, differently from most other methods. The overdamped Langevin equation is now best writtenin vectorial form as ˙ r = 1 γ F ( r ) + √ D w (5)where F ( r ) is a force field and w is a vector of independent white noises. F ( r ) can be expanded in Taylor seriesaround as F ( r ) = F + J r + o ( r ), where F = F ( ) is the force and J = J ( ) the Jacobian at . If we assumethat is an equilibrium point, then F = ; with this assumption, the results in the following become much simplerwithout loss of generality, as this is equivalent to translating the experimental reference frame so that it is centred atthe equilibrium position (we discuss below how to proceed if the equilibrium position is not known). Analogously tothe 1D case explained above, using Eq. 5 the average friction force in the n -th time interval is f n = γ ∆ r n ∆ t n = J r n + σ w n , (6)where w n is an array of independent random numbers with zero mean and variance ∆ t − n . Eq. 6 is again a linearregression model and therefore the MLE estimator of J is given by J ∗ = (cid:2) r T r (cid:3) − r T f , (7)where r = ( r n ) and f = ( f n ) are matrices with N × × F ∗ ( r ) = J ∗ r , where we use the estimated Jacobian [8] (Eq. 7). This is a linear formthat results from the superposition of a conservative harmonic potential (which is characterised by its stiffnesses k ∗ and k ∗ along the principle axes, and the orientation θ ∗ of the principle axes with respect to the Cartesian axes) anda non-conservative rotational force field (which is characterised by its angular velocity Ω ∗ ). Some examples of thisdecomposition are shown in supplementary Fig. S2. It is possible to obtain these two components directly from theJacobian, by separating it into its conservative and non-conservative part as J ∗ = J ∗ c + J ∗ r , making use of the factthat they are respectively symmetric and antisymmetric. The conservative part is J ∗ c = 12 ( J ∗ + J ∗ T0 ) = R ( θ ∗ ) (cid:20) − k ∗ − k ∗ (cid:21) R − ( θ ∗ ) , (8)where R ( θ ) is a rotation matrix that diagonalises J c and whose principal axes correspond to the eigenvectors cor-responding to the principle axes of the harmonic potential and the stiffnesses along these axes correspond to theeigenvalues with a minus sign. The non-conservative (rotational) part is J ∗ r = 12 ( J ∗ − J ∗ T0 ) = (cid:20) − γ Ω γ Ω 0 (cid:21) (9)and, since it is invariant under a rotation of the reference system, can be simply estimated asΩ ∗ = 12 γ (cid:2) J ∗ , − J ∗ , (cid:3) . (10)To demonstrate this 2D version of FORMA at work, we have used it to estimate the transfer of orbital andspin angular momentum to an optically trapped particle. In fact, orbital and spin angular momentum can make atransparent particle orbit, even though the precise angular-momentum-transfer mechanisms can be very complex whenthe beam is focalised, depending on the size, shape and material of the particle, as well as on the size and shape of thebeam [14–18]. We employ the same setup and microparticle as for the results presented in Fig. 2 (see Methods andsupplementary Fig. S1), using a spatial light modulator to generate Laguerre-Gaussian (LG) beams carrying orbitalangular momentum (OAM) to trap the particle and a quarter-wave plate (QWP) to switch their polarisation statefrom linear to positive or negative circular polarisation. The results of the force field reconstruction are presented inFig. 3. In Figs. 3d and 3g, we employ a linearly polarized LG beam with topological charge l = − , − , , , l = 0 corresponds to the standard Gaussian beam already employed in Fig. 2). As already observed in previousexperiments, we also measure that Ω ∗ is proportional to l (Ω ∗ ≈ . − l ), which follows the quantisation of theOAM [19]. We have verified that these values are in good quantitative agreement with those obtained using the morestandard cross-correlation function (CCF) analysis, which is an extension of the ACF analysis that permits one todetect the presence of non-conservative force fields [8]. These non-conservative rotational force fields are very small,produce an almost imperceptible bending of the force lines when comparing LG with LG (insets of Fig. 3g), and,therefore, cannot be detected by directly counting rotations of the particle around the beam axis. To further test ourresults, we have then changed the polarisation state of the beam switching it to positive circular polarisation (Figs. 3eand 3h), which introduces an additional spin angular momentum (SAM), so that Ω ∗ ≈ . − ( l + 1), recovering theSAM quantisation. This relation suggest that the OAM and SAM contributes equally to the rotation of the particle. r n r n − ( J r n ) x ( f n ) x x n y n F ( r ) r eq k x , k y [ p N µ m − ] k x k y LG -2 LG -1 LG LG LG -5.4-3.6-1.801.83.65.4 Ω [ s − ] LG -2 LG -1 LG LG LG FORMAccf LG -2 LG -1 LG LG LG LG LG linear polarisation circular (+) pol. circular ( − ) pol. a b cd e fg h i Figure 3 . FORMA in 2D: Measurement of the non-conservative force-field component. ( a - c ) Schematic of FORMA in 2D for aparticle held in an optical tweezers: ( a ) Samples r n of a particle trajectory (solid line) held in an optical tweezers (the green background illustratesthe depth of the potential) are acquired at times t n . ( b ) FORMA estimates the Jacobian J of the force field from the relation between r n and f n using a 2D MLE. In the schematic, we represent only the estimation of the first row of J , which is related to the x -component of f n ; the completegraph cannot be represented because it is 4D. ( c ) Using this information, FORMA reconstructs the force field around the equilibrium point r eq (see also supplementary Fig. S2). ( d - f ) Stiffnesses k x and k y , and ( g - i ) angular velocity Ω of a Brownian particle optically trapped by ( d , g ) alinearly polarised, ( e , h ) circularly (+) polarised, and ( f , i ) circularly ( − ) polarised Laguerre-Gaussian (LG) beam with l = − , − , , ,
2. ( g - i )The results of FORMA (orange circles) agree well with the results of the CCF analysis (pink triangles). The insets in (b) show the force fields forthe case of a Gaussian beam (LG ), which is purely conservative ( l = 0), and a beam with a charge l = 2 of orbital angular moment (LG ), whichfeatures a non-conservative component that induces a mild bending of the arrows. In all cases, we have acquired trajectories of the motion of aspherical particle with radius R = 0 . µ m in an aqueous medium of viscosity η = 0 . f s = 4504 . − , and used25 windows of 10 samples for the analysis; the error bars are the standard deviations over these 25 measurements. We finally changed the polarisation state of the beam to negative circular polarisation (Figs. 3f and 3i) obtainingΩ ∗ ≈ . − ( l − a priori the positions ofthe equilibria might be unknown, for example when exploring extended potential landscapes. FORMA can be furtherrefined to address this problem and determine the value of F in Eq. 5, which permits one to identify equilibriumpoints when F = . We obtain the following estimator (see Methods for detailed derivation):[ J ∗ F ∗ ] = (cid:0) ˜ r T ˜ r (cid:1) − ˜ r T f , (11)where ˜ r = [ r 1 ] with a column vector constituted of N ones. Having the trajectory of a particle moving in anextended potential landscape, it is possible to use Eq. 11 to simultaneously identify the equilibrium points, classifytheir stability, and characterise their local force field: For each position in the potential landscape, the parts of thetrajectory that fall within a radius a smaller than the characteristic length over which the force filed varies from thisposition are selected and analysed with FORMA; if F ∗ ≈ , then this position is an equilibrium point and J ∗ permitsus to determine the local force field. -0.6 -0.4 -0.2 0 0.2 0.4 0.6 x [ µ m] -4-3-2-10 U ( x ) [ k B T ] x ∗ x ∗ x ∗ x ∗ x ∗ -0.6 -0.4 -0.2 0 0.2 0.4 0.6 x [ µ m] -0.4-0.200.20.4 y [ µ m ]
100 fN ab Figure 4 . Reconstruction of stable and unstable equilibrium points in a multiwell optical potential. ( a ) Multiwell opticalpotential generated by two focused Gaussian beams slightly displaced along the x -direction. FORMA identifies three stable ( x ∗ , x ∗ , x ∗ ;full circles) and two unstable ( x ∗ , x ∗ ; empty circles) equilibrium points, and measures their stiffness (orange solid and dashed lines). Thecorresponding x -potential obtained from the potential method is shown by the green solid line (the green shaded area represents onestandard deviation obtained repeating the experiment 20 times). ( b ) 2D plot of the force field measured with FORMA (arrows) and of thepotential measured with the potential analysis (background colour). The stable and unstable equilibrium points are indicated by the fulland empty circles respectively. We have acquired trajectories of the motion of a spherical particle with radius R = 0 . µ m in an aqueousmedium of viscosity η = 0 . f s = 4504 . − , and used 20 windows of 4 . · samples for the analysis. We have applied this procedure to identify the equilibrium points in a multistable potential, which we realisedfocusing two slightly displaced laser beams obtained using a spatial light modulator (see Methods and supplementaryFig. S1). Similar configurations have been extensively studied as a model system for thermally activated transitionsin a bistable potentials [20, 21]; however, the presence of additional minima other than the two typically expected hasbeen recognized only recently due to their weak elusive nature [22]. The results of the reconstruction are shown inFig. 4. In Fig. 4a, we use the potential analysis to determine the potential (green solid line with error bars denotedby the shaded area) by acquiring a sufficiently long trajectory so that the particle has equilibrated and fully exploredthe region of interest; this potential appears to be bistable with two potential minima (stable equilibrium points) anda potential barrier between them corresponding to an unstable equilibrium point. When we use FORMA, we identifyfive equilibria, and classify them as stable ( x ∗ , x ∗ , x ∗ ; full circles) and unstable ( x ∗ , x ∗ ; empty circles); importantly,we are able to clearly resolve the presence of additional equilibrium points within the potential barrier ( x ∗ , x ∗ , x ∗ ).For each point, FORMA provides the respective stiffness (the corresponding harmonic potentials are shown by theorange lines), which is negative for unstable equilibrium points. This highlights one of the additional key advantagesof the MLE algorithm: It can also measure the properties of the force field around unstable equilibrium points, thanksto the fact that it does not require an uninterrupted series of data points or a complete sampling of the equilibriumdistribution, which are difficult or impossible near an unstable equilibrium. Fig. 4b shows a 2D view of the potential,where we used FORMA to estimate the force field on a 2D grid and to identify the stable and unstable equilibria.The background colour represents the potential reconstructed using the equilibrium distribution, which shows a goodagreement with the results of FORMA, but does not allow to clearly identify the presence of the additional equilibriain the potential barrier. Figure 5 . Reconstruction of the equilibrium points in a speckle pattern. ( a ) The intensity of the speckle (green background,laser wavelength λ = 532 nm) is approximately proportional to the potential depth of the optical potential felt by a particle whose size(particle diameter 1 . ± . µ m) is similar to the speckle characteristic size (2 . µ m) [23, 24]. FORMA identifies several stable (full circles)and unstable (empty circles) equilibrium points, and measures the orientation of the principal axes ( θ ), the stiffnesses along them ( k and k ), and angular velocity (Ω) (see supplementary table S1 for the measured values). We have acquired trajectories of the motion of aspherical particle with radius R = 0 . ± . µ m in an aqueous medium of viscosity η = 0 . f s = 600 s − ;since we cannot expect the particle to spontaneously diffuse over the whole speckle field during the time of the experiment, we have placedthis particle in 25 positions within the speckle field and let it diffuse each time acquiring 2 · samples for the analysis. ( b - d ) Exampleof reconstructed force fields around ( b ) a stable point, ( c ) an unstable point with a significant rotational component (indicated by thearrow), and ( d ) two stable points with a saddle in between; the grey arrows plot the 2D force field measured with FORMA and are scaledby a different factor in each plot. Finally, we can also use FORMA to study larger extended force fields, such as the random optical force fieldsgenerated by speckle patterns [24–28]. Speckles are complex interference patterns with well-defined statistical prop-erties generated by the scattering of coherent light by disordered structures [29]; the equilibrium positions are notknown a priori due to their random appearance, the configuration space is virtually infinite, and there can be a non-conservative component. Thus, the challenge is at least two-fold. First, the configuration space is virtually infiniteand, therefore, cannot be sampled by a single trajectory in any reasonable amount of time. Second, the force fieldcan present a non-conservative component; in fact, while several works have demonstrated micromanipulation withspeckle patterns [24–28], none of them has so far achieved the experimental characterisation of this nonconservativenature of the optical forces. To study this situation, we have employed a speckle light field generated using a secondoptical setup (see Methods and supplementary Fig. S3). A portion of the resulting speckle field is shown by the greenbackground in Fig. 5a. To sample the force field, we have acquired the trajectories of a particle in various regions ofthe speckle field: In each of these trajectories the particle typically explores the regions surrounding several contiguousstable equilibrium points by being metastably trapped in each of them while still being able to cross over the potentialbarriers separating them [23]. These trajectories cannot be used in the potential analysis because they do not providea fair sampling of the position space. Nevertheless, they can be used by FORMA to identify the equilibrium points,which are shown in Fig. 5a by the full circles (stable points) and empty circles (unstable and saddle points), and todetermine the force field around them (see supplementary table S1 for the measured values): For example, in Fig. 5bwe show a stable point, in Fig. 5c an unstable point with a significant rotational component, and in Fig. 5d a seriesof two stable points with a saddle between them in a configuration reminiscent of that explored in Fig. 4.In conclusion, we have introduced FORMA: a new, powerful algorithm to measure microscopic force fields usingthe Brownian motion of a microscopic particle based on a linear MLE. We first introduced the 1D version of FORMA;we quantitatively compared it to other standard methods, showing that it needs less samples, it has smaller relativeerrors, it is more accurate, and it is orders-of-magnitude faster. We then introduced the 2D version of FORMA; weshowed that it can also measure the presence of a non-conservative force field, going beyond what can be done by theother methods. Finally, we applied it to more general force-field landscapes, including situations where the forces aretwo shallow to achieve long-term trapping: we used FORMA to identify the equilibrium points; to classify them asstable, unstable and saddle points; and to characterise their local force fields. Overall, we have shown that, requiringless data, FORMA can be applied in situations that require a fast response such as in real-time applications and in thepresence of time-varying conditions. FORMA can be straighforwardly extended also to measure flow fields and also3D force fields. Thanks to the fact that this algorithm is significantly faster, simpler and more robust than commonlyemployed alternatives, it has the potential to accelerate the development of force transducers capable of measuringand applying forces on microscopic and nanoscopic scales, which are needed in many areas of physics, biology andengineering.
METHODS
Potential analysis.
The potential method [4, 5] relies on the fact that the harmonic potential has the form U ( x ) = kx and the associated position probability distribution of the particle is ρ ( x ) = ρ exp (cid:16) − U ( x ) k B T (cid:17) . Therefore,by sampling ρ ( x ) it is possible to reconstruct U ( x ) = − k B T ln( ρ ( x )) and F ( x ) = − dU ( x ) dx . The value of k is finallyobtained by fitting F ( x ) to a linear function. The potential method requires a series of particle positions acquiredover a time long enough that the system has equilibrated as well as the choice of the number of bins and of the fittingalgorithm to be employed in the analysis: here, we used 100 bins equally spaced between the minimum and maximumvalue of the particle position, and a linear fitting. Power spectral density (PSD) analysis.
The PSD method [4, 6] uses the particle trajectory in the harmonictrap, x ( t ), to calculate the PSD, P ( f ) = D π ( f + f ) − , where f c = (2 πγ ) − k is the harmonic trap cutoff frequency, γ is the particle friction coefficient, and D is its diffusion coefficient. It then fits this function to find the value of k ,and therefore the harmonic force field surrounding the particle. The PSD method requires a time series of correlatedparticle positions at regular time intervals with a sufficiently short timestep ∆ t . It also requires to choose how thePSD is calculated (e.g., use of windowing and binning [6]) and the frequency range over which the fitting is made:here, we performed the PSD fitting over the frequency range between 5 times the minimum measured frequency andhalf the Nyquist frequency without using windowing and binning. Auto-correlation function (ACF) analysis.
The ACF method [4, 7] calculates the ACF of the particle position, C ( τ ) = k B Tk exp (cid:16) − k | τ | γ (cid:17) , where k B is the Boltzmann constant and T is the absolute temperature. It then fits thisfunction to find the value of k , and therefore the harmonic force field surrounding the particle. Like the PSD method,also the ACF method requires a time series of correlated particle positions at regular time intervals. It also requiresto choose over which range to perform the fitting: here, we have performed the fitting over the values of the ACFlarger than 1% its maximum. Experimental setups.
For the single-beam (Figs. 2 and 3) and two-beam (Fig. 4) experiments, we used thestandard optical tweezers shown in supplementary Fig. S1 [4, 28, 30]. An expanded 532-nm-wavelength laser beam(power at the sample 0 . l = − , − , , , − ) polarised (Fig. 3f and 3i). These experiments are performed withspherical silica microparticles with radius R = 0 . µ m in an aqueous medium of viscosity η = 0 . f s = 4504 . − .For the speckle experiments (Fig. 5), we used the SLM in the image plane, as shown in supplementary Fig. S3: The532-nm laser beam is reflected by the SLM, which projects a random phase (with uniform distribution of values in(0 , π )) in every domain of 6 × D i = 1 . ± .
05 mm in the Fourier plane of the first telescope. The mean intensity of the speckle is 0 .
015 mW µ m − .These experiments are performed with spherical polystyrene microparticles with radius R = 0 . ± . µ m in anaqueous medium of viscosity η = 0 . · frames each at a sampling frequency f s = 600 s − . Detailed derivation of FORMA in 2D (Eq. 11).
Here we derive FORMA in its most general form presentedin the article (Eq. 11). The average friction force in the n -th time interval is f n = γ ∆ r n ∆ t n = F + J r n + r k B T γ ∆ t n w n , (12)which can be rewritten explicitly as (cid:20) f x,n f y,n (cid:21) ≈ (cid:20) J , J , F ,x J , J , F ,y (cid:21) xy + r k B T γ ∆ t (cid:20) w x w y (cid:21) , (13)Assuming to have N measurements, we introduce the vectors f = γ ∆ x ∆ t γ ∆ y ∆ t ... ...γ ∆ x n ∆ t n γ ∆ y n ∆ t n ... ...γ ∆ x N ∆ t N γ ∆ y N ∆ t N and ˜ r = x y ... ... ...x n y n ... ... ...x N y N , the MLE is given by Eq. 11, i.e., [ J ∗ F ∗ ] = (cid:20) J ∗ , J ∗ , F ∗ ,x J ∗ , J ∗ , F ∗ ,y (cid:21) = (cid:0) ˜ r T ˜ r (cid:1) − ˜ r T f and the estimated particle diffusivity along each axis can be calculated from the residual error of the MLE D ∗ x = 1 N N X n =1 ∆ t n γ (cid:0) f x,n − J ∗ , x n − J ∗ , y n − F ∗ ,x (cid:1) D ∗ y = 1 N N X n =1 ∆ t n γ (cid:0) f y,n − J ∗ , x n − J ∗ , y n − F ∗ ,y (cid:1) (14) Codes.
We provide the MatLab implementations of the key functionalities of FORMA:1. 1D version of FORMA: function forma1d.m, script test forma1d.m to execute it, and set of test data forma1d.mat.This code estimated the values of k ∗ and D ∗ assuming that the equilibrium position is at x = 0 and implementingEqs. (3) and (4).2. 2D version of FORMA: function forma2d.m, script test forma2d.m to execute it, and set of test data forma2d.mat.This code estimates the value of k ∗ , k ∗ , θ ∗ , and Ω ∗ assuming that the equilibrium position is at r = and imple-menting Eqs. (7), (8), and (10). [1] Aspelmeyer, M., Kippenberg, T. J. & Marquardt, F. Cavity optomechanics. Rev. Mod. Phys. , 1391–1452 (2014).[2] Roca-Cusachs, P., Conte, V. & Trepat, X. Quantifying forces in cell biology. Nat. Cell Biol. , 742–751 (2017).[3] Moerner, W. E. Single-molecule spectroscopy, imaging, and photocontrol: Foundations for super-resolution microscopy(Nobel lecture). Ang. Chemie Int. Ed. , 8067–8093 (2015). [4] Jones, P. H., Marag`o, O. M. & Volpe, G. Optical tweezers: Principles and applications (Cambridge University Press,2015).[5] Florin, E.-L., Pralle, A., Stelzer, E. H. K. & H¨orber, J. K. H. Photonic force microscope calibration by thermal noiseanalysis.
Appl. Phys. A , S75–S78 (1998).[6] Berg-Sørensen, K. & Flyvbjerg, H. Power spectrum analysis for optical tweezers. Rev. Sci. Instrumen. , 594–612 (2004).[7] Bechhoefer, J. & Wilson, S. Faster, cheaper, safer optical tweezers for the undergraduate laboratory. Am. J. Phys. ,393–400 (2002).[8] Volpe, G., Volpe, G. & Petrov, D. Brownian motion in a nonhomogeneous force field and photonic force microscope. Phys.Rev. E , 061118 (2007).[9] Neter, J., Wasserman, W. & Kutner, M. H. Applied linear regression models (Irwin Homewood, IL, 1989).[10] DeGroot, M. H. & Schervish, M. J.
Probability and statistics (Pearson Education, 2012).[11] Crocker, J. C. & Grier, D. G. Methods of digital video microscopy for colloidal studies.
J. Colloid Interface Sci. ,298–310 (1996).[12] Volpe, G. & Volpe, G. Simulation of a Brownian particle in an optical trap.
Am. J. Phys. , 224–230 (2013).[13] Anderson, E. et al. LAPACK Users’ guide (SIAM, 1999).[14] He, H., Friese, M. E. J., Heckenberg, N. R. & Rubinsztein-Dunlop, H. Direct observation of transfer of angular momentumto absorptive particles from a laser beam with a phase singularity. Phys. Rev. Lett. , 826–829 (1995).[15] Simpson, N. B., Dholakia, K., Allen, L. & Padgett, M. J. Mechanical equivalence of spin and orbital angular momentumof light: An optical spanner. Opt. Lett. , 52 (1997).[16] Zhao, Y., Edgar, J. S., Jeffries, G. D. M., McGloin, D. & Chiu, D. T. Spin-to-orbital angular momentum conversion in astrongly focused optical beam. Phys. Rev. Lett. , 073901 (2007).[17] Albaladejo, S., Marqu´es, M. I., Laroche, M. & S´aenz, J. J. Scattering forces from the curl of the spin angular momentumof a light field. Phys. Rev. Lett. , 113602 (2009).[18] Arzola, A. V., J´akl, P., Chv´atal, L. & Zem´anek, P. Rotation, oscillation and hydrodynamic synchronization of opticallytrapped oblate spheroidal microparticles.
Opt. Express , 16207–16221 (2014).[19] Chen, M., Mazilu, M., Arita, Y., Wright, E. M. & Dholakia, K. Dynamics of microparticles trapped in a perfect vortexbeam , 4919–4922.[20] Han, S., Lapointe, J. & Lukens, J. E. Effect of a two-dimensional potential on the rate of thermally induced escape overthe potential barrier. Phys. Rev. B , 6338–6345 (1992).[21] McCann, L. I., Dykman, M. & Golding, B. Thermally activated transitions in a bistable three-dimensional optical trap. Nature , 785–787 (1999).[22] Stilgoe, A. B., Heckenberg, N. R., Nieminen, T. A. & Rubinsztein-Dunlop, H. Phase-transition-like properties of double-beam optical tweezers.
Phys. Rev. Lett. , 248101 (2011).[23] Volpe, G., Volpe, G. & Gigan, S. Brownian motion in a speckle light field: Tunable anomalous diffusion and selectiveoptical manipulation.
Sci. Rep. , 3936 (2014).[24] Volpe, G., Kurz, L., Callegari, A., Volpe, G. & Gigan, S. Speckle optical tweezers: Micromanipulation with random lightfields. Opt. Express , 18159–18167 (2014).[25] Shvedov, V. G. et al. Selective trapping of multiple particles by volume speckle field.
Opt. Express , 3137–3142 (2010).[26] Shvedov, V. G. et al. Laser speckle field as a multiple particle trap.
J. Opt. , 124003 (2010).[27] Hanes, R. D. L., Dalle-Ferrier, C., Schmiedeberg, M., Jenkins, M. C. & Egelhaaf, S. U. Colloids in one dimensional randomenergy landscapes. Soft Matter , 2714–2723 (2012).[28] Pesce, G. et al. Step-by-step guide to the realization of advanced optical tweezers.
J. Opt. Soc. Am. B , B84–B98(2015).[29] Goodman, J. W. Speckle phenomena in optics: Theory and applications (Roberts and Company Publishers, 2007).[30] Grier, D. G. A revolution in optical manipulation.
Nature , 810–816 (2003).[31] Fax´en, H. Die Bewegung einer starren Kugel langs der Achse eines mit zaher Flussigkeit gefullten Rohres.
Arkiv forMatemetik Astronomi och Fysik , 1–28 (1923). Acknowledgements.
We thank Karen Volke-Sepulveda for useful discussions and Antonio A. R. Neves for criticalreading of the manuscript.
Funding
LPG, JD and AVA acknowledge funding from DGAPA-UNAM (grants PAPIIT IA104917 and IN114517).AVA and GV (Giovanni Volpe) acknowledge support from C´atedra Elena Aizen de Moshinsky. GV (Giovanni Volpe)was partially supported by the ERC Starting Grant ComplexSwimmers (Grant No. 677511).
Contributions
GV (Giovanni Volpe) had the original idea for this method while visiting the National AutonomousUniversity of Mexico (UNAM). LPG and AVA performed most of the experiments and analysed the data; JDL per-formed the experiments reported in Fig. 4. GV (Giorgio Volpe) and GV (Giovanni Volpe) performed the simulations.LPG, GV (Giorgio Volpe), AVA and GV (Giovanni Volpe) discussed the data and wrote the draft of the article. Allauthors revised the final version of the article.2
Competing interests
The authors declare that they have no competing financial interests.
Correspondence
Correspondence and requests for materials should be addressed to Alejandro V. Arzola (alejan-dro@fisica.unam.mx) or Giovanni Volpe (email: [email protected]).3
Supplementary Figure S1 . Holographic optical tweezers for standard optical tweezers, transfer of angular momentum,and multiwell potential.
An expanded laser beam (wavelenght 532 nm) is reflected by a spatial light modulator (SLM, HamamatsuX10468-04), through a telescope (lens L1, f = 400 mm, and lens L2, f = 175 mm), where the diffracted field of interest is discriminatedby the iris diaphragm, a half-wave plate (HWP), a polarising beam splitter (PBS), and a second telescope (T2, magnification M = 3); it isthen reflected by a dichroic mirror (DM), and finally goes through a quarter wave plate (QWP) to be focused by an objective (Carl Zeiss,100 × , oil immersion, 1 .
30 NA) within the sample chamber. The trapped particle is visualised through an imaging system consisting of adiode-based white-light illumination, a condenser (20 × , 0 .
75 NA), the 100 × objective, a tube lens ( f = 200 mm), and a CMOS Camera(Basler acA800-510um). The shape of the beam is controlled by means of the SLM: ( a ) blazed grating phase modulation used to generatea single optical tweezers (Fig. 2); ( b ) blazed grating plus a l = 1 azimuthal phase to generate a LG beam (Fig. 3); ( c ) two identicaltweezers with a phase shift π to generate a multiwell potential (Fig. 4). ( d ) The QWP permits us to switch the polarisation state ofthe beam between linearly polarised ( β = 0, Figs. 3d and 3g), circularly (+) polarised ( β = + π/
4, Figs. 3e and 3h), and circularly ( − )polarised ( β = − π/
4, Figs. 3f and 3i). Supplementary Figure S2 . Examples of force fields and their decomposition into conservative and non-conservativecomponents. ( a ) Isotropic harmonic potential where k = k >
0. ( b ) Elliptical harmonic potential where 0 < k < k . ( c ) Saddlepoint where k < k >
0. ( d ) Unstable point where k < k < θ = 45 ◦ . The full and empty circlesindicate the stable and unstable equilibrium points respectively. Supplementary Figure S3 . Speckle optical tweezers.
Optical micromanipulation with an image speckle. An expandedlaser beam (wavelength 532 nm) is reflected by the spatial light modulator (SLM, Hamamatsu X10468-04), goes through atelescope (lens L1, f = 400 mm, and lens L2, f = 100 mm) where it is spatially filtered by an on-axis iris diaphragm (aperture D i = 1 .
54 mm), passes through a dichroic mirror (DM), is focused by lens L3 ( f = 150 mm) on the back-focal plane ofan objective (Olympus, 40 × , 0 .
65 NA). The particle is visualsed through an imaging system similar to that described insupplementary Fig. S1. ( a ) A uniform distributed random phase mask is projected on the SLM to generate the speckle pattern( c ), and ( b ) a ring-like dynamic phase distribution is used to set the initial position of the particle within the area of interest( d ). Once the particle is in the desired initial position, the SLM projects the random phase distribution to generate the specklepattern, and the particle position is recorded as it diffuses for about 6 seconds, then the speckle pattern is changed for thering-like trap to set the initial position, this process is repeated 500 times. The polystyrene particle is trapped in 2D as it ispushed towards the upper coverslip by the beam radiation pressure. Supplementary Table S1 . Equilibrium points identified by FORMA in the data shown in Fig. 5.
For each point, FORMAprovides the position ( x ∗ eq , y ∗ eq ) (from the lower left corner of Fig. a), its stiffnesses k ∗ and k ∗ along the principle axes, the orientation θ ∗ of the principle axes with respect to the Cartesian axes, and the angular velocity Ω ∗ associated to the non-conservative (rotational)component of the force field. N is the number of measurements of the particle displacements used by FORMA for the estimation. Eq. point x ∗ eq ( µ m) y ∗ eq ( µ m) k ∗ (pN µ m − ) k ∗ (pN µ m − ) θ ∗ (rad) Ω ∗ (s − ) N .
79 8 .
08 0 .
14 0 .
07 358 . .
47 74862 4 .
41 7 .
66 0 .
21 0 .
04 37 . − .
07 1754063 1 .
41 7 .
38 0 . − .
09 14 . .
67 8304 2 .
11 7 .
08 0 . − .
08 115 . − .
69 5625 5 .
61 6 .
72 0 .
11 0 .
07 105 . − .
64 12406 5 .
21 6 .
58 0 . − .
10 108 . .
58 10847 1 .
91 6 .
58 0 . − .
11 62 . − .
14 8418 0 .
75 6 .
22 0 .
10 0 .
08 55 . − .
25 416679 4 .
71 6 .
08 0 .
06 0 .
00 357 . .
24 1314410 3 .
11 6 .
02 0 .
07 0 .
01 49 . − .
43 541111 3 .
81 5 . − . − .
06 31 . − .
87 210112 5 .
61 5 .
78 0 . − .
05 123 . .
59 62613 1 .
51 5 .
78 0 . − .
12 43 . .
09 114614 3 .
21 5 . − . − .
05 12 . .
10 236915 3 .
61 5 . − . − .
07 10 . − .
82 182916 6 .
27 5 .
12 0 .
21 0 .
18 113 . − .
00 12387017 4 .
85 5 .
08 0 . − .
04 3 . .
17 615418 1 .
11 5 .
08 0 .
08 0 .
03 351 . .
10 134019 2 .
27 4 .
74 0 .
16 0 .
14 5 . − .
07 48654720 3 .
43 4 .
60 0 . − .
03 89 . − .
86 481921 5 .
43 4 .
44 0 . − .
09 126 . .
14 537922 4 .
71 4 .
38 0 .
08 0 .
00 98 . − .
32 2601323 1 .
41 4 .
38 0 . − .
06 62 . .
55 1494824 3 .
51 4 . − . − .
04 339 . .
89 392025 1 .
05 4 .
20 0 .
12 0 .
06 94 . − .
53 3200626 6 .
01 3 .
98 0 . − .
06 52 . .
14 118827 3 .
95 3 .
98 0 . − .
12 329 . − .
58 459028 5 .
27 3 .
96 0 .
08 0 .
04 125 . − .
74 1754029 2 .
61 3 .
78 0 . − .
13 34 . − .
89 256130 2 .
41 3 . − . − .
11 55 . − .
07 354631 5 .
67 3 .
42 0 . − .
04 331 . − .
80 356432 1 .
11 3 . − . − .
08 16 . − .
75 122033 4 .
31 3 . − . − .
08 355 . − .
03 250534 7 .
23 3 .
32 0 .
22 0 .
14 26 . − .
30 4239935 3 .
19 3 .
24 0 .
10 0 .
05 72 . − .
17 10076336 2 .
31 2 .
98 0 . − .
05 77 . .
18 344037 5 .
07 2 .
92 0 . − .
02 80 . − .
78 172938 5 .
91 2 .
88 0 . − .
04 71 . .
53 141139 4 .
81 2 . − . − .
10 352 . .
58 140740 5 .
71 2 .
58 0 . − .
07 121 . .
57 144941 3 .
71 2 .
48 0 . − .
01 355 . − .
62 790742 1 .
45 2 .
44 0 .
14 0 .
13 21 . − .
64 5780343 3 .
01 2 .
08 0 . − .
10 106 . − .
62 152944 1 .
01 0 .
76 0 .
16 0 .
07 75 . − .
27 1829245 4 .
69 0 .
72 0 . − .
03 104 . − .
22 1662946 4 .
05 0 .
72 0 .
09 0 .
02 39 . .
21 3438447 5 .
15 0 .
56 0 .
14 0 .
06 110 ..
06 110 .. − ..