{'fraction_non_alphanumeric': 0.032336390965852155, 'fraction_numerical': 0.010598936252216142, 'mean_word_length': 4.681410116049923, 'pattern_counts': {'":': 1, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonlinear Intrinsic Variables and State Reconstruction in Multiscale Simulations

Finding informative low-dimensional descriptions of high-dimensional simulation data (like the ones arising in molecular dynamics or kinetic Monte Carlo simulations of physical and chemical processes) is crucial to understanding physical phenomena, and can also dramatically assist in accelerating the simulations themselves. In this paper, we discuss and illustrate the use of nonlinear intrinsic variables (NIV) in the mining of high-dimensional multiscale simulation data. In particular, we focus on the way NIV allows us to functionally merge different simulation ensembles, and different partial observations of these ensembles, as well as to infer variables not explicitly measured. The approach relies on certain simple features of the underlying process variability to filter out measurement noise and systematically recover a unique reference coordinate frame. We illustrate the approach through two distinct sets of atomistic simulations: a stochastic simulation of an enzyme reaction network exhibiting both fast and slow time scales, and a molecular dynamics simulation of alanine dipeptide in explicit water.


I. INTRODUCTION
The last decade has witnessed extensive advances in dimensionality reduction techniques: finding meaningful low-dimensional descriptions of high-dimensional data [1][2][3][4][5] . These developments have the potential to significantly enable the computational exploration of physicochemical problems. If the (high-dimensional) data Y(t) arise from, for example, a molecular dynamics simulation of a macromolecule in solution, or from the stochastic simulation of a complex chemical reaction scheme, the detection of a few good, coarse-grained "reduction coordinates" x(t) can be invaluable in understanding and predicting system behavior.
While the benefits from such reduced descriptions are manifest, a crucial shortcoming of data-driven reduction coordinates is their dependence on the specific data set processed, and not only on the physical model in question. It is well known that, even in the simple linear case of Principal Component Analysis 6 , different data sets on the same low-dimensional hyperplane in the ambient space will lead to different basis vectors -in effect, to different reduction coordinates x. While this can be easily rectified by an affine transformation (see, by analogy, the discussion in Lafon et al. 7 ), the problem becomes exacerbated when the low-dimensional space is curved (a manifold, rather than a hyperplane) and when different data sets are obtained using different instrumental modalities (such as when one wants to merge molecular dynamics data with, for example, spectral information). Clearly, the ability to systematically construct a unique and consistent reduction coordinate set, shared by all measurement ensembles and observation modalities, is invaluable. We will call these coordinates Nonlinear Intrinsic Variables. Embedding data in such a coordinate system allows us to naturally merge different observations of the same system; more importantly, it enables the construction of an empirical mapping between these different observation ensembles, allowing us to complete partial measurements in a test data set from a training data set that consists of different observations. To construct this empirical mapping and the associated observers, accurate interpolation tools must be available in the embedding space; to this end, we will demonstrate the use of a multiscale Laplacian Pyramid approach 8 .
We will illustrate our methodologies with two distinct examples. The first is a simulation of two Goldbeter-Koshland modules in an enzyme kinetics model using the Gillespie Stochastic Simulation Algorithm (SSA) 9 ; in certain parameter regimes, separation of time scales is known to reduce the ODE model of this kinetic scheme to an effective two-dimensional 2 description 10 . Although this example is rather simple, it will serve as an introduction to our techniques and highlight the main features of the algorithms. The second example is a molecular dynamics simulation (in explicit water) of a simple peptide fragment (alanine dipeptide) whose folding dynamics are known to be described through a small set of physical observables 11 . This example will allow us to compare our approach to more common techniques, such as diffusion maps 12 for dimensionality reduction and nearest neighbor interpolation for observation reconstruction. The remainder of the paper is structured as follows: in Section II we present the Nonlinear Intrinsic Variable formulation and the associated inference method. Section III contains our discussion of Laplacian Pyramids that is used for the completion of partial observations. In Section IV, the results of the application of the approach to simulation data from our two illustrative examples are presented and discussed.
We conclude with a summary and our perspective on open issues in Section V.

A. Overview
Let Y(t) be a high-dimensional measured process in R n consisting of n observable variables. We impose two critical assumptions. First, the measured process is assumed to be a manifestation in an observable domain of a low-dimensional diffusion process. Thus, it can be expressed by where f : R d → M is an unknown (possibly nonlinear) function, M ⊂ R n is a d-dimensional manifold, and x(t) is a diffusion process that consists of d underlying variables (with d n).
Second, the dynamics of the diffusion process in each of its underlying variables are described by normalized stochastic differential equations as where a i are unknown drift functions andẇ i (t) are independent white noises. The independence is our second critical assumption.
Given a sequence of samples Y(t), t = 1, . . . , T , we present an empirical method to construct a unique and consistent reduction coordinate set, represented here by x(t) 13 . Because the empirical method we will describe is independent of the observation function f , we refer to the coordinates of x(t) as Nonlinear Intrinsic Variables (NIV). The available samples Y(t) may be the result of different measurement functions f in various observable domains, or they may be partial measurements consisting of merely a subset of the coordinates of the observable domains. The idea is to empirically construct a NIV coordinate system driven entirely by measurements that is invariant to the observation function f (see Figure 1 for a schematic illustration). We remark that the available data should be "rich enough", i.e., consist of a sufficient amount of historical data with adequate variability, in order to obtain the full empirical model.

B. Mahalanobis Distance
Let C(t) be the covariance matrix associated with the measured sample Y(t). In practice, the covariance matrix can be estimated from a short trajectory of samples in time around the sample Y(t) by where µ(t) is the empirical mean of the short trajectory of samples. We define a Riemannian metric between a pair of samples using the associated covariance matrices as this is the Mahalanobis distance (and † denotes a pseudoinverse, as discussed below). As previously described, the covariance matrices convey the local variability of the measurements and are utilized to explore and learn the tangent planes of the observable manifold.
This information is then utilized in (4) to compare a pair of points according to the directions of their respective tangent planes. The Mahalanobis distance is invariant under affine transformations. Thus, by assuming that the observation function f is bi-Lipschitz and smooth, and by using local linearization of the function, i.e., is the Jacobian of f (x(t)) and (t) is the residual consisting of higher-order terms, it was shown by Singer and Coifman 13 that C(t) = J(t)J T (t) and that the Mahalanobis distance approximates the Euclidean distance between the corresponding samples of the underlying process to second order, i.e., This result implies that the Mahalanobis distance is invariant to the measurement function f , and hence, it yields the same distances between samples obtained under different observation functions or even partial observations. We would like to note that, in general, f being bi- Lipschitz implies that f is invertible (on the d-dimensional manifold M). However, in practice, determining whether f contains sufficient information and is "rich enough" to completely determine the underlying process is a non-trivial task. In this work, we exploit the fact that C(t) = J(t)J T (t), which implies that C(t) is an n × n positive semidefinite matrix of rank d, to empirically infer the dimension d. According to the spectrum of the local covariance matrices and their corresponding spectral gaps, we approximate the rank of the matrices. Consistent rank estimates among these local covariance matrices are taken to imply that the measurements are "rich enough", and hence, may be good indicators for the dimension d. Since the dimension d of the underlying process is typically considerably smaller than the dimension of the measured process n, the covariance matrix is singular and non-invertible; thus, we use the pseudo-inverse in (4).

C. Laplace Operator
The Mahalanobis distance described in Section II B enables us to compare observations in terms of the intrinsic variables of the associated underlying diffusion process. In this section, we show how to recover the underlying process itself from the pairwise Euclidean distances through the eigenvectors of a Laplace operator.
Let W be a pairwise affinity matrix (kernel) based on a Gaussian, whose (t, τ )-th element is given by where ε is the kernel scale, which can be set according to Hein and Audibert 14 and Coifman et al. 15 . Based on the kernel, we form a weighted graph, where the measurements Y(t) are 6 the graph nodes and the weight of the edge connecting node Y(t) to node Y(τ ) is W t,τ . In particular, a Gaussian kernel exhibits a notion of locality by defining a neighborhood around are weakly connected to Y(t). In practice, we set ε to be the median of the pairwise distances.
According to the graph interpretation, this implies a well-connected graph because each measurement is effectively connected to half of the other measurements 16 .
Let D be a diagonal matrix whose elements are the row sums of W, and let W norm = D −1/2 WD −1/2 be a normalized kernel that shares its eigenvectors with the normalized graph- The eigenvectors of W norm , denoted ψ j , reveal the underlying structure of the data 12 . Specifically, the i-th coordinate of the j-th eigenvector can be associated with an intrinsic coordinate j of the sample x(i) of the underlying process. The eigenvectors are ordered such that |λ 1 | ≥ |λ 2 | ≥ · · · ≥ |λ n |, where λ j is the eigenvalue associated with eigenvector ψ j . Because W norm ∼ D −1 W, and D −1 W is row-stochastic, λ 1 = 1 and ψ 1 is the diagonal of D 1/2 . The next few eigenvectors can be argued to describe the geometry of the underlying manifold 12 . However, some eigenvectors can be higher harmonics of the same principal direction along the data manifold. This is analogous to how the eigenfunctions cos x and cos 2x of the usual Laplacian in one spatial dimension and with no flux boundary conditions are one-to-one with the values of x for 0 ≤ x ≤ 1; one must check for correlations between the eigenvectors before selecting those that describe the underlying manifold geometry. The above steps to construct the nonlinear intrinsic variables are summarized in Algorithm 1.
Ignoring the higher harmonics, each retained eigenvector then describes an intrinsic variable for the data set of interest. We must normalize the eigenvectors from different data sets so that the resulting embeddings are consistent. We first scale the eigenvectors so that where T is the number of data points, to make the embedding coordinates invariant to the size of the data set. Still, the computed embedding eigenvectors, even for two identical data sets, may differ by a sign. Reconciling the signs for the embeddings of different data sets can be rationally done in several ways and is somewhat problem-specific.
For example, if the mean of the embedding is sufficiently far from 0, we can require ψ i > 0; alternatively, if there is a common region sampled by both data sets, the sign of each eigenvector can be chosen to optimize the consistency of the embeddings of the common region data. We will return to the issue of embedding consistency for different data sets in our 7 concluding discussion; for the moment, we will assume that our different sets sample the same region of data space in a representative enough way such that the correspondence between the sequences of retained eigenvectors for different embeddings is obvious.

III. LAPLACIAN PYRAMIDS FOR DATA EXTENSION
In this work, we are not only interested in extracting the underlying variables x from some (partial) observations Y, but also interested in extending high-dimensional functions on a set of points which lie in a low-dimensional space. More specifically, viewing the am- Let f : Γ → R n be a function that is known on a subset of points S ⊂ Γ. A coarse representation of f is generated using a coarse smoothing operator P 0 . The smoothing operator P 0 is a normalized, coarse Laplacian kernel, defined by where is typically the Euclidean distance, and the parameter σ 0 is set to be large compared to the values of d 2 (i, j). The application of P 0 to f yields a coarse representation of the function, which we denote by f 0 = P 0 (f ).
The difference δ 1 = f − P 0 (f ) is the input for the next iteration of the algorithm, which uses the smoothing operator P 1 , : : scale l: f l = P 0 (f ) + l k=1 P k (δ k ) As l increases, the approximation becomes more refined because P l ∝ e −d 2 (i,j)/2 −l σ 0 uses a Laplacian kernel of a finer width. The iterations stop when the difference between f and f l is smaller than a pre-defined error threshold.
−0.02(x − 4π) 2 + sin(x) + 0.5 sin(3x) + 0.25 sin(9x) 7.5π < x ≤ 10π that contains several scales. The coarse regions of the function (0 ≤ x ≤ 4π) are well approximated by a small number of scales. As the function becomes more oscillatory (π ≤ x ≤ 7.5π and 7.5π ≤ x ≤ 10π), a finer representation, and a larger number of scales l, is required to capture its behavior.
In this work, LP is applied to extend a high-dimensional function f : are constructed and later extended to create the multiscale approximations as defined in Eq.
8. The LP algorithm for the inverse mapping is summarized in Algorithm 2.

A. A Chemical Reaction Network
We first consider a chemical reaction network involving multiple enzyme-substrate interactions 10 . The reaction steps that comprise the network are The " * " denotes an activated form of a species, and the ":" denotes a complex formed between two species; the complexes E : S and S : E are not equivalent. There are 10 species in this reaction system. However, one can write four conservation equations (since total E, S, D, and F are all conserved) to reduce the system to 6 dimensions (which we order as S, E, E : S, S : E, D : S * , F : E * ). We consider a parameter regime in which the ODE approximation of this scheme exhibits a separation of time scales, so that initial conditions quickly approach a two-dimensional manifold. Details about the specific parameter values can be found in Appendix A.
Although the dynamics of chemical reaction networks are typically described by a system of ODEs, the ODEs are only an approximation that holds in the limit of a large number of molecules. When the number of molecules is small, the system is inherently stochastic and its dynamics can be simulated using the Gillespie Stochastic Simulation Algorithm (SSA) 9 ; at intermediate molecule counts, the chemical Langevin approximation 21 becomes useful. We can control the level of noise in our simulation by adjusting the volume V , and therefore, adjusting the number of molecules, in the system. We take the volume small enough so that we can still observe appreciable stochasticity in small simulation bursts, but large enough (in our simulations, we take V = 10 5 ) so that the underlying two-dimensional manifold is (relatively) smooth.
We generate 3000 random initial conditions Y 0 (1), . . . , Y 0 (3000) ∈ R 6 , enforcing that all concentrations must be non-negative. We evolve each point Y 0 (t) forward for 10 time units using the SSA to obtain a point Y(t) ∈ R 6 ; according to the time scales calculated from the linearized ODEs, 10 time units is sufficiently long for the initial points in the ODE system to converge to the two-dimensional manifold, but not long enough for the points to converge to a one-dimensional curve or to the final steady state (see Appendix A for more details). In our stochastic simulations, the initial points appear to converge to an approximate two-dimensional manifold (in expected value, see Figure 3). We consider Y = {Y(t) : t = 1, . . . , 3000} to be representative points "on" this apparent two-dimensional manifold. From each manifold point y ∈ Y, we run 20 short simulation "bursts", each for 0.2 time units. We denote the endpoints from the short simulations as Y burst (y).

13
E : S, and D : S * , i.e., The endpoints of the simulation bursts for the two data sets, Y burst 1 (y) and Y burst 2 (y), are defined analogously. We then estimate the covariances for each point in each data set as whereμ i (y) is the empirical mean of Y burst i (y).
We first demonstrate that NIV produces the same embeddings for Y 1 and Y 2 , even though the two data sets contain information of different chemical species. Figure 4 shows the two-dimensional NIV embeddings for the two different data sets; the embeddings appear visually consistent. We also note that both Y 1 and Y 2 contain points that are projections of Y(1500), . . . , Y(2000). We therefore compute the correlation between the embedding coordinates for these points common to Y 1 and Y 2 . We obtain a correlation of 0.97 and 0.95 for the first and second NIV, respectively, indicating that the two embeddings are in quantitative agreement with each other. We would like to note that both Y 1 and Y 2 are sufficiently high-dimensional ("rich enough") to allows us to recover the common underlying two-dimensional manifold.
We then use NIV together with Laplacian Pyramids to estimate the values of S : E and F : E * for Y 2 . Because Y 1 and Y 2 are measured for different components, there is no simple way to estimate S : E and F : E * directly in the observation space. Instead, we must first embed the data into the NIV space so that we can compute neighbors between the two data sets. We use Y 1 to train an LP function from the two-dimensional NIV embedding to S : E and F : E * . We then use this function to predict the values of S : E and F : E * for Y 2 , using the computed NIV embedding for Y 2 . In this way, we are exploiting the fact that the NIV embedding is intrinsic and consistent between the two data sets, even though the two data sets contain measurements of different chemical species.
The results of the LP prediction are shown in Figure 5. The normalized mean-squared errors between the true and estimated values for S : E and F : E * , defined as are 0.0372 and 0.0287 , respectively. Therefore, we can effectively estimate the unobserved components in the reaction network using NIV together with LP.

B. Alanine Dipeptide
Our second example comes from the molecular dynamics simulation of a small peptide fragment. Alanine dipeptide (Ala2) is often used as a "prototypical" protein caricature for simulation studies 11,[22][23][24][25][26] . We simulate the motion of Ala2 in explicit solvent using the AMBER 10 molecular simulation package 27 with an optimized version 28 of the AMBER ff03 force field 29 . The molecule is solvated with 638 TIP3P water molecules 30 with periodic boundary conditions, and the particle mesh Ewald method is used for long-range electrostatic interactions 31 . The simulation is performed at constant volume and temperature (NVT ensemble), with the temperature being maintained at 300 K with a Langevin thermostat 32 .
Hydrogen bond lengths are fixed using the SHAKE algorithm 33 . The two dihedral angles φ and ψ are known to parameterize the free energy surface, which contains three important minima (labeled A, B, and C, see Figure 6). Our simulations are concentrated around minimum B in the free energy surface, located at φ ≈ −65 • , ψ ≈ 150 • . We start many simulations at 10 • away from the minimum, and allow the simulations to each run for 0.1 ps, while recording the configuration of Ala2 every 1 fs (therefore, each trajectory is 100 points long). Configurations are recorded with all atoms except the hydrogens.
The correlation between the NIV coordinates for the two data sets and the diffusion map (DM) coordinates for the two data sets are shown in Figure 7. The correlation between the two NIV embeddings is higher than the correlation between the two diffusion map embeddings. Therefore, it appears advantageous to use NIV over diffusion maps if one wishes to obtain a consistent embedding and merge data sets from different observation domains (as long as the two main assumptions underpinning the NIV algorithm hold).
We  The NIV embedding for the training data Y all is shown in Figure 8. The embedding is three-dimensional, and visual inspection reveals that each coordinate can be directly linked with one physical variable: the first coordinate describes the flipping of atoms 1 and 3, the second coordinate describes the dihedral angle φ, and the third coordinate describes the dihedral angle ψ. We calculate the correlation between the embedding coordinates for the points in Y all and Y odd that come from the common simulation data points Y(8001), . . . , Y(10000).
The embeddings for the two data sets are found to be fairly consistent, with correlations of 0.97, 0.72, 0.85 for the first, second, and third NIV, respectively. Figure 9 shows the reconstructed position from partial observation versus true position for certain selected atoms. The strong correlation between the true and reconstructed positions is easier to appreciate for atoms that move substantially within the data set (such as atoms 1 and 3). Figure 10 shows molecular structures for the true and reconstructed configurations for selected data points; there is qualitative agreement between the true and reconstructed configurations.
For a brief comparison of LP over other reconstruction techniques, we also reconstruct configurations from the NIV components using simple nearest-neighbor interpolation. The average reconstruction error, scaled by the average bond length within the molecule, is shown in Figure 11; LP arguably outperforms simple nearest neighbor search for all of the atoms.

V. CONCLUSIONS
We have used Nonlinear Intrinsic Variables to analyze two complex atomistic simulations: a stochastic simulation of a chemical reaction network and a molecular dynamics simulation of alanine dipeptide.
In both examples, we were able to uncover the intrinsic variables governing the underlying stochastic process, which are independent of the particular measurement or observation of the system (under the conditions mentioned). The uniqueness of the embedding coordinates allowed us to compare and merge data sets from different measurement functions, and therefore allowed us to use an interpolation/extension scheme (here Laplacian Pyramids) to complete partial observations. Different interpolation techniques (e.g. kriging 34,35 , geometric harmonics 36 , versions of the Nyström extension) can and should be explored, since the performance of such techniques may well be problem dependent, especially for multiscale, complex simulation data.
There are many open questions leading to interesting research directions to be explored.
In this work, we considered data sets that consist of different partial observations, but in which each data set samples the entire underlying manifold in what we loosely referred to as a "representative enough" way. However, NIV could also be used to merge data sets when each data set samples only a portion of the manifold, provided there is enough overlap to "register" the embeddings. Merging data sets that come from different portions of the manifold would not only require scaling the embedding coordinates, but also shifting and possibly permuting the embedding coordinates (in this spirit, see the discussion in Lafon et al. 7 ). The ability to merge data from different regions would then allow us to analyze systems where complete sampling is computationally intractable, such as molecular systems with several high energy barriers separating regions of state space.
Other issues, such as accurately estimating the covariance matrices required for the computation of the Mahalanobis distance, are also of current research interest. It is clearly necessary to link this type of calculation with modern estimation techniques for (multiscale) diffusions [37][38][39] to test the appropriateness of the window sampling lengths selected; this will determine the accuracy of the noise covariance estimation by eliminating the bias due to drift variations. We are confident that the exploration of these open questions will enable the use of our methodology in many interesting applications, such as merging data from molecular simulations at different levels of granularity, or merging simulation data with experimental observations.

Appendix A: Chemical Reaction Network Parameters
We consider the following network of chemical reactions.
In the limit of a large number of molecules, the dynamics of this network is governed by the following ODEs. We choose to eliminate S * , E * , D, and F from the system of ODEs. We therefore obtain