Entropy production factor derivation

Consider a system described by a set of M real, random variables obeying some possibly unknown dynamics. A specific trajectory of the system over a total time T is given by X = {Xⁱ(t)∣t ∈ [0, T]}. Given an ensemble of stationary trajectories, the average EPR, $$ \dot{S}$$, is bounded by^5,6,32

where we have set k_B = 1 throughout and D_KL denotes the Kullback–Leibler divergence which measures the distinguishability between two probability distributions. $$ P\left[\mathbf{X}\right]$$ and $$ P\left[\tilde{\mathbf{X}}\right]$$ are the steady-state probability distribution functionals of observing the path X(t) of length T and the probability of observing its reverse path, respectively. Therefore, the KL divergence in Eq. (1) measures the statistical irreversibility of a signal, and saturates the bound when X contains all relevant, non-equilibrium degrees of freedom.

We further bound the irreversibility itself by assuming the paths obey a Gaussian distribution. Writing the Fourier transform of Xⁱ(t) as xⁱ(ω), where ω is the temporal frequency, and writing the column vector $$ \mathbf{x}\left(\omega \right)={\left(x^1\left(\omega \right),x^2\left(\omega \right),\dots \right)}^T$$:

where x^† denotes the conjugate transpose of the vector x evaluated at the discrete frequencies ω_n = 2πnT⁻¹. C(ω_n) is the covariance matrix in Fourier space with elements $$ C^{ij}\left({\omega }_n\right)=⟨x^i\left({\omega }_n\right)x^j\left(-{\omega }_n\right)⟩T^{-1}$$, and Z is the partition function. The expression for $$ P\left[\tilde{\mathbf{x}}\right]$$ is identical but with C⁻¹(ω_n) → C⁻¹(−ω_n) [see Supplementary Note 1]. Combining Eq. (1) with Eq. (2) and taking T → ∞, we arrive at our main result:

This defines the EPF, $$ \mathcal{E}\left(\omega \right)$$, which measures time reversal symmetry breaking interactions between M ≥ 2 variables, while integrating $$ \mathcal{E}$$ gives $$ \dot{S}$$. $$ \mathcal{E}\left(\omega \right)=D_{\mathrm{KL}}\left(P\left[\mathbf{x}\left(\omega \right)\right]\mid \mid P\left[\tilde{\mathbf{x}}\left(\omega \right)\right]\right)$$ measures the Kullback–Leibler divergence between the joint distribution of M modes at a single frequency ω. While this quantity does not scale with trajectory length, the density of modes near a particular frequency is related to the total trajectory time by Δω = 2πT⁻¹. Since  ± ω modes must be complex conjugates of each other and an overall average phase is prohibited by time translation invariance, asymmetry between these distributions can only be captured by relative phase relationships, quantified by their correlation functions. $$ \mathcal{E}$$ is large when one variable tends to lead another in phase, implying a directed rotation between these variables in the time domain.

As mentioned above, P[x(ω)] describes the dynamics of a non-equilibrium steady state, and no reversal of external protocol is assumed. Further, in writing an expression for $$ P\left[\tilde{\mathbf{x}}\left(\omega \right)\right]$$, we assume that the observables are scalar, time-reversal symmetric quantities, such as the chemical concentrations we analyze below.

The Gaussian assumption we make here makes Eq. (3) exact only for systems obeying linear dynamics. Nevertheless, $$ \mathcal{E}$$ is still defined for non-linear systems, where the integrated $$ \mathcal{E}$$ lower bounds the true $$ \dot{S}$$. To see this, consider projecting complex dynamics onto Gaussian dynamics by choosing a data processing procedure which preserves two point correlations but which removes higher ones. This can be accomplished by multiplying every frequency by an independent random phase — a post-processing procedure which can be applied to individual trajectories. Post-convolution, the integrated EPF is equal to the KL divergence rate between forward and backwards rates. From the data processing inequality, the KL divergence rate of the true fields must be higher, so that the integrated EPF lower bounds the true entropy production rate [see Supplementary Note 2]. In addition to bounding the true $$ \dot{S}$$, we expect the integral of $$ \mathcal{E}$$ to be a good approximation for the wide class of systems where linearization is reasonable. Such Gaussian approximations are starting points in many field theories, with higher order interactions accounted for by adding anharmonic terms in the action of Eq. (2). While this is not our focus here, we expect these additional terms to systematically capture corrections to $$ \dot{S}$$ that do not appear in Eq. (3). As Cⁱⁱ(ω) = Cⁱⁱ(−ω), the only contributions to $$ \mathcal{E}$$ come from the cross-covariances between the random variables of interest. As such, this bound yields exactly 0 for a single variable even though higher order terms may contribute to $$ \dot{S}$$.

This formulation extends naturally to random fields. For M random fields in d spatial dimensions, $$ \mathit{\varphi }=\left\{{\varphi }^i\left(\mathbf{x},t\right)\mid t\in \left[0,T\right],\mathbf{x}\in {\mathbb{R}}^d\right\}$$, the EPR density, $$ \dot{s}\equiv \dot{S}/V$$ where V is the system volume, is:

where C^ij(q, ω) is the dynamic structure factor and $$ \mathcal{E}\left(\mathbf{q},\omega \right)$$ is now a function of wavevector q and frequency ω [see Supplementary Note 1].

Even without an explicit, analytic expression for the structure factor, C, we can estimate $$ \mathcal{E}$$ from data. To use Eq. (4), we consider data of N finite length trajectories of M variables over a time T in d spatial dimensions. Each dimension has a length L_i. We create an estimate of the covariance matrix, $$ \tilde{C}\left(\mathbf{q},\omega \right)$$, from time-series using standard methods [see Methods]. These measurements will inevitably contain noise that is not necessarily time-reversal symmetric, even for an equilibrium system. Noise due to thermal fluctuations and finite trajectory lengths in the estimate of $$ \tilde{C}$$ from a single experiment (N = 1) will systematically bias our estimated $$ \mathcal{E}$$ by $$ \mathrm{\Delta }\mathcal{E}=\frac{M\left(M-1\right)}{2}$$ at each frequency and will thereby introduce bias and variance in our measurement of $$ \dot{s}$$. We can simply remove the bias from our measured $$ \mathcal{E}$$, but to reduce the variance, we smooth $$ \tilde{C}$$ by component-wise convolution with a multivariate Gaussian of width $$ \mathit{\sigma }=\left({\sigma }_{q_1},\dots ,{\sigma }_{q_d},{\sigma }_{\omega }\right)$$ in frequency space, giving $$ \widehat{C}$$. This is equivalent to multiplying each component of the time domain $$ \tilde{C}\left(\mathbf{r},t\right)$$ by a Gaussian, cutting off the noisy tails in the real space covariance functions at large lag times. We then use $$ \widehat{C}$$ in Eq. (4) to create our final estimator for the EPF, $$ \widehat{\mathcal{E}}$$, and thereby the EPR, $$ \widehat{\dot{s}}$$. We calculate and remove the bias in $$ \widehat{\mathcal{E}}$$ and $$ \widehat{\dot{s}}$$ in all results below [see Methods]. Smoothing $$ \tilde{C}$$ with increasingly wide Gaussians in ω and q leads to a systematic decrease in $$ \widehat{\dot{s}}$$ due to reduced amplitudes in $$ \tilde{C}$$ (Supplementary Notes 3 and 4, Supplementary Fig. 1).

To illustrate the information contained in $$ \mathcal{E}$$, its numerical estimation, and the accuracy of $$ \widehat{\dot{s}}$$, we analyze simulations of coupled, one-dimensional Gaussian stochastic fields for which $$ \mathcal{E}$$ and $$ \dot{s}$$ can be calculated analytically. We then study simulations of the reaction-diffusion Brusselator, a prototypical model for non-linear biochemical oscillators, and use $$ \mathcal{E}$$ to study how irreversibility manifests at different time and length scales as the system undergoes a Hopf bifurcation³³.

Driven Gaussian fields

Consider two fields obeying Model A dynamics³¹ with non-equilibrium driving parametrized by α:

where ξ(x, t) is Gaussian white noise with variance $$ ⟨{\xi }^i\left(\mathbf{x},t\right){\xi }^j\left({\mathbf{x}}^{\prime },t^{\prime }\right)⟩={\delta }^{ij}\delta \left(\mathbf{x}-{\mathbf{x}}^{\prime }\right)\delta \left(t-t^{\prime }\right)$$, D is a relaxation constant, and $$ \delta \mathcal{F}/\delta \varphi $$ is the functional derivative with respect to ϕ of the free energy $$ \mathcal{F}$$ given by:so that the fields have units of ℓ^1/2 and r penalizes large amplitudes.

The EPR density, $$ \dot{s}$$, is calculated analytically in two ways. First, we solve Eq. (1) directly using the Onsager–Machlup functional for the path probability functional of $$ \mathit{\eta }\left(x,t\right)={\left(\varphi \left(x,t\right),\psi \left(x,t\right)\right)}^T$$^4,34. Second, the covariance matrices are calculated analytically, used to find $$ \mathcal{E}$$ through Eq. (4), and integrated to find $$ \dot{s}$$. Both cases give the same result for $$ \dot{s}$$. The result for both $$ \mathcal{E}$$ and $$ \dot{s}$$ are [see Supplementary Note 5]:

We see that $$ {\mathcal{E}}^{\mathrm{DGF}}\ge 0$$ and exhibits a peak at (q, ω) = (0, ω₀(0)), where $$ {\omega }_0\left(q\right)=\sqrt{{\left(D\left(r+q^2\right)\right)}^2+{\alpha }^2}$$, indicating that the system is driven at all length scales with a driving frequency of α, dampened by an effective spring constant Dr. In addition, it is clear that multiple combinations of α, r, and D can give the same value for $$ \dot{s}$$ while $$ \mathcal{E}$$ distinguishes between equally dissipative trajectories in the shape and location of its peaks. In this way, $$ \mathcal{E}$$ gives information about the form of the underlying dynamics not present in the total EPR. We note that $$ {\mathcal{E}}^{\mathrm{DGF}}$$ is also recovered using an appropriately modified version of the generalized Harada–Sasa Relation introduced in 34 [see Supplementary Note 6].

We perform simulations to assess how well $$ \mathcal{E}$$ can be extracted from time series data of fields [see methods for details]. The estimated $$ \widehat{\mathcal{E}}$$ shows excellent agreement with Eq. (7) (Fig. 1). Integrating $$ \widehat{\mathcal{E}}$$ gives $$ \widehat{\dot{s}}$$, which also shows good agreement with $$ {\dot{s}}^{\mathrm{DGF}}$$.

Our estimator gives exact results for the driven Gaussian fields because the true path probability functional for these fields is Gaussian. In contrast, the complex patterns seen in nature arise from systems obeying highly non-linear dynamics. For such dynamics, our Gaussian approximation is no longer exact but provides a lower bound on the total irreversibility. To investigate how irreversibility correlates with pattern formation, we study simulations of the Brusselator model for biochemical oscillations³⁵. We begin by describing the various dynamical phases of the equations of motion. Next, we calculate $$ \mathcal{E}$$ and $$ \dot{S}$$ for only the reactions before adding diffusion to study the synchronized oscillations that arise in the one-dimensional reaction-diffusion system.

Reaction-diffusion Brusselator

We use a reversible Brusselator model^30,35–37 with dynamics governed by the reaction equations:

where {A, B, C} are external chemical baths with fixed concentrations {a, b, c}, and all the reactions occur in a volume V (Fig. 2a). The system is in equilibrium when the external chemical baths and reaction rates obey $$ b{k}_2^{+}{k}_3^{+}=c{k}_2^{-}{k}_3^{-}$$. When this equality is violated, the system is driven away from equilibrium and exhibits cycles in the (X, Y) plane. Definingthe Brusselator is at equilibrium when Δμ  =  0 and is driven into a non-equilibrium steady state when Δμ  ≠  0. We vary b and c to change Δμ while keeping the product $$ \left(b{k}_2^{+}{k}_3^{+}\right)\left(c{k}_2^{-}{k}_3^{-}\right)=1$$, keeping the rate at which reactions occur constant for all Δμ³⁸.

As Δμ increases, the macroscopic version of Eq. (8) undergoes dynamical phase transitions. For all Δμ, there exists a steady state (X_ss, Y_ss), the stability of which is determined by the relaxation matrix, R [Supplementary Note 7]. The two eigenvalues of R, λ_±, divide the steady state into four classes³³:

$$ {\lambda }_{\pm }\in {\mathbb{R}}_{<0}\to $$ Stable attractor, no oscillations

$$ {\lambda }_{\pm }\in \mathbb{C},\mathrm{Re}\left[{\lambda }_{\pm }\right]<0\to $$ Stable focus

$$ {\lambda }_{\pm }\in \mathbb{C},\mathrm{Re}\left[{\lambda }_{\pm }\right]>0\to $$ Hopf Bifurcation, limit cycle

$$ {\lambda }_{\pm }\in {\mathbb{R}}_{>0}\to $$ Unstable repeller

The eigenvalues undergo these changes as Δμ changes, allowing us to consider Δμ as a bifurcation parameter. We define Δμ_HB as the value of Δμ where the macroscopic system undergoes the Hopf bifurcation.

Non-equilibrium steady states are traditionally characterized by their circulation in a phase space^39–43. One may then question how it is possible to detect non-equilibrium effects in the Brusselator when the system’s steady state is a stable attractor with no oscillatory component. While this is true for the macroscopic dynamics used to derive λ_±, we simulate a system with finite numbers of molecules subject to fluctuations. These stochastic fluctuations give rise to circulating dynamics, even when the deterministic dynamics do not³⁶. We see persistent circulation in the (X, Y) plane when $$ {\lambda }_{\pm }\in {\mathbb{R}}_{<0}$$, with the vorticity changing sign around Δμ  =  0 (Supplementary Fig. 2).

In order to assess the accuracy of our estimated EPR, $$ \widehat{\dot{S}}$$, we calculate an estimate of the true EPR, $$ {\dot{S}}_{\mathrm{true}}$$, for a simulation of Eq. (8) by calculating the exact entropy produced by each reaction that occurs in the trajectory⁴⁴, and then fitting a line to the cumulative sum (Supplementary Fig. 3, Methods). We find that $$ \widehat{\dot{S}}$$ significantly underestimates $$ {\dot{S}}_{\mathrm{true}}$$ (note the logged axes in Fig. 2c) due to the Brusselator’s hidden dynamics. In the Brusselator, information is lost because the observed trajectories are coarse-grained — they do not distinguish between reactions that take place forward through the second reaction or backwards through the third reaction in Eq. (8). These pathways would be distinguishable if trajectory of B and C were also observable. Our method relies purely on system dynamics to give $$ \widehat{\dot{S}}$$. Eq. (1) is true only if all microscopic details are captured by trajectories X. If X is already coarse-grained, multiple microscopic trajectories will be indistinguishable and Eq. (1) will underestimate the true entropy production rate due to the data processing inequality ^19,45,46.

In order to account for this, we recalculate $$ \dot{S}$$ by considering the rate at which a given transition can occur as the sum over all chemical reactions that give the same dynamics [see Methods]. For example, a transition from (X, Y) → (X − 1, Y + 1) can occur via reaction $$ k_2^{+}$$ or $$ k_3^{-}$$ in the Brusselator, each of which produces a different amount of entropy in general. Looking only in the (X, Y) plane, it is impossible to tell which reaction took place. When calculating the entropy produced by only the observable dynamics, the rate of making the transition (X, Y) → (X − 1, Y + 1) is $$ k_{\mathrm{f}}=k_2^{+}+k_3^{-}$$, while the rate of making the reverse transition is $$ k_{\mathrm{r}}=k_2^{-}+k_3^{+}$$, and the entropy produced is $$ \log \frac{k_{\mathrm{f}}}{k_{\mathrm{r}}}$$. This estimate of the EPR, which we name $$ {\dot{S}}_{\mathrm{blind}}$$, is a coarse-graining of $$ {\dot{S}}_{\mathrm{true}}$$, giving the relation $$ {\dot{S}}_{\mathrm{blind}}\le {\dot{S}}_{\mathrm{true}}$$⁴⁷. We find that $$ {\dot{S}}_{\mathrm{blind}}$$ shows excellent agreement with $$ \widehat{\dot{S}}$$, indicating that the Gaussian approximation provides a good estimate for the observable dynamics even when the system is highly non-linear.

To further benchmark our estimator, we calculate $$ \dot{S}$$ using two alternative methods, one based on the thermodynamic uncertainty relation (TUR)^7,48 and one based on measuring first passage times (MFPT)⁴⁹. The prior method measures a macroscopic current based on a weighted average of a system’s trajectory, j_d, and estimates the EPR using the TUR for diffusive dynamics, $$ \dot{S}\ge 2⟨j_{\mathbf{d}}^2⟩{\left({\tau }_{\mathrm{obs}}\mathrm{Var}\left[j_{\mathbf{d}}\right]\right)}^{-1}$$, where 〈〉 and Var[] denote an ensemble average and variance taken after an observation time τ_obs¹⁷. The latter method requires measuring the MFPT of an observable $$ \mathcal{O}$$ constructed from the system’s dynamics to reach a threshold that depends on a user-defined error tolerance. We choose $$ \mathcal{O}$$ and the threshold based on a drift-diffusion approximation for the winding number of the Brusselator (Supplementary Note 8). Similarly to $$ \widehat{\dot{S}}$$, both of these methods saturate to the true $$ \dot{S}$$ for systems obeying linear dynamics. As such, they also approximate $$ {\dot{S}}_{\mathrm{blind}}$$, but we find that they provide a looser bound than $$ \widehat{\dot{S}}$$ (Supplementary Fig. 4).

Prior to Δμ_HB, both $$ \widehat{\dot{S}}$$ and $$ {\dot{S}}_{\mathrm{blind}}$$ show a shift in their trends, but $$ {\dot{S}}_{\mathrm{true}}$$ does not. The smooth transition is due to the finite system size we employ, and gets sharper as a power law as the system gets larger (Fig. 3a). The power law exponent measured from $$ \widehat{\dot{S}}$$ is nearly linear, consistent with the Gaussian assumption. The exponent differs from that of $$ {\dot{S}}_{\mathrm{blind}}$$ because our Gaussian assumption breaks down at the high values of Δμ where the maximum slope occurs (Fig. 3b).

The Hopf bifurcation for the Brusselator is supercritical²³, meaning the limit cycle grows continuously from the fixed point when Δμ − Δμ_HB ≪ 1. Further from the transition point, the trajectory makes a discontinuous transition. At our resolution in Δμ, this discontinuous transition is what underlies the shift in $$ {\dot{S}}_{\mathrm{blind}}$$ of the Brusselator. This same transition is present in $$ {\dot{S}}_{\mathrm{true}}$$, but is difficult to detect numerically for reasons we explain here. In the deterministic limit, $$ {\dot{S}}_{\mathrm{true}}=\mathrm{\Delta }\mu \left(J^{\mathrm{F}}-J^{\mathrm{R}}\right)$$, where $$ J^{\mathrm{F}}=b⟨x⟩k_2^{+}$$ and $$ J^{\mathrm{R}}=c⟨y⟩k_2^{-}$$ are the forward and reverse fluxes for transforming a B molecule into a C molecule. 〈x〉 is a constant, but by numerically integrating the deterministic version for Eq. (8), we observe a discontinuity in 〈y〉 above the Hopf bifurcation. However, J^F ≫ J^R, obscuring the discontinuity in $$ {\dot{S}}_{\mathrm{true}}$$ (Fig. 3c). Upon coarse-graining, we have $$ {\dot{S}}_{\mathrm{blind}}=\mathrm{\Delta }\mu \left(J_{\mathrm{blind}}^{\mathrm{R}}-J_{\mathrm{blind}}^{\mathrm{F}}\right)$$, with $$ J_{\mathrm{blind}}^{\mathrm{F}}=b⟨x⟩k_2^{+}+{⟨x⟩}^3{k}_3^{-}$$ and $$ J^{\mathrm{R}}=c⟨y⟩k_2^{-}+{⟨x⟩}^2⟨y⟩k_3^{+}$$. These two terms are equal to each other for Δμ < Δμ_HB and diverge continuously when Δμ ⪆ Δμ_HB, followed by the relatively large discontinuity in $$ J_{\mathrm{blind}}^{\mathrm{R}}$$ (Fig. 3c, inset).

One gains further insight into the dynamics through the transition by studying $$ \widehat{\mathcal{E}}$$ (Fig. 2b). For Δμ < Δμ_HB, $$ \widehat{\mathcal{E}}$$ exhibits a single peak that increases in amplitude while decreasing in frequency as Δμ increases. Above Δμ_HB, the peak frequency makes a discontinuous jump, the magnitude of the peak grows rapidly, and additional peaks at integer multiples of the peak frequency appear due to the non-linear shape of the limit cycle attractor. These harmonics are expected for dynamics on a non-circular path. For Δμ  <  Δμ_HB, the magnitude of the peak is independent of system volume, while it gains a linear volume dependence in the limit cycle. The width of the peak is also maximized near the transition, reflecting a superposition of frequencies present in the trajectories (Supplementary Fig. 5).

To investigate how dynamical phase transitions manifest in the irreversibility of spatially extended systems, we simulate a reaction-diffusion Brusselator on a one-dimensional periodic lattice with L compartments, each with volume V, spaced a distance h apart. The full set of reactions are now

where d_j = D_j/h², and D_j is the diffusion constant of chemical species j = {X, Y}. Qualitatively different dynamics occur based on the ratio D_X/D_Y. D_X/D_Y ≪ 1 yields static Turing patterns^10,29. We focus on the D_X/D_Y ≫ 1 regime which exhibits dynamic, excitable waves. All values of {a_i, b_i, c_i} are kept constant in each compartment.

In the steady state, the reaction-diffusion Brusselator has the same dynamics as the well mixed Brusselator, and so it is not surprising that it’s EPR curve as a function of Δμ is similar (Supplementary Fig. 6). However, unlike the well-mixed system, the Hopf bifurcation signals the onset of qualitatively distinct dynamics in the reaction-diffusion system. Prior to the Hopf bifurcation, there are no coherent, spatial patterns in the system’s dynamics (Fig. 4a). Above the Hopf bifurcation, system-spanning waves begin to emerge that synchronize the oscillations across the system (Fig. 4b). Following standard methods^50,51, we define the synchronization order parameter, 0 ≤ r < 1, using

where θ_j(t) is the phase of the oscillator at position x_j and time t, M is the number of oscillators (here, the number of lattice sites in our simulation), and T is the temporal extent of the data [Methods]. ψ denotes the overall phase, and r is close to zero in the asynchronous phase and approaches one as the oscillators synchronize.

Below Δμ_HB, r is low and rapidly approaches one as the system approaches the macroscopic bifurcation point (Fig. 4c). Like $$ \dot{S}$$, this transition occurs more sharply and closer to Δμ_HB as the system size increases, approaching the discontinuous transition to the limit cycle behavior (Fig. 4c, inset)⁵². Throughout these changes, the system is driven further from equilibrium, as reflected in the increasing $$ \widehat{\dot{s}}$$ (Fig. 4d). The shift to collective behavior is not reflected in $$ \dot{s}$$ as it is almost identical to $$ \dot{S}$$ found for the well-mixed Brusselator (Supplementary Fig. 6). Instead, $$ \mathcal{E}$$ carries the signature of the dynamical phase transition. For Δμ < Δμ_HB, $$ \widehat{\mathcal{E}}$$ shows peaks at high wavenumbers, reflecting that irreversibility is occurring incoherently over short length scales. Above Δμ_HB, as the system shows synchronized oscillations, there is an abrupt shift in the peaks of $$ \widehat{\mathcal{E}}$$ to low q, indicating that this collective behavior carries the majority of the irreversibility (Fig. 5b,c). We also infer that the collective behavior is partially composed of traveling waves due to the streaks in $$ \widehat{\mathcal{E}}$$ (Fig. 5b). The slight offset in the transition occurs for high values of Δμ  <  Δμ_HB where small regions synchronize for short periods of time, but system wide oscillations are not observed (Supplementary Fig. 7). Furthermore, the transition moves closer to the macroscopic transition point with increased volume of the individual compartments (Supplementary Fig. 7).

Calculating $$ \mathcal{E}$$ from data

Estimate $$ \mathcal{E}$$ requires estimating frequency-space covariance functions, or cross spectral densities (CSDs). Considering a set of M discrete, real variables measured over time: {Xⁱ(t)}, where t = Δt, …, T, with T = NΔt, and i = 1, …, M indexes the variables, we estimate the CSD using the periodogram,

where $$ x^i\left(\omega \right)=\mathcal{F}\left\{X^i\left(t\right)-⟨X^i\left(t\right)⟩\right\}$$ are the Fourier transforms of the centered variables over the frequencies ω_n = 2πnT⁻¹ for $$ n=\left[-\frac{N}{2},\frac{N}{2}\right]$$.

The periodogram, is known to exhibit a systematic bias and considerable variance in estimating the true CSD. Both of these issues can be resolved by smoothing $$ {\tilde{C}}^{ij}$$ via convolution with a Gaussian with width σ. This is equivalent to multiplying $$ {\tilde{C}}^{ij}$$ in the time domain by a Gaussian of width σ⁻¹. We then define our smoothed CSD as

Once $$ \widehat{C}$$ is calculated, we then use the discrete version of Eq. (3) to estimate $$ \mathcal{E}$$. The extension to higher-dimensional data is done as follows: taking into account the spatial lattice on which the data is taken in Eq. (12), convolving the result with a multivariate Gaussian in Eq. (13), and finally estimate $$ \dot{s}$$ using the discrete version of Eq. (4). The choice of smoothing width, σ, should be guided by the maximum curvature seen in the structure factor, C^ij58.

Bias in $$ \widehat{\mathcal{E}}$$ and $$ \widehat{\dot{S}}$$

Our estimates of $$ \widehat{\mathcal{E}}$$ and $$ \widehat{\dot{S}}$$ are biased. The bias is found by calculating the expected value of $$ \widehat{\dot{S}}$$ for a system in equilibrium. To do this, we assume that the true covariance function is C^ij  =  δ^ij and measurement noise plus finite sampling time and rate gives rise to Gaussian noise in both the real and complex parts of $$ {\tilde{C}}^{ij}\left(\omega \right)$$, obeying the symmetries required for C^ij to be Hermitian. We only cite the results here and refer the reader to Supplementary Note 4 for a full derivation. The bias for random variables is

where M is the number of variables, $$ {\omega }^{\max }$$ is the maximum frequency available, σ is the width of the Gaussian used to smooth $$ \tilde{C}\left(\omega \right)$$, and T is the total time. The bias for random fields iswhere L_i is the length, $$ q_i^{\max }$$ is the maximum wavenumber, and $$ {\sigma }_{q_i}$$ is the width of the Gaussian used to smooth $$ \tilde{C}\left(\mathbf{q},\omega \right)$$ in the i^th spatial dimension.

Simulation details

To simulate the driven Gaussian fields, Eq. (5), we nondimensionalize the system of equations using a time scale τ  =  (Dr)⁻¹ and length scale λ  =  r^−1/2. We use an Euler–Maruyama algorithm to simulate the dynamics of the two fields on a periodic, one-dimensional lattice.

We simulate Eq. (8) using Gillespie’s algorithm⁵⁹ to create a stochastic trajectory through the (X, Y) phase plane with a well-mixed volume of V  =  100. We calculate the true $$ \dot{S}$$ of any specific trajectory z = {m_j∣j = 1, …, N} as follows. For each state $$ m^{\prime }$$, there exists a probability per unit time of transitioning to a new state m via a chemical reaction μ, denoted by $$ W_{m,m^{\prime }}^{\left(\mu \right)}$$. At steady state, the true entropy produced is⁴⁴

Note that ΔS_true is now itself a random variable that depends on the specific trajectory. We estimate $$ ⟨{\dot{S}}_{\mathrm{true}}⟩$$ by fitting a line to an ensemble average of ΔS_true (Supplementary Fig. 3), and compare that to $$ \widehat{\dot{S}}$$. We calculate $$ {\dot{S}}_{\mathrm{blind}}$$ by considering the "rate" at which a transition can occur as the sum over all the rates that give rise to the observed transition in (X, Y), i.e.,where $$ {\sum }_{\left\{{\mu }_j\mid m_{j-1}\to m_j\right\}}$$ denotes a sum over all reaction pathways μ that give rise to the transition m_j−1 → m_j. This procedure coarse-grains ΔS_true, giving ΔS_blind ≤ ΔS_true⁴⁷. ΔS_blind is the maximum entropy production that can be inferred by any method that observes trajectories in (X, Y), but which does not have access to the reaction pathways followed.

To simulate the reaction-diffusion Brusselator, Eq. (10), we take a compartment-based approach⁶⁰ where we treat each chemical species in each compartment as a separate species, and treat diffusion events as additional chemical reaction pathways. We nondimensionalize time by $$ \tau ={\left(k_1^{+}\right)}^{-1}$$ and use a Gillespie algorithm to simulate all reactions on a one-dimensional periodic lattice with L sites.

See Supplementary Tables 1–5 for all simulation parameters used in each figure.

Synchronization order parameter

The synchronization order parameter given in Eq. (11), r, is a function of the oscillator phase at every lattice site j at time t, θ_j(t). In order to calculate θ from our data, we measure the oscillator’s phase with respect to a trajectory’s mean position over time, namely

where the angle is taken over space and time. This phase is then used to calculate r as given in Eq. (11).