Constructing custom thermodynamics using deep learning

One of the most exciting applications of artificial intelligence is automated scientific discovery based on previously amassed data, coupled with restrictions provided by known physical principles, including symmetries and conservation laws. Such automated hypothesis creation and verification can assist scientists in studying complex phenomena, where traditional physical intuition may fail. Here we develop a platform based on a generalized Onsager principle to learn macroscopic dynamical descriptions of arbitrary stochastic dissipative systems directly from observations of their microscopic trajectories. Our method simultaneously constructs reduced thermodynamic coordinates and interprets the dynamics on these coordinates. We demonstrate its effectiveness by studying theoretically and validating experimentally the stretching of long polymer chains in an externally applied field. Specifically, we learn three interpretable thermodynamic coordinates and build a dynamical landscape of polymer stretching, including the identification of stable and transition states and the control of the stretching rate. Our general methodology can be used to address a wide range of scientific and technological applications.

a platform based on a generalized Onsager principle to learn macroscopic dynamical descriptions of arbitrary stochastic dissipative systems directly from observations of their microscopic trajectories.Our method simultaneously constructs reduced thermodynamic coordinates and interprets the dynamics on these coordinates.We demonstrate its effectiveness by studying theoretically and validating experimentally the stretching of long polymer chains in an externally applied field.Specifically, we learn three interpretable thermodynamic coordinates and build a dynamical landscape of polymer stretching, including the identification of stable and transition states and the control of the stretching rate.Our general methodology can be used to address a wide range of scientific and technological applications.
The modern scientific method adopts a universal approach that ensures stable and non-conflicting progression of our understanding of nature: new theories need to be hypothesized and tested on previously amassed data, be compatible with the basic scientific principles and be verifiable by experiments.Unfortunately, there is no general algorithmic recipe to do so in complex systems to facilitate discovery.Hence, up to now, only the most basic physical phenomena -often static, in equilibrium -are described by an intuitive set of equations.Many dynamic, non-equilibrium phenomena, which determine functionality in biology, soft-condensed matter, and chemistry, are instead described via very approximate, empirical laws.The advancement of artificial intelligence and machine learning gives rise to the possibility of a data-driven solution to this challenge [1,2].
In this paper, we develop Stochastic OnsagerNet, an AI platform which can discover an interpretable and closed thermodynamic description of an arbitrary stochastic dissipative dynamical system directly from observations of microscopic trajectories.There are essentially two types of approaches to understand and predict the behavior of dynamical processes from data -unstructured and structured.Unstructured approaches parametrizes dynamical equations by a generic set of building blocks, such as fixed polynomials [3], trainable feature maps or kernels [4,5], and determines the associated parameters that best fit the observations.Physical insights can be incorporated as regularizers in the fitting process [6].Their generality comes at a cost of long-time predictive accuracy, stability [7], and more importantly, interpretability.This is addressed by the class of structured approaches, where physical insights directly guide the design of model architectures.Our approach belongs to this latter category.Previous work in this direction include models based on Hamiltonian or symplectic dynamics [8,9], Poisson systems [10] and quasi-potentials [11].However, to date there lacks a general structured approach to model dissipative, non-equilibrium and noisy dynamics that often arise in soft matter, biophysics and other applications.Our methodology based on the classical Onsager principle [12,13] is tailored to such problems.
Macroscopic thermodynamic descriptions of physical systems are highly sought after for the insights they provide.A prototypical example is the ideal gas law as a macroscopic description of non-interacting gas systems.These guide the design of verification experiments and provides principled ways to manipulate macroscopic behavior.For a general complex dynamical system, however, constructing an intuitive thermodynamic description that enables subsequent analysis and control is a daunting task.Our approach addresses this challenge as follows.For a given microscopic dynamics, we learn a macroscopic thermodynamic description via the simultaneous construction of low-dimensional closure coordinates -ensured to be partially interpretable -and a time evolution law on these coordinates.Unlike general AI approaches, our platform intrinsically limits the search to physically relevant evolution laws.In particular, we ensure compatibility with existing scientific knowledge by constructing our neural network architecture based on a generalized Onsager principle.
We demonstrate our method by learning the stretching dynamics of polymer chains containing up to 900 degrees of freedom, condensing it into a thermodynamic description involving only three macroscopic coordinates that governs polymer stretching dynamics in both computational and experimental data.We build an energy landscape of the macroscopic evolution, revealing the presence of stable and transition states.This can be viewed as a dynamic equation of state.Mastering such an equation allows the design of verification computational experiments, including the interpretation of the thermodynamic coordinates and the control of the stretching rate of the polymers.We extend this further to conduct single-molecule DNA stretching experiments and show that our thermodynamic description can be used to distinguish fast and slow stretching polymers, much beyond current human-labelling capabilities.Furthermore, the predicted fluctuation correlations derived from the free energy landscape agree with experimental data.
Constructing low-dimensional physical models from high-dimensional dynamical data is an active area of research.Data-driven modeling of dynamical processes based on the Onsager's principle was proposed in Yu et al. [7] to study Rayleigh-Bénard convection.Chen et al. [14] combine encoder-decoders and manifold learning to construct latent dynamical models directly from video data, including those of reaction-diffusion processes and pendulum motions.Here, we make several advancements in terms of methodology and applications.First, unlike the deterministic models considered previously, we explicitly capture stochastic fluctuations -an important element of non-equilibrium processes at finite temperatures.In fact, the non-trivial heterogeneity of polymer stretching dynamics studied in this paper is directly caused by thermal fluctuations.Developing our method in the stochastic setting requires extensions of model reduction theory and training algorithms (see theoretical results and model implementation in Methods).Second and more importantly, we go beyond dimensionality reduction [7,14,15] and solve a closure problem: given a priori fixed macroscopic variables of interest (e.g.polymer extension), we construct both the closure coordinates sufficient to govern the evolution of these macroscopic variables, and the dynamics that describes this evolution.Compared with more flexible parametrizations of reduced dynamics [14], our approach inherently limit the search space to those satisfying a generalized Onsager principle, which sacrifices complete generality but facilitates physical interpretation of the closure coordinates and the dynamical landscape.

Results
We now describe our approach.The most complete description of a complex, multi-component system is the coordinates of all the components as a function of time t (X(t)).For the ideal gas it would be the positions and momenta of all molecules and for a magnetic system -the spin state of each atom.An alternative to this expensive microscopic modelling approach is a thermodynamic one, where the full description is replaced by some macroscopic coordinates (Z * (t), with dimensionality much smaller than X(t)).This can be the pressure of an ideal gas, or the magnetization of a magnetic system.The thermodynamic approach links these macroscopic coordinates to other macroscopic coordinates, or closure variables ( Ẑ(t)), and external parameters (volume and temperature for ideal gas, magnetic field and temperature for magnetic systems) via an equation of state.
We propose a generic approach of building such custom thermodynamics for an arbitrary stochastic, dissipative dynamical system from data.We are given macroscopic coordinates Z * whose dynamics we wish to model.For polymer dynamics this can be a single variable -the extension of the polymer chain, see Fig. 1.Then, we learn a set of closure variables Ẑ(t) and simultaneously an evolution law on the combined thermodynamic coordinates Z(t) = (Z * (t), Ẑ(t)) that enables scientific understanding, experimental verification, and control.The evolution equation is a generalization of the classical Onsager principle [12,13] that has been successfully applied to model a variety of non-equilibrium phenomena, including phase separation kinetics, gel dynamics and molecular modeling [16,17].It posits a time evolution law for a given set of coordinates Z(t), where the symmetric positive semi-definite matrix M models energy dissipation and V is a generalized potential.A limitation of the Onsager principle is its inability to capture dynamics far from equilibrium, or with substantial stochastic behavior.To this end, we propose an extension in the form of a generalized stochastic Onsager principle where M (•), W (•) are now functions of the reduced coordinates Z that output d × d matrices.M (•) is symmetric positive semi-definite to conform to stability requirements and Onsager's reciprocal relations [12], while W (•) is anti-symmetric and models conservative forces.The diffusion matrix σ(•) together with the white noise process Ḃ(t) models the thermal fluctuations in the system.Eq. ( 2) forms the basis of our dynamical model in reduced coordinates.We note that alternative generalizations of the Onsager principle have been proposed using large deviations theory [18], but their forms are more complex and hence less amenable to computations.It can be shown that our model has long-time stability through energy dissipation up to the order of thermal fluctuations (Thm.2) and the flexibility to represent many physical stochastic processes, including Langevin and generalized Poisson dynamics (see theoretical results in Methods).Our method departs from classical modelling paradigms, where the unknown equation parameters are few and can be fitted from few experiments.Instead, the unknowns here are functions M, W, V, σ.
We leverage machine learning and represent these functions as trainable deep neural networks, while preserving the required physical constraints (e.g.symmetric positive definiteness of M ).Simultaneously, we generate a set of the closure coordinates by another deep neural network, which combines approximation flexibility of residual networks and approximate feature orthogonality through principal component analysis.This is to be contrasted with generic coarse-graining methods based on volume averages [19], in that we seek a very small set of closure coordinates that are sufficient to describe the motion of the macroscopic states of interest.Our learning-based approach to discover hidden coordinates shares some similarity with recently proposed machine learningbased coarse-graining methods in molecular simulations [15,20] but here we work with a closure problem, thus we may end up with macroscopic dynamics of substantially lower dimensions.We perform end-to-end training of the combined architecture on large-scale microscopic trajectory data to simultaneously learn the reduced coordinates and their dynamics.The overall workflow for creating custom thermodynamics, which we call Stochastic OnsagerNet (S-OnsagerNet), is summarized in Fig. 1.Detailed network architectures and training algorithms are found in Methods.

Training and prediction of polymer stretching dynamics
We first demonstrate our approach by modelling the temporal evolution of polymer extension under elongational forces, which has long been of interest to the polymer physics community [21][22][23][24].Hallmark experimental [25,26] and computational studies [27][28][29] in elongational rheology of dilute polymers have examined the deformation of single DNA molecules in planar elongational flows and revealed the highly heterogeneous stretching dynamics among identical polymer chains.Due to the complex interactions within and stochastic nature of the system, it is challenging to identify macroscopic descriptors of the polymer chain (closure coordinates) and governing equations on these descriptors that are sufficient to determine the outcome of the stretching dynamics.Yet, such a thermodynamic description is essential for understanding the origins of unfolding heterogeneity and paves the way to make desired modifications to the unfolding dynamics.Thus, our data-driven method offers a promising alternative to achieve this goal.We simulate polymer chain stretching in a planar elongational flow.The polymer chain consists of 300 coarse-grained beads connected by rigid rods, resulting in 900 degrees of freedom if we ignore inertial effects (Fig. 2(A,B)).Snapshots of the shape of the chains under stretching conditions are shown in Fig. 2(C), revealing highly heterogeneous dynamics of the chain extension (Fig. 2(C,D)), defined as the projected chain length along the elongational axis of the flow.This is our macroscopic coordinate of interest Z * (t).Our aim is to model its stochastic evolution and understand the origin of its heterogeneity.We train the S-OnsagerNet on this dataset following the workflow in Fig. 1.The network architecture selection and training procedures are found in Methods.Our approach constructs 2 closure coordinates in addition to the chain extension Z * (t), leading to a 3 dimensional dynamical system -following Eq. ( 2) with learned functions M, W, V and σ -that governs the dynamics of stretching.We have empirically chosen the number of macroscopic coordinates: using more than 3 did not substantially improve predictive accuracy, whereas a 2-variable system has modelling limitations due to physical symmetry.The detailed selection procedure of the reduced coordinate dimension is discussed under polymer dynamics analysis in methods.
In Fig. 2(E-M), we test the trained S-OnsagerNet on three unseen, different and representative initial polymer configurations.The selected chains start with similar extension lengths, but subsequently stretch at vastly different rates.Fig. 2(E-G,I-K,M-P) show that the true statistics (black) can be accurately predicted (red).Moreover, the distributions of the time taken to reach a reference extension length are successfully captured (Fig. 2(H,I,M)).

Interpreting learned closure coordinates
Having shown that only two closure coordinates Ẑ = (Z 2 , Z 3 ) are required to characterize the stochastic evolution of the extension length Z * = Z 1 , it is natural to probe the meaning of these discovered coordinates.Here, we utilize an intrinsic property of neural networks -it represents the non-linear reduction functions X → Z as differentiable maps, as if we have learned their analytical forms.We compute via automatic differentiation the perturbations on a generic microscopic configuration X in the directions of ±∂Z 2 /∂X and ±∂Z 3 /∂X, corresponding to conformations with steepest changes in Z 2 and Z 3 respectively.The resulting conformational changes suggest physical interpretation of these coordinates.For example, from Fig. 3(A), we observe that perturbations in the direction of ∂Z 2 /∂X tend to change the end-to-end distance in the elongational axis, or distance between the first and the last bead in the polymer chain along the elongational axis.We confirm this hypothesis by visualizing the correlation of the end-to-end distance and the magnitude of Z 2 in Fig. 3(B,C).A similar analysis reveals the correlation between Z 3 and a degree of foldedness of the chain in the elongational axis of the flow (Fig. 3(D-F)).

Free energy landscape
The constructed potential V can be interpreted as a generalized free energy, allowing us to gain important insights into the dynamical landscape.The local minima of V represent stable or metastable states, while the saddle points correspond to transition states.The differentiable representation of V enables us to probe this landscape.Fig. 4 shows 2-dimensional projections of the 3dimensional free energy V (Z 1 , Z 2 , Z 3 ).We identify the critical points of V by solving ∇V (Z) = 0 using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method.We found two stable fixed points and two saddle points of interest marked in Fig. 4. Using a simultaneously trained PCA-ResNet decoder, we can reconstruct the macroscopic spatial coordinates of the polymer chain at the critical points to identify their physical origin.Up to reflection symmetry in the elongational axis, the stable points correspond to fully stretched states, whereas the saddle points refer to completely folded states.The origin of the heterogeneity in unfolding times is now clear: a rapidly stretching polymer is the one that avoids the saddle point and goes directly to a stable minimum, whereas a slowly stretching chain gets trapped around the stable manifold (attractive part) of the saddle point for a long time, before finally escaping through the unstable manifold (repulsive part) of the saddle towards the stable minimum (see Supplementary Video 1).We confirm this by overlaying a fast and a slow unfolding trajectory with the potential landscape in Fig. 4(C,F).Despite the similarity between the initial chain configurations, as demonstrated by the proximity between the initial points on the potential energy landscape, the chains exhibit different stretching behaviors that can be rationalized by the constructed potential.
Moreover, a Taylor expansion via automatic differentiation of the learned V (Z) captures the leading-order fluctuations near a stable stretched state.We denote by δV a typical energy fluctuation around the stretched state (proportional to temperature).Then, ignoring small terms we find where a 1 = 153.1, a 2 = 205.5, a 3 = 36.96,a 4 = 1.54 and δZ i = Z i − [Z stable,1 ] i is the fluctuations in the thermodynamic variables.Note that the coefficients a j implicitly depend on the strength of the flow.Eq. ( 3) is an effective equation of state, from which we observe the positive correlation of Z 1 and Z 2 .The physical interpretation is that near the stretched state, the chain extension and the end-to-end distance tend to change simultaneously (see Fig. 4(G-I) and polymer dynamics analysis in Methods).

Controlling polymer stretching dynamics
Further, understanding the laws of the custom thermodynamics of polymer chain folding allows us to interact with the dynamics by designing controls over the polymer environment to initiate desired changes in its behavior.To this end, we perform another Taylor expansion of V -this time near a saddle point Z saddle,2 corresponding to a folded state to give another local equation of state where b 1 = 102.96,b 2 = 31.13,b 3 = 24.16 and b 4 = 0.255.Eq. ( 4) suggests that to escape this saddle point leading to polymer unfolding, it is most effective to increase end-to-end distance (Z 2 ) while decreasing foldedness (Z 3 ) in a proportional way.This leads to a data-driven control protocol in Fig. 5(A,B).We choose the external elongational flow as the only control parameter (in real experiments it corresponds to switching on and off the flow of fluid or the electric field [26]).We start with a polymer configuration near the saddle point of the energy landscape, corresponding to a folded state.Without any intervention, the subsequent unfolding is expected to follow a slow trajectory, staying near the folded state for a long time.From our landscape analysis above, the most effective escape from the saddle point is along its unstable manifoldapproximately corresponding to increasing |δZ 2 − b 4 δZ 3 |.Thus, to speed up unfolding we can design the following control strategy: we turn off the external elongational flow, so the polymer drifts randomly under Brownian forces around the saddle point.We track its reduced coordinates, and once we observe sufficient alignment with the unstable manifold, we turn on the externally applied elongational flow.We observe in Fig. 5(C) that this simple control system speeds up the unfolding dynamics substantially.We can also increase the unfolding time by reversing this protocol (Fig. 5(D,E)).These control strategies based on the learned thermodynamic description have notable advantages over classical model-free control regimes (e.g.reinforcement learning), which may require large exploration times or small, finite state spaces [30].

Experimental validation
Remarkably, some qualitative predictions of the constructed thermodynamic description are confirmed by physical experiments.Not only can one show that the constructed dynamical landscape allows for fine-grained classification of stretching behavior on simulation data (see polymer dynamics analysis in Methods), we demonstrate in Fig. 6 that this applies directly to physical experiments.Here, we perform single-molecule experiments to observe the stretching trajectories of DNA molecules in a planar elongational field (see data preparation in Methods).We select two samples that initially appear similar (Fig. 6(E,G)), making it impossible to visually distinguish them in terms of stretching behavior.We then cast them into the learned thermodynamic coordinates Z, which when super-imposed on the free energy landscape reveals that the Z 2 coordinates differ subtly, leading to different predicted stretching statistics (Fig. 6(I)).This substantially improves upon human-level labelling, which can only occur much later in the dynamical evolution (Fig. 6(H,J)).Furthermore, we show in Fig. 6(J) that the effective equation of state Eq. ( 3) that captures the correlations of Z 1 and Z 2 around the stretched state also applies to experimental data from two sources, including the current experiments and previously available data [31].These results demonstrate the promise of the current approach in enabling physical understanding and control of real polymer dynamics.

Modelling spatial epidemics
To further demonstrate general applicability, we employ our method to derive macroscopic dynamics of spatial epidemics.The classical spatial SIR model [32] governs microscopic evolution of infective and susceptible individuals on a spatial domain (See Extended Data Figure 1(A-C)).Using these microscopic trajectories, we construct a thermodynamic model that accurately models the evolution of the spatial averages of infective and susceptible individuals (Extended Data Figure 1(D,E)) with an additional learned closure coordinate.Moreover, following the same approach before, we can interpret this coordinate as the spatial overlap of infective and susceptible individuals (Extended Data Figure 2), thereby rationalizing the dynamical landscape (Extended Data Figure 3), where this overlap determines the onset and outcome of disease spread.Details are found in Methods under spatial epidemics analysis.

Discussion
The potential applicability of our method goes beyond polymer and epidemic dynamics, and includes general complex dissipative processes such as protein folding [33], self-assembly [34,35] and glassy systems [36,37].Despite the importance of potential energy landscapes for functional properties material systems, the challenge in constructing them has limited current approaches to systems with small degrees of freedom, and/or requiring judicious selection of system descriptors based on expert knowledge [38].The method described in this work offers the potential of automating this process, creating pathways towards a multitude of opportunities for understanding and control over various complex systems and their scientific applications.There are many worthwhile future research directions to further improve the robustness and generality of the proposed method.Here, we inherently constrain the search space to macroscopic dynamics that conform to the generalized stochastic Onsager principle.Thus, it naturally has limited ability to model systems that may not readily admit such a description, such as chaotic systems.Moreover, the present model reduction and stability theory require the thermal noise to be small compared with the dissipative and conservative forces.While this is the case for the polymer dynamics studied here, our theory needs to be expanded in order to handle highly stochastic cases.In terms of training methodology, the current trial-and-error selection of the dimension of the closure variables can be made more systematic, e.g. by building on manifold learning approaches [14].Another potential improvement is the data sampling process.We observed in our numerical experiments that accurate construction of the dynamical landscape requires the trajectory data to sufficiently sample the regions of interest (stable and transition states).An adaptive or active learning algorithm [39,40] that couples data sampling and Constructing Custom Thermodynamics Using Deep Learning S-OnsagerNet training can be developed to improve upon the current random sampling strategy.On the scientific problem of polymer dynamics, we have only considered motion under a single stretching force.It is worthwhile to extend our study to varying stretching conditions in order to build a more comprehensive picture of polymer stretching.More broadly, one may apply our approach to learn macroscopic thermodynamics of other systems of scientific interest.

Theoretical results
In this section, we collect a number of theoretical results concerning the stochastic OnsagerNet (S-OnsagerNet) approach.We first show that if a high dimensional stochastic dynamical system satisfies the generalized stochastic Onsager principle (GSOP), then, any well-behaved reduction into a lower dimensional system will result in one that obeys approximately the GSOP introduced in Eq. ( 2) (Theorem 1).An immediate consequence is that our model reduction approach is theoretically justified for a wide variety of dissipative and conservative systems, including molecular dynamics [41], stochastic Hamiltonian systems [42], and stochastic Lotka-Volterra model [43].Next, we prove that dynamics described by the GSOP satisfy an energy dissipation law (see Theorem 2) and thus our machine learning approach produces stable dynamics at sufficiently low temperatures.
For convenience, we show that there are two equivalent forms of the GSOP.The formulation of the GSOP we use to construct our neural networks is Now, assuming Observe that Eq. ( 5) can be rewritten in the form A similar construction shows that we can also rewrite Eq. ( 6) in the form of Eq. ( 5) for any M , W , σ assuming similar invertibility conditions.Thus, they are in fact equivalent.While the numerical implementation is based on Eq. ( 5), the form Eq. ( 6) is also useful, and in the following we refer to both as GSOP.Now, we demonstrate the general applicability of the GSOP in the context of model reduction.We consider a microscopic (high dimensional) dynamics satisfying a GSOP of the form where X(t) ∈ R D , M 1 (•), W 1 (•) are symmetric positive semi-definite and antisymmetric matrix valued functions respectively, Σ 1 (•) is the D×p 1 dimensional diffusion matrix, B 1 is a p 1 -dimensional Brownian motion and ϵ 1 is a positive parameter related to temperature (e.g.ϵ 1 ∝ k B T with k B being the Boltzmann's constant and T the temperature).
We now show that many high dimensional systems of physical interest indeed satisfy a version of GSOP.We first consider the well-known Langevin dynamics, which has been used to model many stochastic dynamical systems e.g.molecular dynamics [44].
Example 1 The Langevin equation can be written in the form of Eq. ( 6), where the dot denotes a time derivative, ẋ is the velocity, ẍ is the acceleration, U (x) is the particle interaction potential, and so −∇U (x) is the potential force; γ 1 is the damping constant (units of reciprocal time), R(t) = Ḃ(t) is a delta-correlated stationary Gaussian process with zero-mean, satisfying If we set ẋ = v, the Langevin equation can be written as , the Langevin equation can be written in the form of the GSOP as follows: Another important class of dynamical systems are those described by Poisson brackets [45], whose stochastic extension encompasses many applications, including the stochastic Lotka-Volterra models and variants [43].In the following, we show that these dynamical systems can also be written in the form of the GSOP.
can be written in the form of Eq. ( 7), where is the Hamiltonian of the system, F is an arbitrary function depending on the system variables.The Poisson bracket {•, •} and the dissipation bracket [•, •] are defined as where M is symmetric positive semi-definite.
By the definition of {F, H} and taking W = 0 −In In 0 , we have Hence, Eq. ( 9) can be written as Next, let us consider the reduction of a microscopic dynamical system satisfying a GSOP (X(t)) into a macroscopic dynamical system (Z(t)).This is achieved by a differentiable reduction function ϕ : R D → R d such that Z(t) ≈ ϕ(X(t)).Moreover, we consider a differentiable reconstruction function ψ : R d → R D such that X(t) ≈ ψ(Z(t)).Our main result is that Z(t) also satisfies an approximate GSOP.In other words, the GSOP family is approximately invariant under dimensionality reduction, or coordinate transformation in general.
In the following, we adopt the notation We will also adopt the following technical assumptions.
1. Growth condition: there exists L > 0 such that, for all x ∈ R D and z ∈ R d , 2. Lipschitz condition: there exists L > 0 such that, for all x ∈ R D , and z ∈ R d , the function ( and ψ(z) satisfy the Lipchitz condition with constant L. 3. Approximate reconstruction: there exists ϵ 0 > 0 such that, sup x∈Ω |x − ψ(ϕ(x))| < ϵ 0 where Ω ⊂ R D is a domain such that X(t) ∈ Ω for all t ∈ [0, T ] almost surely.
Here, |•| is the Euclidean norm for a vector and Frobenius norm for a matrix.
Theorem 1 Let X(t) satisfy Eq. (7) and Z(t) satisfy Eq. (5) with Then, for each u ∈ C 2 (R d ) (twice continuously differentiable function) there exists a constant C > 0, independent of ϵ 0 and ϵ 1 , such that Proof Suppose Y (t) satisfy the following equation: we will prove By Itô's formula, Eq. ( 7) and Assump. 1, we obtain where ∇ 2 ϕ i is the Hessian matrix of ϕ i and Tr is the trace of a square matrix.Now, by definition, Y (t) satisfies the following SDE Subtracting equations ( 11) and ( 12), and integrating on [0, t], we get where By Cauchy-Schwarz inequality, and Itô isometry, we have By Assump. 1, there exists positive constants C 1 and C 2 such that We employ a similar above argument of the third term of Eq. ( 13) and get Here, we have used the Lipschitz conditions and the boundness of the first moment of f 1 (X) and g i (X), which is implied by the growth condition in Assump.1.This shows that ϕ(X(t)) and Y (t) are close in the mean-square sense.By Gronwall's inequality, we get Now, we employ a similar argument to show that u(ϕ(X(t))) and u(Y (t)) are close for any sufficiently smooth u.We apply Itô formula to u(ϕ(X(t))) and u(Y (t)) to get As before, by Assump. 1, there exist positive constants C 6 and C 7 such that Combining equations ( 14) and ( 15), we obtain By Jensen's inequality, we can get According to the Eq. ( 5) and ( 10), we can get Z(t) and Y (t) has the some distribution, i.e.
Finally, by triangle inequality, we have This completes the proof.□ This results demonstrate the validity of the GSOP as a dimensionality reduction method.In short, it says that if the microscopic dynamics satisfies a GSOP, then the macroscopic dynamics will also satisfy a GSOP approximately.Since a large amount of conservative and dissipative microscopic physical systems are shown to satisfy the GSOP, the S-OnsagerNet approach based on the GSOP is a principled model reduction ansatz for physical processes.
Next, we show the stability of a solution of the GSOP.More precisely, we prove in Theorem 2 below that the mean of the potential is non-increasing in t for sufficiently low temperatures (small |σ|).Consequently, the S-OnsagerNet produces dissipative dynamical systems that enjoy long term stability.

Theorem 2
The solution of Eq. (5) satisfies the dissipation law Here, where M 1/2 is the non-negative square-root of M .If we assume further that there exists a positive constant α such that Proof By Itô formula, we obtain where ∇V • M ∇V = ∥∇V ∥ 2 M and ∇V • W ∇V = 0 are used.Integrating Eq. ( 16) from 0 to t and taking expectation, we obtain Finally, according to the proposed condition, it is easy to arrive at thus the mean of the energy potential is non-increasing in t.This completes the proof.□ Theorem 2 shows that energy is being dissipated if the temperature of the ambient reservoir is sufficiently low.Accordingly, if the free energy V has compact level sets, then the dynamics at low temperatures will be confined on average to these compact sets and is thus stable.This contrasts with unstructured methods that may learn dynamics that are accurate for short times but induce instability at long times.We numerically verify this energy dissipation law in the learned polymer dynamics in Supplementary Figure 1.

Model implementation
We provide in this section the detailed implementation of S-OnsagerNet.
Model architecture.We begin with discussing the architecture design of the neural network approximators.Following the acquisition of the closure coordinates Z(t) = (Z * (t), Ẑ(t)), the S-OnsagerNet architecture implements equation (5).To ensure the symmetric positive definiteness of M (Z) and the anti-symmetry of W (Z), we use a neural network to approximate A(•) : R d → R d 2 with dimension d 2 .Then, we take the lower-triangular part as L 1 (Z) and the upper-triangular part as L 2 (Z).M (Z) and W (Z) are represented by where α is a positive constant and I is an identity matrix.
The energy function V (•) is lower bounded, so we use the following structure where U (Z) is a neural network with d-dimensional input and m-dimensional output, {γ ij } is a trainable matrix, and β is a positive parameter.
In the architecture used for the polymer dynamics application in this paper, we set α = 0.1 and utilize a neural network with 2 hidden layers with 20 neurons each and the tanh activation function to approximate M (z) and W (z). To parameterize the potential V (z), we decompose it into a sum of squares of the output layer (size m = 50) of 1 hidden layer neural network with 128 hidden neurons and the ReQUr activation function [7,46].This is to ensure that the potential satisfies the correct growth conditions as outlined in Assump.1.
For the diffusion matrix σ(z), since it has no symmetry constraints other than a growth condition, we use a fully connected neural network to approximate it.In our polymer dynamics application, we found empirically that a diagonal, z-independent diffusion matrix (corresponding to a linear network with zero weight and trainable diagonal bias) performed the best, but our algorithm can handle general architectures for σ(z).
Closure coordinate construction.We now provide details of the procedure to construct closure coordinates Ẑ(t) using the time series observation data of chain configuration coordinates at {t k } Nt k=1 , with 0 = t 0 < t 1 ... < t Nt = T .The available data are {X j } M j=1 with X j = {X(t i ) (j) } Nt i=1 ∈ M D×Nt , where M is the number of trajectories and X(t i ) (j) is the j th observation trajectory at t = t i .We obtained 610 observational trajectories, and for each trajectory, the number of time snapshots is 1001, i.e.M = 610 and N t = 1001.We reshape the observation data as X = [X 1 , X 2 , ..., X M ], where X ∈ M D×NtM .We re-center the data, such that the mean of the training data is 0 for each time snapshot, and set it as X.The covariance matrix of X is S = Cov(X).Denote its eigenvalues by Λ = diag(λ 1 , λ 2 , ..., λ D ) (arranged in non-increasing order) and corresponding eigenvectors as We use the following PCA-ResNet encoder to find the closure coordinates where and NN e (•) is a fully-connected neural network with input dimension D and output dimension d − 1.We can reconstruct the high-dimensional coordinates via the decoder where ) is another neural network with input dimension d − 1 and output dimension D. Note that without NN e , NN d , this amounts to a principal component analysis (PCA) based coordinate reduction.The combination of PCA and neural networks combines approximate feature orthogonality and approximation flexibility.
We construct the reconstruction loss function as |X − X| 2 .We set the reconstruction error of PCA as E pca = |X − P † d P d X| 2 .To make the reconstruction error near but less than the reconstruction error of PCA alone, we add the regularization term ReLU(log|X − X| 2 − log(E pca )) in the loss function, where is the reconstruction error of PCA and ReLU is the rectified linear unit, i.e.ReLU(u) = max(0, u).Thus, the combined reconstruction loss function is where ρ 1 is a regularization parameter.The regularization term penalizes the loss if we observe a reconstruction error that is higher than PCA.
Training algorithm based on maximum likelihood estimation.After constructing the structure of the drift term f (Z(t)) = −(M (Z(t) + W (Z(t)))∇V (Z(t)) and diffusion term σ(Z(t)), we consider how to construct the loss function to learn the stochastic dynamics.In deterministic dynamical systems, we can use the mean square error to learn f given the trajectory observation data.However, to deal with stochastic dynamics (in particular, learning the diffusion matrix σ), we have to devise more general methods based on maximum likelihood estimation.
Here, {ξ i } Nt−1 i=0 are independent random vectors following the standard normal distribution.
In the training data set, we have access to the microscopic coordinates X(t), from which we construct the reduced coordinates {(Z(t i ) (j) , Z(t i,j=0,1 via the (to be trained) reduction function ϕ.Given the neural networks f θ and σ θ (we use the subscript θ to denote all trainable parameters) constructed previously, the conditional probability is given by p where we use the short form σ θ = σ θ (Z(t i ) (j) ) and det denotes the determinant of a matrix.
Taking the logarithm of the above equation, we obtain As a result, we may obtain the loss function by maximizing the log-likelihood of the training data ) .
The total loss is then where ρ is a parameter to balance the accuracy of the learned dynamics and the error from reconstruction.In our computation, we first train the loss function (19) with ρ = 0.01.After some training steps, we fix the encoder part (17) and the decoder part (18) of the neural network, and train loss MLE (ρ = 0) to fine-tune accuracy of the stochastic dynamics.We use the Adam optimizer for training.The overall implementation, including the network architectures and loss computation, is shown in Supplementary Figure 2.

Data preparation
Simulation data.We used a Brownian dynamics approach to simulate linear, touching-bead chains as polymer chains in a planar elongational flow (left of Fig. 2).Each polymer chain consisted of N = 300 (D = 3N ) beads with diameter r at positions of bead i r i , connected by N − 1 rigid rods of length b = r = 10 nm.The governing stochastic differential equation was obtained by considering the following forces acting upon the system: excluded volume, constraint, Brownian and hydrodynamic.The excluded volume potential characterizes the short-range repulsions between beads and can be described by where µ = 35 pN has been demonstrated to result in a low frequency of chain crossings [47].The constraint force is given by where b i is the unit vector of bond i and T i is the tension in rod i that imposes constant bond length.The Brownian forces are random forces that satisfy the fluctuation-dissipation theorem, represented as where δ ij is the Kronecker delta, I is the identity matrix and ∆t is the simulation time step.By neglecting hydrodynamic interactions between the beads, the drag force on the ith bead is where ζ ≈ 3πηr is the drag coefficient of a single bead, η is the solvent viscosity and u(r i ) is the unperturbed solvent velocity.Due to the small mass of the beads, it is common to neglect inertial effects and set the sum of forces on the beads to be zero.Hence, the Langevin equation that describes the motion of each bead along the chain is We employed a predictor-corrector scheme [48] to determine the position of each bead at every time step.The enforcement of rigid rod constraints leads to a system of nonlinear equations for the rod tensions T i , which we solved for using Newton's method.
For each simulation run, the polymer chain was allowed to equilibrate for 10 4 τ d , with τ d = b 2 ζ/k B T being the characteristic rod diffusion time.During equilibration, the chain would adopt random configurations as governed by thermal fluctuations.At t = 0, the chain was subjected to a planar elongational flow of the form u(r i ) = ε(x − ŷ) • r i , where ε is the strain rate and x and ŷ are unit vectors parallel to the x and y axes respectively.The simulations were run until t = 10 4 τ d , using a time step of ∆t = 5 × 10 −4 τ d .To generate training data, we simulated 610 stretching trajectories.To test the predictions, we simulated 500 stretching trajectories each for three different initial chain configurations, which were deliberately selected from three vastly different trajectories.For each trajectory, we obtain the (x, y, z) positional coordinates r i of N beads every 10 τ d .Given the observation data, we can get the chain extension (Fig. 2).We note that time reported hereon is in units of τ d .Although the data set in this work is generated based on known equations, it should be highlighted that our machine learning approach for constructing the reduced dynamical model and all resulting consequences are independent of the microscopic model used in the simulations, with only the positional coordinates as inputs into the S-OnsagerNet.In other words, the approach used is purely data-driven and can therefore be generally applied to other non-equilibrium problems.
Experimental data.We provide the details of the experimental validation of the S-OnsagerNet results.
Experiments on electrokinetic stretching of DNA.To collect the experimental data leading to Fig. 6, we needed to create an automated single molecule stretching trap.We describe the essential features of the experiments, while the reader can find additional details about the trap and the material preparation methodology in Soh et al. [49].The polymer samples used for this study were T 4 phage double-stranded DNA (165.6kbp, Nippon Gene), chosen for high monodispersity and ready availability.The DNA was diluted in buffer solution, and fluorescently labelled (YOYO-1) to aid visualization.
The trap was based on a microfluidic cross-slot channel device with a wide central chamber, splitting into 40 µm wide channels in each of the four cardinal directions (Supplementary Figure 3).Each of the four channels was terminated by a macroscopic reservoir, in which we pipetted the DNA sample and inserted platinum wire electrodes.These electrodes were connected to a computercontrolled analog voltage source, such that the North and South electrodes were grounded, but East and West were positively biased V 0 = +30V .This electric quadrupole arrangement set up a potential well in the North-South direction and a potential hill in the East-West direction.The saddle point was located in the central chamber.
In aqueous solution, DNA is naturally negatively charged, and in the microfluidic device, molecules drifted electrokinetically from the North and South reservoirs towards East and West.At the saddle point, the East-West potential hill could be exploited to stretch a DNA molecule, but being an unstable equilibrium point, molecules would approach its location and slow down, but could not remain there for long observation times without external intervention.
To actively trap and observe single DNA molecules at the saddle point, the microfluidic electrokinetic device was placed on an inverted fluorescence microscope (Nikon Eclipse Ti2U) with a 60X oil-immersion objective (1.4NA).A live fluorescence image was captured by a sCMOS camera (Teledyne Photometrics Prime 95b) at 50ms intervals and sent to a desktop computer that analyzed the incoming images in real time.The image analysis included a fast cleanup step (detailed in the next section) that removed background noise and stray 'passer-by' molecules/fragments, followed by a step to calculate the intensity centroid of the target molecule and its projected length.The displacement x of the molecule's centroid from the saddle point was input to a proportional feedback loop that output a voltage tilt ∆V ∝ x, which was then superimposed on the East-West electrode biases, such that ∆V = V 0 ± Gx.Setting G = 2.2 V/µm confined the DNA centroid to within 1 µm of the saddle point even while a molecule stretched.This feedback process for tracking and trapping DNA molecules was automated through a custom LabVIEW program.
In a single stretch experiment, the platform was programmed to actively search for a molecule to trap as they flowed through the central chamber.Once a molecule was trapped, V 0 was set to zero temporarily to allow the molecule to relax into an un-stretched, equilibrium state over 10 s (chosen based on experience from preliminary experiments).After this relaxation period, recording began with images being streamed to the computer's solid-state drive.With every image, the associated centroid coordinates, projected length, voltages, and other parameters were logged in real time.V 0 was then reset to +30V, which re-established the East-West potential hill and stretched the DNA molecule.By monitoring the projected length history, the platform could recognize when the molecule was fully stretched, and in response stop the recording and release the molecule to escape naturally towards East or West.Using this protocol, we were able to capture a diverse set of stretching trajectories, including the dumbbell and folded conformations that were selected for analysis in Fig. 6.
Real-time image processing.A DNA molecule deforms over time, and the purpose of this section is to elucidate our method used to extract sequential image snaps that can be used to capture the DNA unfolding.In order to capture the exact location of the DNA molecule, we devised an algorithm to classify pixels according to the density of a pixel's immediate surrounding.The algorithm exploited the fact that the targeted molecule would very likely be centered in the tracking Region-Of-Interest (ROI), and this was used to discriminate against stray particles during tracking.To make this algorithm fast enough for real time, the algorithm operated only on the binarized version of the raw image, after application of a threshold.The result was a binary image containing the targeted DNA pixels at the exclusion of pixels belonging to noise speckles and stray particles.See Supplementary Figure 4 as an illustrated example of the real time image processing pipeline.
One would imagine that calculating the intensity centroid is a natural method for defining the molecule location and tracking it across image frames in time.However, doing so with the raw DNA images is problematic for the two following reasons: First: Background noise from the camera biases the centroid towards the center of the frame.This occurred even though we subtracted a pre-recorded background frame.The problematic noise in this case arises from pixel shot noise (electronic, photon) and pixel readout, and manifests as low-intensity sparkles in the image.One possible solution is to increase the excitation intensity, but this caused the DNA to photocleave much more readily, too often breaking into fragments before stretching out fully.
Second: While a trapped DNA molecule was being stretched, stray molecules continued to drift into the camera's field of view, which biased the centroid towards these stray molecules.This problem was partially solved by calculating the centroid from a smaller ROI that was just large enough to cover the fluctuating motion of a single molecule.This ROI was software-based and centered on the molecule centroid for tracking, i.e. implemented as a subset of the image array.The reduced ROI size avoided most of the stray particles, but it was still common for some to enter the tracking ROI during a trap-stretch cycle.
Data selection.Among the stretching trajectories collected, a few were selected for further analysis.Specifically, examples of molecules adopting the "dumbbell" and "folded" conformations during the stretching process with similar initial extension lengths were chosen in order to benchmark to human labeling.The trajectories of the conformations as the molecules were stretched were plotted in the reduced coordinates space (Fig. 6).When in the fully stretched stable state, the DNA molecules still undergo conformational fluctuations due to Brownian motion, which result in fluctuations in the reduced coordinate space.Images of molecules in this stable state were analyzed to obtain the fluctuations in Z 1 , Z 2 and Z 3 , as plotted in Fig. 6.The method to extract the learned thermodynamic coordinates Z 2 , Z 3 from experimental images are described next.
Extracting reduced coordinates from filtered images.Different from the simulated data described earlier, given a filtered image, it is more challenging to obtain in a robust way a sequence of ordered coordinates representing the location of each part of the DNA molecule.However, identifying the two end points is easier.Thus, we first find the coordinates of the two end points and also the centre of mass (weighted by intensity) of all illuminated pixels (see Supplementary Figure 5).Then, we estimate Z 2 , Z 3 by computing the end-toend distance and foldedness as defined in Fig. 3, which are shown to strongly correlate with Z 2 and Z 3 respectively.These computations only require the relative position vectors of the two end points with respect to the centre of mass (r 1 = (r 1,x , r 1,y , 0) and r N = (r N,x , r N,y , 0)).The third coordinate is set to 0 because the DNA molecule is confined to have limited motion in this Constructing Custom Thermodynamics Using Deep Learning direction.The center of mass r cm is computed by where Ī = m,n I m,n , u m,n is the spatial position of the (m, n)-pixel, I n,m is the intensity of the image at pixel (m, n).This allows for the computation of r 1 and r N (see Fig. [reffig:filter low]).According to the definition of end-toend distance and foldedness and their linear correlation with Z 2 , Z 3 (Fig. 3), we may compute the reduced coordinates as where L is the extension length scale of the experimental data, and the parameters C 1 , C 2 , C 3 and C 4 are scaling parameters to account for the change in molecule configuration properties (e.g.length scales) from simulation to experimental data.In particular • C 1 is obtained by dividing the average extension of unfolding simulated data by that of experimental data.• C 2 is the relationship between Z 2 and end-to-end distance (see Fig. 3 and is obtained by dividing the average value of Z 2 in simulated data by that of the end-to-end distance.• C 3 and C 4 is obtained by the relationship between the Z 3 and foldedness. We use least squares method to get C 3 and C 4 (see Fig. 3 (E)).
The method of extraction differs in the plot of Fig. 6 (J).Here, we only consider the stretched state, for which we can extract 300 coordinates (assuming the third coordinate z = 0) equally spaced along the stretched polymer.A scaling along the x-axis is performed so that the average full extension of the experimental polymer data matches that of the simulation data.These coordinates are then fed into the trained PCA-ResNet to extract the Z 2 and Z 3 values.

Polymer dynamics analysis
In this section, we detail the from the modelling of the polymer unfolding problem using the S-OnsagerNet.
Accurate Prediction of the statistics of unfolding.In our training process, we have 610 training trajectories and 110 testing trajectories.Note that although the dataset is generated by known yet complex microscopic equations, our approach does not require, nor rely on, the knowledge of these equations.The chain extension evolution of the training data with different initial chain are shown in Supplementary Figure 6(A) (black) and the test data are shown in Supplementary Figure 6(E) (black).After training the model, we predict the extension of the polymer (red).We can see the extension of the polymer can be predicted well.We also compute the mean (Supplementary Figure 6(B,F)), standard derivation (Supplementary Figure 6(C,G)) and the distribution of unfolding time (Supplementary Figure 6(D,H)) of the training and test results.We also compute the error of the mean (relative L 2 error), standard derivation, and the probability distribution of unfolding time of the training and test results in Supplementary Table 1.We observe that our model successfully captures the statistical behavior of a polymer stretching with only a 3-variable dynamical system.
Interpreting the learned closure coordinates.Supplementary Figure 7 shows the evolution of chain extension and the two learned closure coordinates with time for the training data.The trajectories are colored by unfolding time.Based on the unfolding times, we observe that chains with similar initial extensions (determined by the y-intercept, i.e., Z 1 at t = 0) can take vastly different times to stretch, hence it is not sufficient to consider only chain extension for the purposes of predicting dynamics.However, the successful prediction of the statistics of unfolding dynamics implies that Z 2 and Z 3 capture crucial information of the system.Thus, we seek to gain some physical understanding of the learned closure coordinates Z 2 and Z 3 .
Here, we provide details on the analysis of the closure coordinates (Z 2 , Z 3 ) that characterize the stochastic evolution of the extension length (Z 1 ), Recall that Z 2 , Z 3 are deterministic functions of the microscopic configuration X, and the functions are approximated by a trained PCA-ResNet.A useful property of we can exploit is differentiability of the neural network.We can ask: given a certain configuration X, what kind of small perturbations to X will most drastically increase or decrease the value of Z 2 and Z 3 ?In other words, we can consider perturbations in the directions of ±∂Z 2 /∂X and ±∂Z 3 /∂X respectively.We first analyze the closure coordinate Z 2 .From Fig. 3(A), we observe that perturbations in the direction of Z 2 tend to change the end-toend distance in the elongational axis (distance between the first and the last bead in the polymer chain along the elongational direction, i.e. |r N,x − r 1,x |), where r j,i is the i th coordinate of the j th bead in the chain (cyan: Z 2 = 0.453, blue: Z 2 = 0.194, and black: Z 2 = 0.323).We confirm this hypothesis by visualizing the correlation of the end-to-end distance and the magnitude of Z 2 in Fig. 3(B,C).Generally, we observe that as |Z 2 | decreases, the distance between the chain ends (marked by red points in the figure) decreases.Thus, we can interpret the first learned closure coordinate as an indicator of end-to-end distance.
We proceed with a similar analysis for the other closure coordinate Z 3 .Fig. 3(D) shows a given chain configuration and perturbations in the positive and negative directions of ∂Z3 ∂X (cyan: Z 3 = 7.903, blue: Z 3 = 1.595, and black: Z 3 = 4.749).Here, we observe that the end-to-end distance is largely unchanged, but the degree of foldedness of the chain in the elongational axis of the flow (x direction) appears to change.This leads us to hypothesize that the second learned coordinate represents a degree of foldedness with respect to the elongational flow.As a measure of the degree of foldedness of a chain, we compute |r 1,x + r N,x |.During data pre-processing, the chain is centered such that its center-of-mass is 0. Hence, if |r 1,x + r N,x | is small, the polymer is symmetric around 0 in the elongational x direction and tends to be in the elongated state.If |r 1,x + r N,x | is large, the polymer is likely to be in the folded state.We plot |r 1,x + r N,x | as a function of |Z 3 | for all configurations in the training data set in Fig. 3(E).The strong correlation between |r 1,x + r N,x | and |Z 3 | supports our interpretation that the second learned closure coordinate is an indicator of the degree of foldedness in the elongational direction.To demonstrate this, we plot in Fig. 3(F) visualizations of different chains with a range of |Z 3 | values, with the chain ends marked by red points.As |Z 3 | decreases, we observe the chains generally shift from the folded to the elongated state.We note that the degree of foldedness is sufficiently described solely by considering the projection of chain coordinates onto the elongational axis of the flow, as the flow is stable in the compressional axis and thus the degree of foldedness is primarily relevant to the unstable elongational axis that drives the unfolding process.
Advancing classification methods for polymer stretching.With the new understanding of the closure variables, we now consider how our results improve the current understanding of polymer chain dynamics.In the landmark experimental study of dilute polymer chains stretching under elongational flow, the molecules were categorized into seven different conformations and the dynamics of dominant conformations were analyzed [25].Specifically, it was found that chains in the "folded" conformation (which is one of the seven categories) took the longest time to stretch, while chains in the dumbbell conformation stretched relatively quickly (Supplementary Figure 8).
Our analysis shows that instead of a categorical labelling, it is perhaps more useful to characterize the stretching dynamics of a polymer by three numbers representing the generalized coordinates: Z 1 (extension length), Z 2 (related to end-to-end distance), and Z 3 (related to foldedness in the flow direction).Our results show that these are sufficient to predict the dynamics of the chain extension.We show that this characterization is largely consistent with previous categorical ones, but improves upon them in some intermediate cases.In Supplementary Figure 9, we plot the values of the low dimensional coordinates |Z 2 | and |Z 3 | of different chains with initial dumbbell and folded configurations at selected chain extension Z 1 values, colored by the predicted unfolding times.We observe segregation between the folded and dumbbell configurations in the Z 2 − Z 3 space, indicating that the qualitative differences between different conformations can be captured by our characterization.Generally, the folded chains take a longer time to stretch compared to the dumbbell chains.This is consistent with experimental and computational observations reported in the literature [25][26][27].However, we highlight that the region with high |Z 2 | and low |Z 3 | values encompasses a mix of folded and dumbbell chains with similar unfolding times.Therefore, while the broad categorization scheme allows for coarse discrimination of the stretching dynamics, the qualitative classification does not allow for finer predictions.Instead of classifying the stretching trajectories based on qualitative, human judgement of chain conformation during the process, we present a robust, quantitative approach to interpreting the stretching dynamics that involves consideration of the initial chain configuration in reduced dimensions.
Free energy landscape analysis.We now provide details on the analysis of the free energy landscape.We begin with an important technical note.Our learned GSOP following Eq.( 5) is in general not guaranteed to be a gradient system, unless W (Z) = 0 and M (Z) = I.However, since the drift term f (Z) = (f 1 (Z), f 2 (Z), f 3 (Z)) T = −(M (Z) + W (Z))∇V (Z), the stationary points of V are also critical points of the dynamics (f (Z) = 0).The saddle points of V are the saddle foci of the non-gradient dynamics Ż = f (Z).For simplicity, we refer to a saddle point of V and saddle focus of f interchangeably.
We compute the critical points by numerically solving ∇V (Z) = 0 with the BFGS method from different initial conditions.We obtain four critical points: (247.15,1.701, 0.193) (yellow point), (247.29,−1.700, −0.166) (blue point), (87.858, −0.041, − 4.269) (cyan triangle) and (85.887, −0.050, 4.126) (magenta triangle), which are shown in Fig. 4. The Jacobian matrix of the drift term is then computed We calculate the eigenvalues and eigenvectors of J at the four critical points.It has three eigenvalues with negative real parts at the yellow and blue points, so these are the stable points (stable nodes).We name them Z stable,1 and Z stable,2 .For the cyan triangle and magenta diamond points, J has one eigenvalue with positive real part, and two with negative real parts.These points are saddle points (saddle focus) of index 1.We call these two points Z saddle,1 and Z saddle,2 .The direction of the eigenvector corresponding to the eigenvalue with positive real part is the unstable manifold, along which the trajectories escape from the saddle.On the other hand, the span of the real and imaginary parts of the other eigenvector (one must be a complex conjugate of the other) constructs the stable manifold, along which trajectories are attracted to the saddle.
The local behavior of the learned potential V (Z) around its critical points characterize the typical fluctuations around them.There lies important physical meaning that we can probe.For example, let us exploit automatic differentiation of neural networks and expand V (Z) in a Taylor series around Z stable,1 (fully stretched state), corresponding to (247.15, 1.701, 0.193).Neglecting small terms, we obtain the approximate formula where δZ i = Z i − [Z stable,1 ] i is the fluctuations in the thermodynamic variables.Now, assuming that the distribution of states around Z stable follows a Constructing Custom Thermodynamics Using Deep Learning Boltzmann distribution Z ∼ exp[−V (Z)/k B T ] (here we assume that M ∼ I) where |σ| 2 ∼ k B T , the typical small fluctuations δV is proportional to k B T .In other words, Eq. ( 20) is approximate "isotherms" that captures the form of typical fluctuations.For example, we can infer from the formula that typical fluctuations of the extension length (Z 1 ) and the end-to-end distance (Z 2 ) are highly correlated.This is sensible, since in the fully stretched state, these two quantities are expected to change simultaneously.In Fig. 4(G,H,I), we confirm these correlations.
Similarly, one can also expand V (Z) around a saddle point (fully folded state) Z saddle,2 .We obtain the formula δV ≈ 102.96δZHere, we immediately observe that to escape the saddle point (lowering the energy), one should increase end-to-end distance (Z 2 ) and decrease foldedness (Z 3 ).This approximately aligns with the unstable manifold described above, and forms the basis of our control protocols described in the main text.
Selection of the dimension of the reduced coordinates.In this part, we describe how we arrived at the selection of a 3-variable reduced coordinate space.We tested various different number of reduction dimensions (d = 2, 3, 4) and the relative errors of predictions are summarized in Supplementary Table 1.We observed that going to a higher dimension (d = 4) did not result in noticeable gains.In fact, increasing dimension may cause increasing optimization error and hence worse results in standard deviation.Going to a lower dimension (d = 2) resulted in increased prediction error in general.
Interestingly, we can formulate a physical argument that suggests that a 2dimensional system (with only 1 additional closure coordinate) is not a suitable reduced model for the polymer dynamics we study.The argument is based on index theory for 2-dimensional dynamical systems [50].
Let us assume for the sake of contradiction that the projected free energy landscape of our learned 3-dimensional system into Z 1 − Z 2 plane is the free energy of a system in 2-dimensional (Supplementary Figure 10).Since Z 1 is the polymer extension, we expect that 1.There exists 2 stable states at large Z 1 corresponding to the fully extended state.There are 2 of them due to reflection symmetry in the flow direction.Around this Z 1 value, all trajectories in the reduced space should converge to one of the stable steady states.2. There cannot be saddle points with the same Z 1 value, since it is close to the maximal extended chain length.
These conditions are enough to imply a contradiction in the following way.
The simple concept we use from index theory is the definition of the index of a closed curve in the phase space of a 2-dimensional dynamical system (See Strogatz [50], Chapter 6.8 for details).Let Γ be a closed curve in R 2 , and consider a dynamical system ż = f The index of Γ with respect to the dynamics f is defined as Intuitively, this is the sum of the angles of the force vectors across the curve Γ.
From index theory, we know that the index of a closed curve must equal to the sum of indices of critical points it encloses.Moreover, the index of a stable critical point is +1, and the index of a saddle is −1.Now, we draw a curve enclosing the two stable points as shown in Supplementary Figure 10.By condition 1 above, the vector fields should point inwards towards the interior of the curve, thus we can show (from Eq. ( 21)) that the index of this curve is +1.However, it can only enclose stable critical points due to condition 2, and consequently its index is at least +2.Thus we arrive at a contradiction, showing that the dynamical landscape in Supplementary Figure 10 is not possible in 2 dimensions.
The only remaining possibility is that Z 2 does not distinguish the two symmetric stable extended states.However, in this case to allow a saddle point (which is the key to inducing heterogeneity in unfolding) we must have another stable steady state at a different Z 1 value.This is unlikely from a physical viewpoint for non-self-interacting polymers, since we expect the only asymptotically stable state to be the stretched, fully extended state.
The impact of dataset size on predictive accuracy of S-OnsagerNet.In order to study the impact of the number of training data on the computational accuracy of the model, the original training dataset containing 610 trajectories was split into two subsets containing 25% and 75% of the data.A third subset was produced by selecting at random 50% of the data from the original dataset.The three datasets and the learned potential landscapes are presented in Supplementary Figure 11.The 25% and 75% datasets have no common trajectories, whilst the 50% dataset contains some trajectories from each of the other two datasets.The quantitative results are reported in Supplementary Table 2, where the trained models are used to predict 500 unseen trajectories with fast, medium and slow unfolding times, suggest that with a smaller amount of training data, the S-OnsagerNet model loses some predictive accuracy.However, the more important factor is the diversity of the data contained in the training datasets.The full dataset has a high proportion of trajectories with fast, followed by middle (slower), and then slow unfolding times, however it is the most balanced of all the training datasets (25%, 50%, 75%, 100%), in addition to containing the most data.Training with the 25% dataset, which contains the largest proportion of trajectories with fast unfolding times results in the lowest L 2 error, but the model over-fits to that type of trajectories, and has the largest L 2 error for the medium and slow unfolding times.On the other hand, the 50% dataset is the most balanced of the reduced datasets, which is reflected in the relative L 2 errors.It can also be observed in Supplementary Figure 11 (bottom) that the potential landscapes resulting from training with smaller datasets do not contain a clear stagnation (saddle) point.Overall, for training datasets with different number of trajectories and their respective proportions of fast, medium and slow trajectories, the results, both in terms of prediction error and potential landscape characteristics suggest that the model is relatively robust, and the amount of training data cannot be reduced much without affecting predictive performance.

Spatial epidemics analysis
In this section, we provide details of our analysis method for an alternative application of S-OnsagerNet -modelling the macroscopic dynamics of the spread of epidemics.This highlights the general applicability of our method.
We focus on the most well-known model for disease spread in a spatial domain -the spatial SIR model [32].Let us consider a two-dimensional square domain (representing a city, say) discretized into n × n sectors.We use I i,j and S i,j represent the number (density) of infective and susceptible individuals at spatial location (i, j).The basic mechanism of the model is as follows: each infective individual may infect a susceptible individual in the same spatial location.At the same time, each infective recover (or is removed) at a rate, after which they are no longer infective.Finally, both infective and susceptible individuals move on the spatial domain randomly.Mathematically, the temporal evolution of I, S (understood as n × n matrices, or length n 2 vectors) are governed by the following dynamics Ṡi,j = − βI i,j S i,j + δ δ 2 x İi,j =βI i,j S i,j − γI i,j + δ δ 2 x As usual, the dot denotes time derivative.The parameter β is a measure of the transmission efficiency of the disease from infectives to susceptibles, and 1/γ is the life expectancy (or expected recovery time) of an infective.The constant δ is the diffusion coefficient, and this term in the equation models the spatial movement of individuals as a diffusion process over the domain.The last terms of the equation models the stochastic fluctuations of the number densities, with σ as the noise intensity.The parameters δ x and δ y are the spatial discretization sizes in the two spatial directions.In our simulations, we take n = 16, β = 0.3, γ = 0.13, δ = 0.5, σ = 0.03 and δ x = δ y = 2 3 .Eq. ( 22) governs the microscopic dynamics of disease spread, and non-trivial outcomes can result from different initial spatial configurations and parameters (e.g.infection rate, recovery rate).See Extended Data Figure 1.
While this microscopic model and its variants has been subject to intense study (see Murray [32] and references therein), a macroscopic understanding of the dynamics of disease spread is challenging due to the complex spatial interactions.For example, one may be interested to model the dynamics of average (or total) number of infective and susceptible individuals over the spatial domain.Observe in Extended Data Figure 1 that configurations with identical initial spatial averages of infective and susceptible individuals can have drastically different subsequent evolution.More precisely, one spatial configuration Extended Data Figure 1(A) can lead to initial disease spread (epidemic), where the mean number of infected individuals initially increases sharply, whereas another spatial configuration Extended Data Figure 1(B) leads to the disease dying out monotonically.Thus, it is of interest to develop a thermodynamic description to capture and elucidate the driving factors of such variations.
Modelling the thermodynamics of the spatial SIR model using S-OnsagerNet.Following our framework, we now set the macroscopic variables of interest as the respective spatial averages Z 1 = δ x δ y n i,j I i,j and Z 2 = δ x δ y n i,j S i,j .Recall from Extended Data Figure 1 that these variables alone are insufficient to determine their subsequent evolution.Our goal is to learn closure variable(s), and a stochastic dynamics that describes the evolution of these variables.
Training data is generated through integrating Eq. ( 22) using the Euler-Maruyama method with time step size dt = 0.03.The initial conditions are selected as I i,j (0) = 5e −(xi−x0) 2 −(yj −y0) 2 and S i,j (0) = 5e −(xi−1) 2 −(yj −1) 2 + 5e −(x−1) 2 −(y+2) 2 , where x 0 and y 0 are randomly generated from a uniform distribution in the unit square, and x i = −5+iδ x and y j = −5+jδ y .This initial condition corresponds to the scenario where the initial susceptible population is fixed, but the initial infective population is a cluster that is uniformly and randomly distributed in the domain.Then, we carry out the S-OnsagerNet workflow as shown in Fig. 1, with β = 0.01, α = 0.001, m = 50 and d = 3 (i.e. one closure coordinate).We employ a PCA-encoder network to obtain the closure coordinate.The training process involves initially training the PCAencoder and the S-OnsagerNet simultaneously.Subsequently, we fix the former and continue training the latter to a desired error tolerance.Note that in the SIR model case, we do not utilize the decoder network as there is no need to obtain a reconstruction of I and S.
Capturing the stochastic dynamics.First, we show that with just one additional learned closure coordinate Z 3 , we can capture the statistics of the macroscopic dynamics of the spatial averages of infected and susceptible individuals.The results are shown in Extended Data Figure 1(D,E), where true mean and standard deviation of Z 1 and Z 2 are obtained from Eq. ( 22), while the predicted results are derived from the S-OnsagerNet.Four representative test initial conditions are shown: two with disease spread and the other two without.We observe that we can successfully capture the macroscopic dynamics of disease spread with just one additional closure coordinate.
Interpreting the closure coordinate.Next, since we have shown that only one closure coordinate is required for a thermodynamic description, it is natural to probe its physical meaning.We use the same technique described in the polymer stretching case, where we investigate the effect of perturbation of a microscopic state in the directions that cause the sharpest changes in Z 3 .The microscopic state we probe is chosen as a pair of partially overlapping clusters of infective and susceptible individuals (Extended Data Figure 2(A)).We observe that the perturbations induced by dZ 3 /dX (where X = (I, S)) correspond to increasing/decreasing the overlap of clusters of susceptible and infected individuals.Hence, this suggests that Z 3 is a macroscopic descriptor that correlates with such effective spatial overlap.We confirm this hypothesis via a scatter plot in Extended Data Figure 2(B), where the overlap is defined by IS mean = δ x δ y i,j I i,j S i,j .Thus, we can interpret the closure coordinate Z 3 as an indicator of the overlap of clusters of susceptible and infected individuals.This is physically sensible, since the measure of overlap of spatial clusters should determine the outcome of an epidemic.Nevertheless, we must emphasize that the learned Z 3 is a quantitative measure and can be applied to more complex configurations than a pair of clusters, for which one may not be able to easily define a notion of effective overlap by empirical observation.
Free energy landscape.Finally, we study the dynamical landscape of the learned S-OnsagerNet model.Using the same projection technique in the polymer case, we plot 2D projections of the learned 3D free energy landscape in Extended Data Figure 3, overlaid with two representative trajectories (with disease spread in blue and disease dying out in red).
A number of interesting features can be gleaned from the landscape.First, we can clearly see the origins of the divergence of the two different types of trajectories.While they have identical initial Z 1 (average infective) and Z 2 (average susceptible) values, their initial Z 3 (infective/susceptible overlap) values differ (Extended Data Figure 3(B,C)).In particular, the initial disease spread seen in the red and blue trajectories is approximately in accordance to the steepest descent of the energy landscape.That is, the initial Z 3 value is a determining factor for the onset of epidemics.Second, we compute using the steady states of the dynamics corresponding by solving ∇V (Z) = 0. Instead of isolated steady states as in the case of polymer dynamics, we find a 1D manifold of stable steady states in the Z 2 − Z 3 plane (see Extended Data Figure 3(C,F)).This implies in particular that the terminal state (the remaining number of susceptible individuals) is not unique, but rather depends on the initial configuration, and in particular on the initial overlap described by the learned coordinate Z 3 .This rationalizes the observed heterogeneity in the terminal configurations as shown in Extended Data Figure 1.

Statistics & reproducibility
The train-test splits in this paper are performed by random uniform subsampling.Repeats of numerical experiments are performed by running the same code with different random seeds.No data were excluded from the analyses, and the investigators were not blinded to allocation during experiments and outcome assessment.

Figure legends/Captions
Fig. 1 Overall workflow of the proposed approach.Given a complex system described by X, the goal is to model the behavior of macroscopic coordinates of interest Z * .We construct closure coordinates Ẑ and closed (dynamical) equation on the combined reduced coordinates Z = (Z * , Ẑ).The classical ideal-gas law is an illustration of this process (top panel), where k B is the Boltzmann constant.For general non-equilibrium, dynamic systems (bottom panel), carrying out this workflow from theoretical analysis is challenging.Our machine learning method (mid panel) addresses this by simultaneously constructing the closure coordinates using PCA-ResNet (see Methods), and governing equations on reduced coordinates using the S-OnsagerNet with drift term −(M (Z) + W (Z))∇V (Z) and noise term σ(Z) (see Eq. ( 2)).The statistics of chain extension projected along the elongational axis are recorded.Times are reported in units of the characteristic rod diffusion time τ d (see data preparation in Methods).Different initial conditions (colors) are chosen to have similar initial extension but varying unfolding times.Identical initial configurations also have different unfolding dynamics due to thermal fluctuations.In (E-P), we show that S-OnsagerNet can capture this heterogeneity.We plot in (E) 500 trajectories of polymer extensions (Z 1 ) from the same initial condition, together with their mean (F) and standard deviation (G).The probability density functions (PDF) of unfolding times is shown in (H).(I-L) and (M-P) show results for two other unseen initial chain configurations.Fig. 4 Learned potential energy landscape.We plot V projected onto Z 1 − Z 2 (A,D), Z 1 − Z 3 (B,E) and Z 2 − Z 3 (C,F) planes (inset: stable and unstable directions of the saddle points).Projection is computed via minimization (e.g.V (Z 1 , Z 2 ) = min Z 3 V (Z 1 , Z 2 , Z 3 )), which at low temperatures closely approximates marginalizing with respect to the Boltzmann distribution.The stable and saddle points are marked on the energy landscape, and their corresponding reconstructed fully extended and folded states are shown.A pair of each exists due to reflection symmetry in the flow direction.Example "fast" (red) and "slow" (blue) trajectories from the training data set are overlaid on the landscape.The fast trajectory avoids the saddle points and goes directly towards a stable minimum, whereas the slow trajectory gets trapped for long times near saddle 2, before finally escaping through its unstable manifold.For (B,E), the stable manifolds of the saddles closely align with Z 2 , and hence are not visible due to minimization (marginalization).(G-I) Scatter plot together with predicted isotherms (solid lines) capturing typical fluctuations around a fully stretched state Z stable,1 .The insets are magnified views of the fluctuating trajectories around the stretched state over the energy landscape.Fig. 5 Data-driven control of the stretching dynamics.(A) Control protocol to speed up unfolding.Projected along the extension direction, the polymer must overcome energy barriers to transition from the folded to the stretched state (left).In the control protocol, the flow is turned off if the reduced coordinate of the polymer is near the saddle point corresponding to the folded state, and the polymer drifts under the Brownian motion (middle).Then the flow is turned back on if the reduced coordinates of the polymer becomes sufficiently aligned (green shaded region in (B)) with the unstable manifold of the saddle or when the the equilibration time reaches 100τ d (right).Without any control, a folded state near the saddle point will unfold slowly (grey lines in (C)); with control, the chains unfold more rapidly (green lines in (C)).For the slowest 10 trajectories shown, their mean unfolding time was reduced by 14.14%.(D,E) Reversed control to slow down unfolding by turning on the flow when the reduced coordinates become aligned with the stable manifold instead (blue shaded region in (D)).The mean unfolding time increased by 14.96% (blue lines in (E)).
Fig. 6 Analysis of experimental data.(A-B) Schematic and photo of experimental setup, consisting of a microfluidic cross-slot device and platinum electrodes in the reservoirs.Via the electrodes, positive voltage levels V 1 and V 2 are applied to the West and East reservoirs (W and E, respectively), whilst the North and South reservoirs (N and S, respectively) are kept at 0V.During trapping, the negatively charged DNA will thus flow from N/S to W/E, until a molecule is trapped at the center of the channel -blue dot in the schematic.(C) Snapshots of a DNA molecule stretching.(D,F) Processed experimental images at various percentages of the unfolding time (t uf ).The selected trajectories have similar initial configurations and are visually indistinguishable in terms of unfolding dynamics.(H) Learned potential landscape and predicted slow and fast trajectories using S-OnsagerNet with the initial configurations at t = 5%t uf from (D,F).We note slight differences in the initial Z 2 values only (inset).(E,G) Reconstructed high dimensional configurations of selected simulated trajectories with similar low dimensional coordinates as the experimental configurations in (D,F).The S-OnsagerNet is capable of distinguishing between the manually classified dumbbell and folded states.(I) Predicted probability density functions of unfolding time with initial condition of fast and slow experimental trajectories at t = 5%t uf using S-OnsagerNet.(J) Fluctuations in Z 1 and Z 2 around the stable (stretched) state from experimental images.Data was obtained from "Soh et al. [49]" (triangle) and "Soh et al. [31]" (square).Reduced coordinates are constructed according to the procedure outlined in data preparation in Methods.Colors indicate predicted energy levels according to the learned potential.We observe that fluctuations in Z 1 and Z 2 are highly correlated and agree well with that predicted by the effective equation of state.

Fig. 2
Fig.2Simulation setup, data visualization, and predicted vs true stretching trajectories and their statistics.(A,B) The polymer chain is represented by a bead-rod model with bead diameter r (in units of b = 10nm) and rigid bonds, subjected to hydrodynamic and Brownian forces.(C,D) The statistics of chain extension projected along the elongational axis are recorded.Times are reported in units of the characteristic rod diffusion time τ d (see data preparation in Methods).Different initial conditions (colors) are chosen to have similar initial extension but varying unfolding times.Identical initial configurations also have different unfolding dynamics due to thermal fluctuations.In (E-P), we show that S-OnsagerNet can capture this heterogeneity.We plot in (E) 500 trajectories of polymer extensions (Z 1 ) from the same initial condition, together with their mean (F) and standard deviation (G).The probability density functions (PDF) of unfolding times is shown in (H).(I-L) and (M-P) show results for two other unseen initial chain configurations.

Fig. 3 2 ∂Z 2 ∂X
Fig. 3 Physical interpretation of learned closure coordinates.(A) Perturbation of the polymer chain X ± ε 2 ∂Z 2 ∂X from a given configuration, ε 2 = 100/∥ ∂Z 2 ∂X ∥ 2 .(B) Plot of projected end-to-end distance |r N,x − r 1,x | as a function of |Z 2 | for the training data.(C) Configurations of different polymer chains with decreasing |Z 2 | values.As |Z 2 | decreases, the projected distance between the chain ends decreases.(D) Perturbation of the polymer chain X ± ε 3 ∂Z 3 ∂X from a given configuration, ε 3 = 260/∥ ∂Z 3 ∂X ∥ 2 .(E) Plot of degree of foldedness |r 1,x + r N,x | as a function of |Z 3 | for the training data.(F) Configurations of different polymer chains with decreasing |Z 3 | values.As |Z 3 | decreases, the chains tend towards a more stretched state.