Estimating time-dependent entropy production from non-equilibrium trajectories

The rate of entropy production provides a useful quantitative measure of a non-equilibrium system and estimating it directly from time-series data from experiments is highly desirable. Several approaches have been considered for stationary dynamics, some of which are based on a variational characterization of the entropy production rate. However, the issue of obtaining it in the case of non-stationary dynamics remains largely unexplored. Here, we solve this open problem by demonstrating that the variational approaches can be generalized to give the exact value of the entropy production rate even for non-stationary dynamics. On the basis of this result, we develop an efficient algorithm that estimates the entropy production rate continuously in time by using machine learning techniques and validate our numerical estimates using analytically tractable Langevin models in experimentally relevant parameter regimes. Our method only requires time-series data for the system of interest without any prior knowledge of the system’s parameters. While methods for estimating the entropy production rate of a stationary process are relatively well established, this is still a challenge in non-stationary conditions. Here, the authors propose a scheme to infer the exact value of the time-dependent entropy production rate as well as entropy production along with single realizations directly from trajectory data.

T he entropy production rate is an important quantitative measure of a non-equilibrium process and knowing its value is indicative of useful information about the system such as heat dissipated 1,2 , efficiency (if the non-equilibrium system in question is an engine [3][4][5] ) as well as free energy differences 6,7 (if the non-equilibrium process interpolates between two equilibrium states). In particular, the entropy production rate often characterizes the energy consumption of nonequilibrium systems 8 . It also provides useful information for systems with hidden degrees of freedom 9,10 , or interacting subsystems where information-theoretic quantities play a key role [11][12][13][14] .
The entropy production rate can be directly obtained from the system's phase-space trajectory if the underlying dynamical equations of the system are known [15][16][17][18] . This is not the case however for the vast majority of systems, such as biological systems [19][20][21] , and consequently, there has been a lot of interest in developing new methods for estimating the entropy production rate directly from trajectory data [22][23][24][25][26][27][28][29][30][31][32][33] . Some of these techniques involve the estimation of the probability distribution and currents over the phase-space 22,26 , which requires huge amounts of data. Some other techniques are invasive and require perturbing the system 1,2 , which may not always be easy to implement.
An alternative strategy is to set lower bounds on the entropy production rate [34][35][36][37][38] by measuring experimentally accessible quantities. One class of these bounds, for example, those based on the thermodynamic uncertainty relation (TUR) [38][39][40][41][42] , 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 [26][27][28][29] . Recently a similar variational scheme using neural networks was also proposed 30 . As compared to other trajectory-based entropy estimation methods, these inference schemes do not involve the estimation of any kind of empirical distributions over the phase-space and are hence known to work better in higher dimensional systems 26 . In addition, it is proven that such an optimization problem gives the exact value of the entropy production rate in a stationary state if short-time currents are used [27][28][29][30] . The short-time TUR has also been experimentally tested in colloidal particle systems recently 43 . However, whether these existing schemes work well for non-stationary states has not been explored as yet.
Non-stationary dynamics ubiquitously appear in biological phenomena such as in adaptive responses to environmental change 44 and spontaneous oscillations 45 , all of which are inevitably accompanied by energy dissipation. However, for a nonstationary system, it has only been possible to place bounds on the time-dependent entropy produced during a finite time interval under specific 46,47 or more general 48 conditions. In addition, there is no guarantee that these bounds can be saturated by any quantity related to the entropy production of the system. Hence there is no established scheme that has been proven to work for obtaining the exact entropy production rate under timedependent conditions.
Here, we address this problem by proposing a class of variational inference schemes that can give the exact value of the timedependent entropy production rate under non-stationary conditions as well as the entropy production along with single realizations. These schemes, which can be directly implemented on time-series data obtained from experiments, involve maximization over an objective function that consists of a single projected current determined from the data. We demonstrate that this objective function can either be of the form dictated by the recently proposed short-time TUR [27][28][29] or the form recently suggested in 30 , or a variation of these. The collection of these schemes works for both diffusive systems described by overdamped Langevin equations as well as finite-state-space systems described by master equations and work for both transients as well as stationary states.
We implement these variational schemes by means of an efficient algorithm that estimates the entropy production continuously in time by modeling the time-dependent projection coefficients with a feedforward neural network and by carrying out gradient ascent using machine learning techniques. This algorithm can in principle be directly used on real experimental data. As a proof of concept, here we consider time-series data generated by two models; one of a colloidal particle in a timevarying trap and the other of a biological model that describes biochemical reactions affected by a time-dependent input signal, for both of which we can obtain exact solutions for the timedependent entropy production rate as well as the entropy production along single trajectories. We then demonstrate that our proposed scheme indeed works by comparing the numerical implementation to our theoretical predictions (see Fig. 1).

Results
Short-time variational representations of the entropy production rate. The central results we obtain, summarized in Fig. 1, are applicable to experimental data from any non-equilibrium system, at least in principle, described by an overdamped Langevin equation or a Markov jump process even without knowing any details of the equations involved. Here, we use the model of a generic overdamped Langevin dynamics in d-dimensions in order to introduce the notations. We consider an equation of the form: where A(x, t) is the drift vector, and B(x, t) is a d × d matrix, and η(t) represents a Gaussian white noise satisfying hη i ðtÞη j ðt 0 Þi ¼ δ ij δðt À t 0 Þ. Note that we adopt the Ito-convention for the multiplicative noise. The corresponding Fokker-Planck equation satisfied by the probability density p(x, t) reads ∂ t pðx; tÞ ¼ À∇ jðx; tÞ; ð2Þ where D is the diffusion matrix defined by and j(x, t) is the probability current. Equations of the form Eq. (2) can, for example, be used to describe the motion of colloidal particles in optical traps [49][50][51][52] . In some of these cases, the Fokker-Planck equation can also be solved exactly to obtain the instantaneous probability density p(x, t). Whenever j(x, t) ≠ 0, the system is out of equilibrium. How far the system is from equilibrium can be quantified using the average rate of the entropy production at a given instant σ(t), which can be formally obtained as the integral 53 where F(x, t) is the thermodynamic force defined as Note that the Boltzmann's constant is set to unity k B = 1 throughout this paper. Further, the entropy production along a stochastic trajectory denoted as S[x( ⋅ ), t] can be obtained as the integral of the single-step entropy production where ∘ denotes the Stratonovich product. This quantity is related ARTICLE COMMUNICATIONS PHYSICS | https://doi.org/10.1038/s42005-021-00787-x to the average entropy production rate as σ(t) = 〈dS(t)/dt〉, where 〈 ⋯ 〉 denotes the ensemble average. Similar expressions can be obtained for any Markov jump processes if the underlying dynamical equations are specified 17 .
In the following, we discuss two variational representations that can estimate σ(t), F(x, t), and S[x( ⋅ ), t] in non-stationary systems, without requiring prior knowledge of the dynamical equation. We also construct a third simpler variant and comment on the pros and cons of these different representations for inference.
TUR representation. The first method is based on the TUR 26,38-42 , which provides a lower bound for the entropy production rate in terms of the first two cumulants of nonequilibrium current fluctuations directly measured from the trajectory. It was shown recently that the TUR provides not only a bound, but even an exact estimate of the entropy production rate for stationary overdamped Langevin dynamics by taking the short-time limit of the current [27][28][29] . Crucially, the proof in Ref. 28 is also valid for non-stationary dynamics.
This gives a variational representation of the entropy production rate, given by the estimator where J d is the (single-step) generalized current given by The expectation and the variance are taken with respect to the joint probability density p(x(t), x(t + dt)). In the ideal short-time limit dt → 0, the estimator gives the exact value, i.e., σ TUR (t) = σ(t) holds 28 . The optimal current that maximizes the objective function is proportional to the entropy production along a trajectory, J Ã d ¼ cdS, and the corresponding coefficient field is d * (x) = cF(x, t), where the constant factor c can be removed by calculating 2hJ d i=VarðJ d Þ ¼ 1=c.
NEEP representation. The second variational scheme is the Neural Estimator for Entropy Production (NEEP) proposed in Ref. 30 . In this study, we define the estimator σ NEEP in the form of a variational representation of the entropy production rate as where the optimal current is the entropy production itself, Again, in the ideal short-time limit, σ NEEP (t) = σ(t) holds. Eq. (9) is a slight modification of the variational formula obtained in Ref. 30 ; we have added the third term so that the maximized expression itself gives the entropy production rate. Although it was derived for stationary states there, it can be shown that such an assumption is not necessary in the short-time limit. We provide proof of our formula using a dual representation of the Kullback-Leibler divergence [54][55][56] in Supplementary Note 2. The box on the right displays the steps in our inference scheme: we train the model function d(x, t|θ) with parameters θ to get the optimal values θ * , and use them for estimating the (single-step) entropy production Δ b S. b Estimated entropy production along a single trajectory. The dashed green line is the estimated entropy production, and the solid black line is the true entropy production calculated analytically. The estimation is conducted for the trajectory depicted in panel (a) after training the model function using 10 5 trajectories. The blue circles (variance-based estimator (Eq. (14)) and green triangles (simple dual representation (Eq. (10)) are the estimated entropy production rate using (c) 10 4 or (d) 10 5 trajectories, and the black line is the true value. The model function is trained with the simple dual representation in both cases. As is evident, the variance-based estimator reduces the statistical error significantly. In (c) and (d), the mean of ten independent trials are plotted for the estimated values, and the error bars correspond to the standard deviation divided by ffiffiffiffiffi 10 p . In contrast to the TUR representation, NEEP requires the convergence of exponential averages of current fluctuations, but it provides an exact estimate of the entropy production rate not only for diffusive Langevin dynamics but also for any Markov jump process. Since there are some differences in the estimation procedure for these cases 28,30 , we focus on Langevin dynamics in the following, while its use in Markov jump processes is discussed in Supplementary Note 2.
Simple dual representation. For Langevin dynamics, we also derive a new representation, named the simple dual representation σ Simple by simplifying he ÀJ d i in the NEEP estimator as Here, the expansion of he ÀJ d i in terms of the first two moments is exact only for Langevin dynamics and hence this representation cannot be used for Markov jump processes to obtain σ (however, as shown in 57 , the equivalence of the TUR objective function to the objective function in the above representation continues to hold in the long-time limit). The tightness of the simple dual and TUR bounds can be compared as follows: In Langevin dynamics, for any fixed choice of J d , where we used the inequality 2a 2 b ≥ 2a À b 2 for any a and b > 0. Since a tighter bound is advantageous for the estimation 56,58 , σ TUR would be more effective for estimating the entropy production rate for the Langevin case.
On the other hand, σ NEEP and σ Simple have an advantage over σ TUR in estimating the thermodynamic force F(x, t), since the optimal coefficient field is the thermodynamic force itself for these estimators. In contrast, σ TUR needs to cancel the constant factor c by calculating 2hJ d i=VarðJ d Þ ¼ 1=c, which can increase the statistical error due to the fluctuations of the single-step current (see Supplementary Note 2 for further discussions and numerical results). In the next section, we propose a continuoustime inference scheme that estimates in one shot, the timedependent thermodynamic force for the entire time range of interest. This results in an accurate estimate with less error than the fluctuations of the single-step current. σ NEEP and σ Simple are more effective for this purpose, since the correction of the constant factor c, whose expression is based on the single-step current, negates the benefit of the continuous-time inference for σ TUR . In Table 1, we provide a summary of the three variational representations.
We note that the variational representations are exact only when all the degrees of freedom are observed; otherwise, they give a lower bound on the entropy production rate. This can be understood as an additional constraint on the optimization space. For example, when the i-th variable is not observed, it is equivalent to dropping x i from the argument of d(x) and setting d i = 0. We also note that the variational representations are exact to order dt; in practice, we use a short but finite dt. The only variational representation which can give the exact value with any finite dt is σ NEEP , under the condition that the dynamics are stationary 30 .
An algorithm for non-stationary inference. The central idea of our inference scheme is depicted in Fig. 1a. Equations (8), (9), and (10) all give the exact value of σ(t) in principle in the Langevin case, but here, we elaborate on how we implement them in practice. We first prepare an ensemble of finite-length trajectories, which are sampled from a non-equilibrium and non-stationary dynamics with Δt as the sampling interval: Here, i represents the index of trajectories, N is the number of trajectories, and M is the number of transitions. The subscript (i) will be often omitted for simplicity. Then, we estimate the entropy production rate σ(t) using the ensemble of single transitions fx t ; x tþΔt g i at time t. σ(t) is obtained by finding the optimal current that maximizes the objective function which is itself estimated using the data. Hereafter, we use the hat symbol for quantities estimated from the data: for example, b σ Simple ðtÞ is the estimated objective function of the simple dual representation. We also use the notation b σðtÞ when the explanation is not dependent on the particular choice of the representation. The time interval for estimating b σðtÞ is set to be equal to the sampling interval Δt for simplicity, but they can be different in practice, i.e., transitions {x t , x t+nΔt } with some integer n≥1 can be used to estimate b σðtÞ for example. Concretely, we can model the coefficient field with a parametric function d(x|θ) and conduct the gradient ascent for the parameters θ. As will be explained, we use a feedforward neural network for the model function, where θ represents, for example, weights and biases associated with nodes in the neural network. In this study, we further optimize the coefficient field continuously in time, i.e., optimize a model function d(x, t|θ) which includes time t as an argument. The objective function to maximize is then given by The optimal model function d(x, t|θ * ) that maximizes the objective function is expected to approximate well the thermodynamic force F(x, t) (or c(t)F(x, t) if σ TUR is used) at least at jΔt(j = 0, 1, . . . ), and even at interpolating times if Δt is sufficiently small. Here, θ * denotes the set of optimal parameters obtained by the gradient ascent, and we often use d * to denote the optimal model function d(x, t|θ * ) hereafter. This continuous-time inference scheme is a generalization of the instantaneous-time inference scheme. Instead of optimizing a time-independent model function d(x|θ) in terms of b σðjΔtÞ with a fixed index j, the continuous-time scheme needs to perform only one optimization of the sum Eq. (13). This makes it much more data efficient in utilizing the synergy between ensembles of single transitions at different times. This also ensures that we can get the Table 1 A summary of the comparison between the different variational representations considered in this work. σ NEEP is the Neural Estimator for Entropy Production (NEEP) representation (Eq. (9)), σ Simple is the simple dual representation (Eq. (10)), and σ TUR is the thermodynamic uncertainty relation (TUR) based representation (Eq. (8)). They have different applicabilities to Markov jump processes and Langevin dynamics. The optimal coefficient field d * (x) that maximizes a variational representation is equivalent to or proportional to the thermodynamic force F(x, t). The TUR representation is the tightest as shown in Eq. (11).

Rep.
Markov jump Langevin Optimal field Tightness σ NEEP Yes Yes smooth change of the thermodynamic force, interpolating discrete-time transition data.
Variance-based estimator. In principle, all the three variational representations work as an estimator of the entropy production rate as well. However, as we detail in Supplementary Note 2, once we have obtained an estimate of the thermodynamic force d * ≃ F (taking into account the correction term for b σ TUR ) by training the model function, it is possible to use a variance-based estimator of the entropy production rate, which can considerably reduce the statistical error. This is due to the fact that d hJ d i fluctuates around hJ d i more than d VarðJ d Þ does around Var(J d ) for any choice of d, for small dt (see Supplementary Note 2 for the derivation). The above advantage in using the variance as an estimator, instead of the mean, would normally be masked by noise in the estimation of d * . However, if the coefficient field is trained by b σ Simple or b σ NEEP with the continuoustime inference scheme, remarkably, d * is obtained with an accuracy beyond the statistical error of d hJ d i since it takes the extra constraint of time continuity into account. This results in the error of d VarðJ d Ã Þ being smaller than that of d hJ d Ã i, because of the difference in how the leading-order terms of their statistical fluctuation scale with dt. We note that b σ TUR is not appropriate for this purpose, since in this case, d * should be multiplied by 2 d hJ d i= d VarðJ d Þ to obtain an estimate of the thermodynamic force, which increases the statistical error to the same level as d hJ d i. In numerical experiments, we mainly use b σ Simple for training the coefficient field to demonstrate the validity of this new representation and use the variance-based estimator for estimating the entropy production rate.
We adopt the data splitting scheme 28,30 for training the model function to avoid the underfitting and overfitting of the model function to trajectory data. Concretely, we use only half the number of trajectories for training the model function, while we use the other half for evaluating the model function and estimating the entropy production. In this scheme, the value of the objective function calculated with the latter half (we call it test value) quantifies the generalization capability of the trained model function. Thus, we can compare two model functions, and expect that the model function with the higher test value gives a better estimate. We denote the optimal parameters that maximize the test value during the gradient ascent as θ * . Hyperparameter values are obtained similarly. Further details, including a pseudo code, are provided in Supplementary Note 1.
Numerical results. We demonstrate the effectiveness of our inference scheme with the following two linear Langevin models: (i) a one-dimensional breathing parabola model, and (ii) a twodimensional adaptation model. In both models, non-stationary dynamics are repeatedly simulated with the same protocol, and a number of trajectories are sampled. We estimate the entropy production rate solely on the basis of the trajectories and compare the results with the analytical solutions (see Supplementary Note 3 for the analytical calculations). Here, these linear models are adopted only to facilitate comparison with analytical solutions, and there is no hindrance to applying our method to nonlinear systems as well 28 .
We first consider the breathing parabola model that describes a one-dimensional colloidal system in a harmonic-trap Vðx; tÞ ¼ κðtÞ 2 x 2 , where κ(t) is the time-dependent stiffness of the trap. This is a well-studied model in stochastic thermodynamics 49,50,59 and has been used to experimentally realize microscopic heat engines consisting of a single colloidal particle as the working substance 60,61 . The dynamics can be accurately described by the following overdamped Langevin equation: Here γ is the viscous drag, and η is Gaussian white noise. We consider the case that the system is initially in equilibrium and driven out of equilibrium as the potential changes with time.
Explicitly, we consider a protocol, κ(t) = γα/(1 + αt), where the parameters α, γ as well as the diffusion constant D are chosen such that they correspond to the experimental parameter set used in 60 (see Supplementary Note 3). In Fig. 1, we illustrate the central results of this paper for the breathing parabola model. We consider multiple realizations of the process of time duration τ obs as time-series data (Fig. 1a). The inference takes this as input and produces as output the entropy production at the level of an individual trajectory b SðtÞ for any single choice of realization (Fig. 1b), as well as the average entropy production rate b σðtÞ (Fig. 1c, d). Here, the entropy production along a single trajectory b SðtÞ is estimated by summing up the estimated single-step entropy production: while the true entropy production S(t) is calculated by summing up the true single-step entropy production: Note that their dependence on the realization x( ⋅ ) is omitted in this notation for simplicity. Specifically, we model the coefficient field d(x, t|θ) by a feedforward neural network, and conduct the stochastic gradient ascent using an ensemble of single transitions extracted from 10 4 or 10 5 trajectories (see Supplementary Note 1 for the details of the implementation) with Δt = 10 −2 s and τ obs = 1s. We note that, in recent experiments with colloidal systems, a few thousand of realizations of the trajectories have been realized with sampling intervals as small as Δt = 10 −6 s 62 , and trajectory lengths as long as many tens of seconds 60,61 .
A feedforward neural network is adopted because it is suitable for expressing the non-trivial functional form of the thermodynamic force F(x, t) 30,63 , and for continuous interpolation of discrete transition data 64 . In Fig. 1b, the entropy production is estimated along a single trajectory. We can confirm the good agreement with the analytical value. In Fig. 1c, d, the entropy production rate is estimated using 10 4 and 10 5 trajectories. In both cases, the simple dual representation is used to train the model function on half the number of trajectories. On the other half, we use both the simple dual representation as well as the variance-based estimator in Eq. (14) for the estimation, in order to compare their relative merits. We see, quite surprisingly, that the variance-based estimator performs better than the simple dual representation and has much less statistical error. Since the simple dual representation is essentially just a weighted sum of the mean and variance, this implies that the error in it is due to the noise in the mean, as also explained above (and in Supplementary Note 2).
Another advantage of our method is that it also spatially resolves the thermodynamic force F(x, t), which would be hard to compute otherwise. To demonstrate this point, we further analyze a two-dimensional model that has been used to study the adaptive behavior of living systems 21,44,65,66 . The model consists of the output activity a, the feedback controller m, and the input signal l, which we treat as a deterministic protocol. The dynamics of a and m are described by the following coupled Langevin equations: where η a and η m are independent Gaussian white noises, aðmðtÞ; lðtÞÞ is the stationary value of a given the instantaneous value of m and l, and a linear function aðmðtÞ; lðtÞÞ ¼ αmðtÞ À βlðtÞ is adopted in this study.
We consider dynamics after the switching of the input as described in Fig. 2a. For separation of time scales τ m ≫ τ a , the activity responds to the signal for a while before relaxing to a signal-independent value, which is called adaptation 44 . Adaptation plays an important role in living systems for maintaining their sensitivity and fitness in time-varying environments. Specifically, this model studies E. coli chemotaxis 21,44,65,66 as an example. In this case, the activity regulates the motion of E. coli to move in the direction of a higher concentration of input molecules by sensing the change in the concentration as described in Fig. 2a.
In this setup, the system is initially in a non-equilibrium stationary state (for t < 0), and the signal change at t = 0 drives the system to a different non-equilibrium stationary state. We show the results of the estimation of the entropy production rate  14)). The black line is the true entropy production rate. The labels d1, d2, and d3 are the time instances when we also estimate the thermodynamic force as shown in (c). c Analytical solutions for the thermodynamic force at the instances d1, d2, d3, and (d) estimates over 10 4 trajectories. Here the horizontal axis is the direction of a, the vertical axis is that of m, and an arrow representing the size of 100 is shown at the top of each figure for reference. The brighter the color, the larger the thermodynamic force as shown in the color bar. Note that in this particular case, the thermodynamic force becomes weaker as time evolves, and hence the size of the vectors reduces. The system parameters are set as follows: the time constants τ a = 0.02, τ m = 0.2, the coefficients of the mean activity function α = 2.7, β = 1, the strengths of the white noise Δ a = 0.005(t < 0), 0.5(t ≥ 0), Δ m = 0.005, and the inhibitory input l(t) = 0(t < 0), 0.01(t ≥ 0), which are taken from realistic parameters of E. coli chemotaxis 65,66 . The trajectories of length τ obs = 0.1 are generated with the time interval Δt = 10 −4 . The simple dual representation (Eq. (10)) is used for training the model function. and the thermodynamic force in Fig. 2b, c, respectively. Because of the perturbation at t = 0, the non-equilibrium properties change sharply at the beginning. Nonetheless, the model function d(x, t|θ * ) estimates the thermodynamic force well for the whole time interval (Fig. 2c), and thus the entropy production rate as well (Fig. 2b). In particular, we plot the result of a single trial in Fig. 2b, which means that the statistical error is negligible with only 10 4 trajectories. We note that the entropy production rate is orders of magnitude higher than that of the breathing parabola model. The results of Figs. 1 and 2 demonstrate the effectiveness of our method in estimating a wide range of entropy production values accurately. In the numerical experiments, we have used Δt = 10 −4 s. We note that sampling resolutions in the range Δt = 10 −6 s to 10 −3 s have been shown to be feasible in realistic biological experiments 67 . We also note that an order of 10 3 realizations are typical in DNA pulling experiments 68 .
The thermodynamic force in Fig. 2c has information about the spatial trend of the dynamics as well as the associated dissipation since it is proportional to the mean local velocity F(x, t) ∝ j(x, t)/ p(x, t) when the diffusion constant is homogeneous in space. At the beginning of the dynamics (t = 0), the state of the system tends to expand outside, reflecting the sudden increase of the noise intensity Δ a . Then, the stationary current around the distribution gradually emerges as the system relaxes to the new stationary state. Interestingly, the thermodynamic force aligns along the m-axis at t = 0.01, and thus the dynamics of a becomes dissipationless. The dissipation associated with the jumps of a tends to be small for the whole time interval, which might have some biological implications as discussed in Refs. 21,66 .
So far, we have shown that our inference scheme estimates the entropy production very well in ideal data sets. Next, we demonstrate the practical effectiveness of our algorithm by considering the dependence of the inference scheme on (i) the sampling interval, (ii) the number of trajectories, (iii) measurement noise, and (iv) time-synchronization error. The analysis is carried out in the adaptation model, for times t = 0 and t = 0.009, at which the degrees of non-stationarity are different. The results are summarized in Fig. 3. In most of the cases, we find that the estimation error defined by b σðtÞ À σðtÞ =σðtÞ is higher at t = 0 when the system is highly non-stationary. In Fig. 3a, b, we demonstrate the effect of the sampling interval Δt on the estimation. For both the t values, we find that the estimation error does not significantly depend on the sampling interval Δt in the range 10 −5 to 10 −3 , which demonstrates the robustness of our method against Δt.
In Fig. 3c, d, we consider the dependence of the estimated entropy production rate on N-the number of trajectories used for the estimation. We find that roughly 10 3 trajectories are required to get an estimate that is within 0.25 error of the true value for t = 0.009. On the other hand, we need at least 10 4 trajectories at t = 0 to get an estimate within the same accuracy. This is because the system is highly non-stationary at t = 0 and thus the benefit of the continuous-time inference decreases.
In Fig. 3e, f, the effect of measurement noise is studied. Here, the measurement noise is added to trajectory data as follows: where Λ is the strength of the noise, and η is a Gaussian white noise satisfying hη i a η j b i ¼ δ a;b δ i;j . The strength Λ is compared to Λ 0 = 0.03 which is around the standard deviation of the variable m in the stationary state at t > 0. We find that the estimate becomes lower in value as the strength Λ increases, while a larger time interval for the generalized current can mitigate this effect. This result can be explained by the fact that the measurement noise effectively increases the diffusion matrix, and its effect becomes small as Δt increases since the Langevin noise scales as / ffiffiffiffiffi Δt p while the contribution from the measurement noise is independent of Δt. Since the bias in d VarðJ d Þ is the major source of the estimation error, we expect that the use of a bias-corrected estimator 31,69 will reduce this error. Indeed, we do find that the bias-corrected estimator, star symbols in Fig. 3e, f, significantly reduces the estimation error (see Supplementary Note 1 for the details).
Finally, in Fig. 3g, h, the effect of synchronization error is studied. We introduce the synchronization error by starting the sampling of each trajectory att and regarding the sampled trajectories as the states at t = 0, Δt, 2Δt, . . . (actual time series is t ¼t;t þ Δt; :::). Here,t is a stochastic variable defined bỹ where uni(0, Π) returns the value x uniformly randomly from 0 < x < Π, the brackets are the floor function, and Δt 0 ¼ 10 À4 is used independent of Δt. The strength Π is compared to Π 0 which approximately satisfies σ(Π 0 ) ≈ σ(0)/2. We find that the estimate becomes an averaged value in the time direction, and the time interval dependence is small in this case.
In conclusion, we find that our inference scheme is robust to deviations from an ideal dataset for experimentally feasible parameter values and even steep rates of change of the entropy production over short-time intervals.

Conclusion
The main contribution of this work is the insight that variational schemes can be used to estimate the exact entropy production rate of a non-stationary system under arbitrary conditions, given the constraints of Markovianity. The different variational representations of the entropy production rate: σ NEEP , σ Simple , and σ TUR , as well as their close relation to each other, are clarified in terms of the range of applicability, the optimal coefficient field, and the tightness of the bound in each case, as summarized in Table 1.
Our second main contribution is the algorithm we develop to implement the variational schemes, by means of continuous-time inference, namely using the constraint that d * has to be continuous in time, to infer it in one shot for the full-time range of interest. In addition, we find that the variance-based estimator of the entropy production rate, performs significantly better than other estimators, in the case when our algorithm is optimized to take full advantage of the continuous-time inference. We expect that this property will be of practical use in estimating entropy production for non-stationary systems. The continuous-time inference is enabled by the representation ability of the neural network and can be implemented without any prior assumptions on the functional form of the thermodynamic force F(x, t). Our work shows that the neural network can effectively learn the field even if it is time-dependent, thus opening up possibilities for future applications to non-stationary systems.
Our studies regarding the practical effectiveness of our scheme when applied to data that might conceivably contain one of several sources of noise, indicate that these tools could also be applied to the study of biological 19 or active matter systems 70 . It will also be interesting to test whether these results can be used to infer new information from existing empirical data from molecular motors such as kinesin 71 or F 1 -ATPase 72,73 . The thermodynamics of cooling or warming up in classical systems 74 or the study of quantum systems being monitored by a sequence of measurements [75][76][77][78] are other promising areas to which these results can be applied.  20)). (b)(d)(f)(h) are the error analysis of (a)(c)(e)(g) between the estimated and the true entropy production rate (:¼ jb σðtÞ À σðtÞj=σðtÞ) at two time instances t = 0 and 0.009. The strength of the measurement noise Λ is shown in units of Λ 0 , the standard deviation of m in the stationary state at t > 0. The strength of the synchronization error Π is shown in units of Π 0 , the approximate time instance when the entropy production rate becomes half of its initial value. The black line is the true entropy production rate as used in Fig. 2. a, b The estimation is robust against the choice of the sampling interval. c, d The number of trajectories required for the convergence is large near the initial highly non-stationary time. e, f As for the strength of the measurement noise Λ increases, the estimate reduces because of the effective increase of the diffusion matrix. A larger time interval as well as correcting the bias (via Eq. (S4) of Supplementary Note 1) substantially mitigate this effect. g, h The estimate becomes an averaged value in the time direction. In contrast to (e)(f), the time interval dependence is small. For (a)-(h), the simple dual representation (Eq. (10)) is used for training the model function, and the variance-based estimator (Eq. (14)) is used for the estimation. The mean of ten independent trials are plotted, and the error bars correspond to the standard deviation divided by ffiffiffiffiffi 10 p . The system parameters are the same as those in Fig. 2

Code availability
Computer codes implementing our algorithm and interactive demo programs are available online at https://github.com/tsuboshun/LearnEntropy.