Quantitative analysis of non-equilibrium systems from short-time experimental data

Estimating entropy production directly from experimental trajectories is of great current interest but often requires a large amount of data or knowledge of the underlying dynamics. In this paper, we propose a minimal strategy using the short-time Thermodynamic Uncertainty Relation (TUR) by means of which we can simultaneously and quantitatively infer the thermodynamic force field acting on the system and the (potentially exact) rate of entropy production from experimental short-time trajectory data. We benchmark this scheme first for an experimental study of a colloidal particle system where exact analytical results are known, prior to studying the case of a colloidal particle in a hydrodynamical flow field, where neither analytical nor numerical results are available. In the latter case, we build an effective model of the system based on our results. In both cases, we also demonstrate that our results match with those obtained from another recently introduced scheme. Thermal fluctuations play a crucial role in non-equilibrium phenomena at microscopic length scales, making it challenging to analyse and interpret experimental data. Here, the authors demonstrate that the short-time thermodynamic uncertainty relation inference scheme can estimate the entropy production rate for a colloidal particle in time-varying potentials and with background flows determined by the presence of a microbubble.

N on-equilibrium thermodynamics at microscopic length scales is dominated by a fascinating range of phenomena 1 , where thermal fluctuations play a crucial role. These phenomena can now be observed in great detail experimentally, due to the availability and scope of current microscopic manipulation techniques 2 . The interpretation and quantitative analysis of the experimentally available data are however lagging behind these advances, mostly due to the fact that the vast majority of these systems are too complicated to model without making several approximations 3 , despite having far fewer degrees of freedom than their macroscopic counterparts. Even when it is possible to build such simplified models, these are still usually too complicated to solve except sometimes by numerical analysis of specific systems, which however lack general insights 4 . There could also be other factors making the system hard to solve, such as the presence of a background flow, for which the spatial dependence of the flow velocity needs to be known by means of solving the corresponding Navier-Stokes equation; usually a difficult task, especially for unsteady flows 5,6 . In the face of all these challenges, a relevant question is whether it is at all possible to gain any precise quantitative information about a complex nonequilibrium system directly from experimental data, bypassing the first step of either having a known model to compare with or building in simplifying assumptions about the system.
Not surprisingly, this question has aroused a lot of recent interest. Broadly speaking, measurements from experiments can be used to obtain general information about the system, such as identifying that detailed balance is broken and hence the system is out-of-equilibrium [7][8][9] (not always obvious for microscopic systems such as at the cellular level), or to obtain more specific properties of the system such as the rate of dissipation of energy (equivalently the rate of entropy production) [10][11][12][13][14][15][16][17] , the average phase-space velocity field 7,18,19 related to the so-called thermodynamic force field 20,21 or the microscopic forces driving the system 16,19,22 . The motivation for such studies is that if quantitative information about the system can be directly obtained from experimentally observed quantities, then this understanding can be used for building more realistic and experimentally validated models of the system of interest 7,23,24 .
A very informative quantity about a non-equilibrium system is the rate of entropy production. This quantity not only signalswhen it is non-zero-that the system is out of equilibrium, but also provides a quantitative measure of how far from equilibrium a system is and the irreversibility of the dynamics [25][26][27] . In the context of microscopic machines 28 , a quantification of the amount of energy dissipated directly provides information about engine efficiencies [29][30][31] and prescriptions for obtaining optimal operating conditions 32 . The value of the entropy production rate can also be used to obtain information-theoretic quantities of interest 33 , or even information about hidden degrees of freedom 34 . The entropy production rate is also a very robust quantity to measure from the experimental point of view, since it is not so strongly affected by conversion-factor errors in measuring particle positions, as we remark later.
The entropy production rate can be obtained directly from experimental data, at least for systems where it is understood that the underlying dynamics is Markovian, by several means. These include utilizing the Harada-Sasa equality 10 that involves a spectral analysis of trajectory data 35,36 , determining the average steady-state current and steady-state probability distribution from the data 11 , determining the time-irreversibility of the dynamics 27,[37][38][39][40][41] and relatedly determining estimators for the ratio of forward and backward processes directly from the data 14,42,43 . Recent approaches 16,19 also advocate inferring first the microscopic force field from which the entropy production rate can be deduced.
An alternative strategy to direct estimation, is to set lower bounds on the entropy production rate 44-48 by measuring experimentally accessible quantities. One class of these bounds, for example, those based on the thermodynamic uncertainty relation (TUR) 3,48-51 , have been further developed into variational inference schemes, which translate the task of identifying entropy production to an optimization problem over the space of a single projected fluctuating current in the system 15,[52][53][54] . Recently, a similar variational scheme using neural networks was also proposed 55 . As compared to some of the other trajectorybased entropy estimation methods, these inference schemes do not involve the estimation of probability distributions over the phase-space. Rather they usually only involve means and variances of measured currents. Hence they are known to work better in higher dimensional systems 15 . In addition, it is proven that such an optimization problem gives the exact value of the entropy production rate in a steady state as well as the exact value of the thermodynamic force field in the phase space of the degrees of freedom we can measure, if short-time currents are used [52][53][54][55] . However, these methods have not yet been tested against experimental data to the best of our knowledge.
Here we test the short-time TUR scheme against the challenges posed by experimental setups involving colloidal particles in time-varying potentials with (possible) background flows. In order to benchmark the scheme, we first test it in a setup where the entropy production rate of the system can be analytically predicted for any set of parameters. For this setup, we test our predictions against both analytical results as well as another recently proposed numerical scheme, namely stochastic force inference (SFI) 19 . After this benchmarking exercise, we apply our scheme to a modified system for which the underlying model is both unknown and hard to estimate. Though there is no theoretical value to compare within this case, the short-time TUR's predictions are again in perfect agreement with that predicted by the SFI technique 19 . These results provide a motivation for modeling this system in terms of coupled Langevin equations with two free parameters. We demonstrate that such a model does indeed capture the experimental observations, hence demonstrating the usefulness of these schemes in modeling complex scenarios.

Results and discussions
Model. Our results apply to systems with continuous state-space but a finite-number of degrees of freedom, described by overdamped Langevin equations of the type Here μ ¼ 1; ; d is the number of degrees of freedom of the system and we use ⋅ to refer to the Ito convention. F μ (X) is a function of X, but not an explicit function of time t, ξ μ is dÀ dimensional white-in-time noise such that ξ μ ðtÞξ ν ðt 0 Þ ¼ δ μν δðt À t 0 Þ, where Á h i denotes averaging over the statistics of the noise. The corresponding Fokker-Planck equation for the probability distribution function P is given by: where the repeated indices are summed over. In the steady state ∂ t P = 0. The total rate of entropy production σ can be obtained as 11,39 , is called the thermodynamic force field 15 . Overdamped Langevin equations are excellent descriptions for colloidal particle systems. Even for systems where the Langevin equation is not known, the fact that such a description exists in principle is all that is needed in order to apply Eq. (3a) and obtain σ by determining the current and steady-state probability density directly from the time-series data 11,15 . Another approach is to first infer the terms in the Langevin equation, F μ and D 16,19 and use Eq. (3a) to obtain σ. These methods can be applied directly on data obtained from tracking the system or even by using tracking-free methods in image space 16 .
Short-time TUR approach. In this paper, we demonstrate an alternative method for the simultaneous determination of both the entropy production rate as well as the thermodynamic force field F μ from experimental data, using the recently introduced short-time thermodynamic inference relation [52][53][54] . Our method is built on an exact result obtained in [52][53][54] : where k B is the Boltzmann constant and J is a weighted scalar current constructed from the non-equilibrium stationary state as shown below. The notation 〈⋅〉 stands for an ensemble average. The current that maximizes the term within the square brackets is J ∝ ΔS tot . Here Δt is the short-time interval over which the mean and variance of the current is evaluated 52 . In this work, it also coincides with the sampling rate of the trajectory. As for the ordinary TUR 56 , our result too holds for any X that is even under time reversal. The equality in (4) holds only when X includes all degrees of freedom of the system. If not, then the RHS of (4) gives a lower bound. The proof presented for Eq. (4) in ref. 52 was based on exact results for non-trivial models. It was shown that Eq. (4) is a consequence of fluctuations of ΔS tot becoming Gaussian, in the Δt → 0 limit. Later in refs. 53 and 54 , Eq. (4) was rigorously proved for overdamped diffusive processes.
Let us now discretize X in time with time interval Δt: We use latin indices as superscripts for the discrete time labels and the Einstein summation convention is applied to the greek indices. For a given function d(X) we can now define a time-discretised scalar function constructed from the steady-state current, Any such current, when substituted in the expression inside the square brackets of Eq. (4) can be shown to give a lower bound σ L which is ≤σ. In addition, for a special value of d = d*, J ∝ ΔS tot and σ L = σ. The algorithm we use, which obtains this d* and σ through a maximization procedure is as follows: 1. We first obtain a time-series of experimental data: X k . 2. To be able to perform the maximization we use a set of basis functions ψ m (X), m = 1, …, M, in the space spanned by X such that where w m 2 R d and are the parameters to be optimized. We use two sets of basis functions: Gaussian and linear and generate all our results in both these bases, for comparison.
3. Maximize Eq. (4) to obtain σ. This maximization is done using a numerical optimizer: We start with an initial guess for w m , calculate the time-series J k , construct the function within the square brackets in (4) and then maximize over w m to obtain σ and also the set of values w Ã m such that d Ã ¼ ∑ M m¼1 w Ã m ψ m ðXÞ maximizes Eq. (4). The maximizing current J* is constructed from d* using Eq. (5) and in addition can be shown to be proportional to ΔS tot 52 .
Furthermore, the thermodynamic force is proportional to d* that maximizes (4) 52-54 , i.e., Hence, by solving an optimization problem, where the RHS of Eq. (4) is maximized in the space of all currents we can obtain σ as the optimal value as well as its conjugate thermodynamic force field, F ¼ c d Ã where the proportionality constant can be fixed by using Var(J*) = 2〈J*〉 at Δt → 0 52 as, c ¼ 2hJ Ã i Var ðJ Ã Þ . We note that, for any set of basis functions ψ m (X), m = 1, …, M which give an adequate representation of d(X), an analytic solution to the maximization problem is known 54 . This solution gives a deterministic estimate of σ as, where Further, the optimal coefficients can be directly computed without any optimization as Repeated indices are summed over as before. Numerically, this involves inversion of the matrix Ξ. On the one hand, if d; N, and M are not very large, this deterministic scheme is faster compared to a numerical optimization algorithm, and does not get stuck in local maxima. On the other hand, numerical optimization schemes can in principle simultaneously handle the optimization of parameters of the basis functions. This is discussed in some detail in ref. 53 . In addition, numerical optimization schemes have also been extended to systems driven in a time-dependent manner 57 , where it is as yet unclear how the deterministic scheme will perform. In this work, we implement the numerical optimization scheme using a particle-swarm optimizer. We provide a brief introduction to the algorithm in the "Methods" section. We note that refs. 53,54 have already demonstrated the feasibility of the scheme described here with numerical data. Here we test this scheme instead on controlled experimental setups.
Colloidal particle in a stochastically shaken trap. To test the inference scheme we first apply it to an experimental problem for which the rate of entropy production is known from theory 58-61 -a colloidal particle in a stochastically shaken optical trap. This model was first experimentally studied in 62 . We study it again in order to understand the limitations posed by experimental setups for our inference scheme as well as test and benchmark our scheme for a system where the results are known.
We trap a polystyrene particle in an optical trap; further details of how the experiment is performed may be found in the "Methods" section. We modulate the position of the center of the trap λ(t) along a fixed direction x on the trapping plane perpendicular to the beam propagation (+z). The modulation is a Gaussian Ornstein-Uhlenbeck noise with zero mean and covariance λð0ÞλðsÞ ¼ Aτ 0 expðÀjsj=τ 0 Þ, i.e., where η is Gaussian, has zero-mean and is white-in-time. The correlation time τ 0 is held fixed for all our experiments. Note that Aτ 0 can be interpreted as an effective temperature 63 . The dynamics of the colloidal particle is well described by an overdamped Langevin equation, where K is the spring constant of the harmonic trap, γ is the drag coefficient, ξ is the thermal noise, D = k B T/γ is the diffusion coefficient of the particle, and T the temperature of the medium. The noise ξ is also Gaussian, zero-mean, and white-in-time and mutually independent from the noise η in Eq. (10). Equations (11) and (10) together define the model we call the Stochastic Sliding parabola. Starting from arbitrary initial conditions for x and λ, the system reaches a non-equilibrium steady state, with the probability distribution function and current given respectively by 58 where the dimensionless parameters θ and δ are defined as, The rate of entropy production and the thermodynamic force field for this model are, In Fig. 1, we compare the above exact results to the outcome of the inference algorithm applied to numerically generated data for this model. Different sets of time-series data were generated by varying the noise amplitude ratio θ by varying A, keeping the other parameters fixed. In Fig. 1a, we show the trajectories of the system in the (x, λ) space for three different θ values. In Fig. 1b, we see that the inference algorithm predicts a value σ L which is lower than the true value σ in the beginning, but gets very close to the true value, after a relatively modest number of steps.
As we run the algorithm longer, σ L saturates to something very close to the actual value. The inference algorithm also simultaneously gives an optimal force field d*(x) which is very similar to the thermodynamic Force field F μ ðxÞ expected from theory (see Supplementary Fig. S1 in Supplementary Note 1). From Eq. (14), it is clear that σ increases linearly with θ or equivalently the parameter A. Figure 1c illustrates that the inference algorithm captures this behavior accurately. Since we are limited by the minimal resolution of the time series in probing the Δt → 0 limit of Eq. (4), the inferred value of entropy production is in general different from the exact value by an O½Δt term. For this model, we can also compute this correction analytically as (using expressions previously obtained in ref. 61  where σ Δt is the result one gets from Eq. (4) for a fixed value of Δt. Notice that the O½Δt correction increases with the value of θ. The inferred values of σ indeed lie between these two limits. Next, we tested the algorithm on experimentally generated data for the same model. In the experiments, we varied A ranging from 0.1 to 0.35 in units of 0:6 10 À6 À Á 2 m 2 s À1 (corresponds to θ varying from 0.22 to 0.77), while the other experimental parameters such as the trap stiffness, as well as the bath temperature, were assumed to be constant for the entire length of the experiment. In reality, however, the laser used to trap the particle is prone to power fluctuations, and there can also be minor changes in the bath temperature due to heating caused by the long exposure to the laser. For large values of θ, we also expect non-harmonic effects to be significant, due to the particle exploring the peripheral regions of the trap 64 . We comment in the following paragraph on the implications of these fluctuations for our results. An immediate consequence is however that the theoretically predicted values Eq. (12a) can only be used as a reference. We benchmark our results instead by comparing them with values obtained by the application of the stochastic force inference technique (SFI) scheme recently proposed in 19 , which gives an independent estimate of both σ as well as the force fields. Experiments for individual parameter sets were carried out for a duration of 100 s, with a sampling rate of 10 kHz for the particle position. Only about 2/3rd of the available experimental data was used and the remaining 1/3rd was discarded because of the presence of uncontrolled experimental errors in them. For the analyzed data, each of the 100 s long data sets were further divided into 12.5 s long patches, upon which the inference algorithm was then tested. In Fig. 2, we demonstrate the results of the analysis of the experimental data. The dark-blue dashed line corresponds to the theoretically predicted value, Eq. (12a), of the entropy production rate for the model given by Eqs. (10) and (11) with the given parameters. The blue line is the entropy production for a slightly modified model obtained by analyzing the data obtained from SFI and calculating the drift and diffusion terms from it (see Supplementary Note 2 for details). The region between the red dashed lines corresponds to the error bar set by the variation in model parameters (namely the drift and diffusion coefficients) in different experiments as quantified by the SFI analysis. The data points are the results of our inference algorithm as well as the SFI scheme. As is evident, our inference scheme predicts exactly the same or very similar values for σ as the SFI algorithm, for all values of A.
The prediction of the inference scheme and SFI matches also for the thermodynamic force (Fig. 3a, b). Namely, the optimal current d*(x), which we get as an outcome of our inference algorithm, also matchesF ðxÞ, which is F ðxÞ estimated from the trajectory data by means of the SFI technique 19 . We conclude that our inference algorithm infers the correct entropy production value, as well as the correct thermodynamic force field, for the experimental data, since we get the same results when using a completely independent and different technique.
A colloidal particle trapped near a microbubble. We now demonstrate how our scheme performs in estimating the rate of entropy generation for the case where the mechanical force on the colloidal particle is not known. For this purpose, we study a particle trapped in the vicinity of a microscopic bubble of size 20-22 μm. We have already used this experimental setup to study one or more microbubbles with colloidal particles moving in the liquid in a different context 65,66 . The microbubbles are nucleated on a liquid-glass interface. The surface is pre-coated by linear patterns of a MB-based soft oxometalate (SOM) material. We focus a laser beam on any region along this pattern, the SOM material gets intensely heated, and a microbubble forms. The top of the bubble is colder than its bottom where it is anchored to the interface. As the surface tension is a function of temperature, the variation of the surface tension along the surface of the bubble sets up a Marangoni stress, driving a flow along the surface of the bubble. Marangoni flow around freely floating bubbles under a temperature gradient have been studied both experimentally 67 and analytically 68 . The additional complexity here is the presence of the bottom surface on which the flow must satisfy no-slip boundary conditions. The flow around the bubble in this setup is not yet known in detail although an approximate description, valid if we are not too close to the bubble, has been developed 66 , as shown in Fig. 4. This flow drags the trapped colloidal particle and changes its steady-state probability distribution (Fig. 5a, b). Since the flow streamlines are directed towards the bubble, we expect that these will confine the trapped particle more than the case without the bubble. This is indeed the case as we show later. We expect that the underlying description of the particle is still an overdamped Langevin equation, including a flow velocity field u(x). However, the quantification of this flow field is rather difficult, even numerically, as argued above. As a result, we have a system where the details of the microscopic description and forces are unknown. Our inference scheme, on the other hand, is easily applicable even in this context.
At the level of the non-equilibrium trajectories of the system, we see that there is a qualitative difference from the case without the bubble. First, we see that the particle is more confined in the trap along the x direction, when there is a bubble in the vicinity (see Fig. 5b) as mentioned earlier. This confinement is caused both by the flow towards the bubble (as shown in Fig. 4), which gets balanced at the confined position by the opposing force of Fig. 2 The short-time inference scheme tested on experimental data. Test of our inference algorithm on experimental data for different values of the parameter A (or θ) where A is the amplitude of the noise in the Ornstein-Uhlenbeck process defined in Eq. (10) and θ is defined in Eq. (13). The dark-blue dashed line corresponds to the theoretical value given by Eq. (14). The squares and triangles corresponds to the entropy production rate σ estimated from the experimental data using our thermodynamic uncertainty relation (TUR)-based inference scheme (Eq. (4)) with a Gaussian basis and a linear basis, and using Δt = 0.1 ms. The error bars correspond to averages over eight independent realizations of duration 12.5 s. The circles correspond to σ estimated using the stochastic force inference scheme (SFI) 19 for the whole 100 s data set, and the error bars for these correspond to a self-consistent estimate of the inference error that the SFI provides 19 . The blue line corresponds to σ predicted by a model obtained from SFI (Supplementary Note 2), and the red dashed lines correspond to error bars for this SFI-based model. The parameters used in the experiment are as follows: the corner-frequency of the harmonic trap, f c = 135 ± 10 Hz, relaxation time of the Ornstein-Uhlenbeck process τ 0 = 0.0025 s, temperature of the aqueous medium T = 298 K, and the rate at which the trajectory is sampled: Δt = 0.0001 s. the confining potential, as well as the reduced fluctuations close to the bubble due to proximity effects 69 . Further statistical analyses also reveal weaker non-equilibrium currents (see Supplementary Note 3 and Supplementary Fig. S2 for details). Consistent with these observations, on applying the inference algorithm, we observe that the value of σ is substantially reduced in the presence of the bubble. The corresponding entropy production rates estimated are σ = 244.68 k B s −1 for the no-bubble case and σ = 7.66 k B s −1 for the case with the bubble. We also find that the thermodynamic force, estimated using the inference scheme, is significantly reduced along the x direction, and the force field is less tilted along that direction as compared to the case without the bubble, as shown in Fig. 5c.
To further analyze the effect of the bubble, we performed another experiment, where we trapped the particles at different distances from the bubble. As we go a distance d~1.5r (r is the radius of the microbubble) from the surface of the bubble, we see that the inferred value of σ gets closer to the value the system would have had in the absence of the bubble. This is demonstrated in Fig. 6.
The significance of the inferred value of σ has to be discussed in the light of these findings. In the case without the bubble, it is exactly the total heat dissipated to the environment as a consequence of maintaining the system in a non-equilibrium steady state (by shaking the trap). In the case with the bubble, however, this is not the case. We present a possible mathematical description of this situation as an overdamped Langevin equation with space-dependent diffusion and damping terms in an unknown flow field u(x). Since the trap constrains the particle motion on scales that are at least two orders of magnitude smaller than the distance to the bubble, u(x) is further assumed to be a constant u d at a distance d from the surface of the bubble. σ calculated from this model, reproduces the values we find from the experimental data, independent of u d , and purely as a consequence of the space-dependent diffusion and damping term, and the two fitting parameters a and b. This is demonstrated in Fig. 6. As we discuss in Supplementary Note 4, however, there is another component of the entropy production, related to the Fig. 3 The thermodynamic force fields obtained from the inference scheme. Thermodynamic force field obtained as the optimal field d*(x, λ) using Eq. (7) (shown in black) compared toF ðx; λÞ (shown in blue) which is the thermodynamic force field obtained using the stochastic force inference technique 19   work that the flow does against the confining potential 42,70 . This component, which does indeed depend on the value of u d , is not estimated by our inference scheme, due to the fact that u d is a field (corresponding to the velocities of the molecules of the thermal bath) which is odd under time reversal, for which the TUR does not hold 3,56,[71][72][73] . Hence, we expect that the values of σ we find close to the bubble are underestimates of the true value. We elaborate on this point in Supplementary Note 4.
Mathematical model. The colloidal system in the presence of the bubble and consequently the flow u d , can be simulated using the following equations: where, Here the parameters a and b can be tuned to match the experimental data. Particularly, 1/b stands for a characteristic length scale below which the flows created by the bubble are significant. When the distance of the trapped particle from the bubble is much greater than 1/b, we expect that the expressions will match the case without the bubble. Using a trial and error approach, we obtained the fit parameters as a = 282.743 and b = 1/3 μm −1 . We remark that, as we did for the case without the bubble, the SFI technique could be used to model this case as well, since it explicitly gives the drift and diffusion terms. These are however particularly susceptible to erroneous estimates of a conversion factor, which is needed to obtain the particle trajectory data in units of nm. We expand on this issue in the "Methods" section as well as in Supplementary Note 2. This error can be thought of as assigning wrong units to the affected phase-space coordinates. Since σ is a sum over all phase-space coordinates, its evaluation is not affected by such an error unlike other quantities such as forces, diffusion terms, and the thermodynamic force. Another way to understand this is to note that σ quantifies the irreversibility of the dynamics, which again is clearly not affected by a choice of units. The colloidal system in the presence of the bubble. a The microbubble-colloidal particle system. b System trajectories without (red) and with (green) the bubble in the neighborhood of the colloidal particle. We see that the colloidal particle is strongly confined in the presence of the bubble. c The thermodynamic force field computed as the optimal field d*(x) (Eq. (7)) without the bubble (red) and in the presence of the bubble (green). The corresponding entropy production rates estimated are σ = 244.68 k B s −1 for the no-bubble case and σ = 7.66 k B s −1 for the case with the bubble. The parameters used in the experiment are as follows: the corner-frequency of the harmonic trap, f c = 57 ± 3 Hz, relaxation time of the Ornstein-Uhlenbeck process τ 0 = 0.025 s, temperature of the aqueous medium T = 298 K, the rate at which the trajectory is sampled: Δt = 0.0001 s and the amplitude of the Ornstein-Uhlenbeck noise A ¼ 0:3 ð0:6 10 À6 Þ 2 m 2 s À1 . Hence our model for the setup with the bubble only tries to reproduce the value of σ as a function of distance.

Conclusion
In conclusion, we have experimentally tested a simple and effective method, based on the thermodynamic uncertainty relation [52][53][54] for inferring both the rate of entropy production σ and the corresponding thermodynamic force fields, in microscopic systems in non-equilibrium steady states. We have confirmed that an entirely independent method, SFI 19 , gives the same answers in all the situations we have studied, hence adding weight to the physical significance of our findings. We have also carried out an extensive investigation of the convergence properties of our code as several parameters or hyper parameters are varied, as well as a comparison with the SFI algorithm (see Supplementary Notes 5, which includes Supplementary Figs. S3-S5 and Supplementary Notes 6, which includes Figs. S6 and S7). Our short-time inference scheme does not need any model in order to be applicable. However, we can use our findings to come up with plausible models, which give the same σ values for a range of parameters, even in cases where modeling the system from the first principles is complicated. In this regard, it would also be interesting to perform a systematic study of different algorithmic schemes available to model a complex nonequilibrium systems, with a focus on the advantages and disadvantages when applied to experimental data.
Experimental systems that would be particularly interesting to study are molecular motors or other cellular processes. Recently, ref. 74 tried to quantify the activity of a cell by measuring the power spectral density of the fluctuations of the position of a phagocytosed micron-sized bead inside a cell. As it is possible to also trap such beads inside a cell with optical tweezers 74 , this too could be a very interesting system to study. Finally, in other recent work 57 , it has been demonstrated that inference schemes of this kind can also be made to work for non-stationary non-equilibrium states, further diversifying the scope of this class of techniques.

Methods Experiment
A single colloidal particle in a stochastically shaken trap. The experimental setup consists of a sample chamber placed on a motorized xyz-scanning microscope stage, which contains an aqueous dispersion of spherical polystyrene particles (Sigma-Aldrich) of radius r = 1.5 μm. The sample chamber consists of two standard glass cover-slips (of refractive index~1.52) on top of one another. The thickness of the chamber is kept at 100 μm by applying double-sided sticky tape in between the cover-slips. The aqueous immersion is made out of double distilled water at room temperature, which acts as a thermal bath. A single polystyrene particle is confined by an optical trap, which is created by tightly focusing a Gaussian laser beam of wavelength 1064 nm by means of a high-numerical-aperture oil-immersion objective (×100, NA = 1.3) in a standard inverted microscope (Olympus IX71). The trap is kept fixed at a height, h = 12 μm from the lower surface of the chamber in order to avoid spatial variation in the viscous drag due to the presence of the wall. The corner-frequency of the trap (f c ) is set to be 135 Hz. For the first set of experiments, the center of the trap is modulated (λ(t)) using an acousto-optic deflector, along a fixed direction x in the trapping plane, perpendicular to the beam propagation (+z). Thus, the modulation may be represented as a Gaussian Ornstein-Uhlenbeck noise with zero mean and covariance λðsÞλðtÞ ¼ Aτ 0 expðjt À sj=τ 0 Þ. The correlation time τ 0 is held fixed for all our experiments. We determine the barycenter (x, y) displacement of the trapped particle by recording its back-scattered intensity from a detection laser (wavelength 785 nm, copropagated with the trapping beam) in the back-focal plane interferometry configuration. The measurement is carried out using a balanced-detection system comprising of high-speed photo-diodes 75 , with sampling rate of 10 kHz and final spatial resolution of 10 nm. In all cases, the trap parameters including the conversion factors for the trajectory data were calibrated by fitting the probability distribution of the particle position in thermal equilibrium to the Boltzmann distribution PðxÞ ¼ ð2πDτÞ À1=2 expðÀx 2 =2DτÞ, where τ is the relaxation time in the trap, given by τ = 1/(2πf c ). We assume that the diffusion constant D has the room temperature value, D = 1.645 × 10 −13 m 2 s −1 . It is important to note that we assume that the trap parameters, as well as the conversion factor for the trajectory data, are unaffected when the Ornstein-Uhlenbeck modulation is turned on. However, in practice, the trap parameters can indeed be altered by small amounts over long durations of measurement (−100 s), primarily due to the power fluctuations of the trapping laser. Further, the probe particle also moves in the y and z directions in the trap, which we have not measured here. These factors led to issues that prevented us from producing an exact replica of the theoretical model in the experiment. However, we have taken into account these limitations in our analysis as detailed in Supplementary Note 2.
In the second set of experiments, i.e., for those with the microbubble, we employ a coverslip that is pre-coated by a polyoxometalate material 76,77 absorbing at 1064 nm as one of the surfaces of the sample chamber (typically bottom surface), and proceed to focus a second 1064 nm laser on the absorbing region. A microbubble is thus nucleated-the size of which is controlled by the power of the 1064 nm laser 76 . Typically, we employ bubbles of size between 20 and 22 μm. Note that the sample chamber also contains the aqueous immersion of polystyrene particles. We trap a polystyrene probe particle at different distances from the bubble surface, and modulate the trap center in a manner similar to the experiments without the bubble. The particle is trapped at a axial height corresponding to the bubble radius. The other experimental procedures remain identical to the first set of experiments. However, an important additional step here is the determination of the distance of the particle from the bubble surface. This we accomplish by using the pixels-todistance calibration provided in the image acquisition software for the camera attached to the microscope, which we verify by measuring the diameters of the polystyrene particles in the dispersion (the standard deviation of which is around 3% as specified by the manufacturer), and achieve very good consistency. Note that we obtain a 2-d cross-section of the bubble as is demonstrated in Fig. 5a, and are thus able to determine the surface-surface separation between the bubble and the particle with an accuracy of around 5%. During the experiment, we also ensure that the bubble diameter remains constant by adjusting the power of the nucleating laser-indeed the bubble diameter is seen to remain almost constant for the 100 s that we need to collect data for one run of the experiment.
As opposed to the previous setting, the particle is now trapped in a region where it experiences: (1) a temperature gradient created by the laser beam used to generate the microbubble, (2) the microscopic flow generated due to the bubble, which affects particle trajectory, and (3) Faxen-like corrections to the viscous drag coefficient of the encompassing fluid due to the proximity to a wall 78 -which is the bubble surface in this case. Now, since the trap parameters-particularly the conversion factor for the trajectory data-are determined for the equilibrium setting, it is clear that there could be significant deviations from those in the non-equilibrium configuration produced due to the presence of the bubble in close proximity of the trapped probe particle. It is also clear that the previous procedure for obtaining the trap parameters will not work in this case. This severely hampers any exact modeling and as a result we concentrate on getting only the value of σ and its variation as a function of the distance to the bubble, both of which are robust against the above errors.
Numerical algorithm. Our aim is to maximize a cost function C, which is a function of a set of parameters w. We use a particle-swarm optimization algorithm 79 to achieve this. A domain is chosen and N p particles are initialized in that domain. The kth particle follows Newtonian dynamics given by: Here ω k and V k are the position and velocity vector of the kth particle and A k is a stochastic function that depends on the position of all the particles. Different variants of this algorithm use different A. The simplest-the one that we use-is called the Original PSO. Let us first define the following:

•
The kth particle carries an additional vector P k which is equal to ω k for which the value of the function C as observed by the kth particle was maximum in its history.
• At any point of time let G denote the position of the particle in the whole swarm for which the function has the maximum value.
The function A is given by Here the Greek indices run over the dimension of space. W 1 and W 2 are two weights. The two terms in Eq. (19) push the particle in two different directions: one towards the point in history where the particle found the function to be a maxima and the other towards the point where the swarm finds the maximum value of the function at this point of time. These are multiplied by two random vectors U 1 and U 2 of dimension same as the dimension of space. Each of the components are independent, uniformly distributed (between zero and unity), random numbers. We keep track of the highest value of the function seen by the swarm and also the location of that point. There are two major advantages to this over standard gradient ascent algorithms: one, it does not require evaluation of the gradient of the function and two, it can be parallellized straightforwardly. All the numerical results reported in this paper are obtained using this algorithm. We implement this optimization scheme using open-sourced PYSWARM package in Python 80 , with a default choice for the hyper parameters.
Implementation of the algorithm. Here we describe how we applied this algorithm to numerical/experimental data. We generate numerical data using first order Euler integration of Eqs. (10) and (11) with a time step of Δt = 0.0001. In either case, we generate many copies of trajectories of length 12.5 s, and construct the cost function in Eq. (4) using Eq. (5) and Eq. (6). We have tried out two different choices of basis functions to construct d(X). The first one is a Gaussian basis in which we represent d(X) as, Making use of the spatial symmetry of the problem, we assume d(X) to be an antisymmetric function, with d( − X) =−d(X), and that reduces the dimensionality of the problem by a factor of 2. Here M is the number of Gaussian functions, and b i are the variance of the Gaussian in the x and λ direction. The centers of the Gaussian (x m , λ m ) are placed equally spaced in a rectangular region enclosing the data. Both M and b i are hyper parameters. We used M = 16 and b 2 x=λ ¼ fx=λg max =30. Secondly, we have also tried a linear basis (motivated by the prior knowledge of the linearity of the system) where we take In Supplementary Note 6 which includes Supplementary Figs. S6 and S7, we study the dependence of the output of the algorithm, for both basis functions, on Δt, length of the time-series data, N p as well as M.
Since we have used a finite amount of data to construct the cost function, it will be prone to statistical errors. Therefore, we independently maximize the cost function for different 12.5 s long data sets, and take their mean value as the optimized estimate of σ. We show the value of sigma inferred (σ L ) as a function of the number of steps in the optimization algorithm for different 12.5 s long data sets in Fig. 7.

Data availability
The data used to produce the results in this paper is available at https://doi.org/10.6084/ m9.figshare.14176664. Fig. 7 Inference for different copies of trajectories. The plot shows the inferred value of the entropy production rate σ L as a function of the number of steps in the optimization algorithm, using the linear basis. The data is shown for different data sets all spanning 12.5 s in total length, generated numerically for the same parameter choices as in Fig. 1b. The black dashed line corresponds to the theoretical estimate of the entropy production rate σ for this parameter choice.