Exploring chemistry and catalysis by biasing skewed distributions via deep learning

Domain: chemistry
Published: 2026-02-21
Editorial outcome: minor_revision

Abstract
The automated discovery of chemical and catalytic reactions remains a major challenge in computational chemistry, particularly in complex systems where conventional methods struggle to identify optimal searching directions. Here, we propose Loxodynamics, a machine-learning-driven approach for reaction exploration via biased molecular dynamics. By leveraging the skewness of local probability distributions, Loxodynamics dynamically determines low-energy barrier directions, efficiently guiding the system toward metastable states. The core of our framework is Skewencoder, an Autoencoder augmented with a skewness-based loss function that extracts reaction coordinates from minimal sampling data. Through iterative sample-and-search cycles, the system adaptively maps the free energy surface, capturing finite-temperature effects critical to complex reactive environments. We validate our method across model potentials, gas-phase reactions (SN2 and Diels-Alder), and catalytic alcohol dehydration in acidic chabazite under operando conditions. Loxodynamics provides a systematic framework for reaction discovery that addresses key limitations of conventional techniques, notably by allowing acceleration without elevated temperatures or a priori knowledge of collective variables.

Introduction
Exploring chemical reactivity from first principles is a grand challenge in computational chemistry1–3. This is fundamental across all fields of chemistry and particularly critical in catalysis, where identifying multiple products, intermediates, and reaction mechanisms is key. Even simple systems can generate a diverse range of possible outcomes, and the complexity further increases when reactivity is influenced by environmental factors such as interfaces, solvents, and catalytic conditions3,4.

The conventional workflow to tackle this problem typically involves inferring the structures of reactants and products based on known or anticipated reaction outcomes or relying predominantly on chemical intuition.5–18. Once potential products or intermediates have been hypothesized, a guess structure is constructed, and its energy is minimized through geometry optimization.

While this approach has achieved significant success, it faces two main limitations. First, designing the guessed product structure heavily relies on chemical intuition, which may not always yield realistic predictions. To mitigate this, various methods have been developed to explore the potential energy surface (PES). Many of these techniques leverage the local curvature of the PES to identify the minimum energy pathway connecting reactants, transition states, and products19–24. These methods have proven effective but are generally restricted to small, low-complexity systems where complex dynamic effects can be neglected.

Alternatively, for more complex scenarios, conceptual chemical knowledge can guide the identification of potential intermediates, such as by searching for a minimum energy path (MEP) connecting them. If the two endpoints of an elementary reaction are known, interpolation and string methods can be applied to locate the corresponding MEP6,8,25. Conceptual chemical insight can also be exploited to rapidly identify potential intermediate candidates that are connected to the starting structure via an elementary reaction—particularly when the type of reaction mechanism relevant to the system under study is known. These heuristic approaches have proven highly effective in characterizing products and intermediates in well-established processes17,26. However, their applicability is often constrained by the requirement for prior knowledge of the reaction pathways. The use of graph theory and connectivity analysis offers an efficient yet mathematically elegant alternative for exploring chemical reaction networks27, with promising applications in catalysis15,27.

For a detailed overview and classification of these reaction exploration approaches, we refer the reader to the excellent reviews by Reiher and coworkers2,28,29.

The second key challenge is that PES exploration methods very often rely on structure optimization techniques that overlook finite-temperature effects, thereby neglecting the broader configurational space accessed by the system. This limitation increases the risk of missing the true lowest-energy minimum corresponding to the reaction product or intermediate of interest.

To incorporate finite-temperature effects in reactivity studies, molecular dynamics (MD) simulations can be employed. However, standard MD operates on timescales much shorter than those of activated processes like chemical reactions, making it challenging to capture such rare events.

To overcome this limitation, various enhanced sampling methods have been developed to promote transitions across high-energy barriers separating metastable states. Those strategies involve increasing the system’s temperature30–32, or adding a bias potential to the global PES of the system33–35, in order to boost the likelihood of barrier crossing. These approaches often employ multiple replicas, as in parallel tempering, or constrain molecules within a small fictitious volume36 while applying elevated temperatures. While effective in facilitating exploration, these methods can also introduce artifacts, leading to the identification of high-energy pathways and reaction products that may not be experimentally relevant under realistic conditions.

Another widely used approach involves applying bias potentials along predefined collective variables (CVs) to accelerate transitions between states. This category includes umbrella sampling37, metadynamics38, and On-the-fly Probability Enhanced Sampling (OPES)39. While highly effective for sampling known reaction pathways, these methods rely heavily on prior knowledge of metastable states, needed to design the CVs. This dependence makes them less suitable for reaction discovery.

In general, designing a fully agnostic CV for reactivity exploration with a high sampling efficiency is impractical without prior knowledge of the expected products. Moreover, even with this information, selecting the correct bias direction along the CV is essential to guide the system along the lowest-energy reaction pathway and efficiently drive the transition to the products.

Here, we propose a method that overcomes these limitations by automatically constructing a low-dimensional CV while simultaneously determining the optimal direction toward the most probable reaction products. By following pathways with the lowest energy barriers, this approach enables efficient and agnostic exploration of likely reaction channels. To achieve this, we draw inspiration from two key concepts: the Anharmonic Downward Distortion Following (ADDF) method19,21 and the Artificial Force Induced Reaction (AFIR) method40–43, combined with pulling-based exploration techniques used in steered MD44,45 and ratcheted-and-pawl MD46 simulations.

The ADDF method explores reaction pathways by following anharmonic distortions by mapping on the PES based on the concept of hypersphere search19,47, while the AFIR method identifies associative reactions using artificial forces between reactants. Both have been successful in mapping chemical reactivity, but neglect temperature effects.

Conversely, dynamic approaches like steered MD and ratcheted-and-pawl MD enable exploration of activated processes but require prior knowledge of initial and final states to define an efficient CV. Additionally, the pulling direction must be predefined by the user.

Our approach unifies these concepts by automatically identifying the CV along which anharmonicity is maximal, thereby guiding the exploration along the slowest ascending directions of the PES. Instead of explicitly evaluating high-dimensional PES directions, the method first samples the local minimum of the reactant state. By analyzing asymmetries in the local distribution—specifically its skewness—it infers the optimal biasing direction. This eliminates reliance on predefined knowledge, enabling automated reaction pathway discovery and efficient sampling of the most probable reaction channels.

Results
Loxodynamics

The essence of our method can be illustrated, for the sake of simplicity, by considering the dynamics of a finite mass particle in a one-dimensional double-well potential (Fig. 1) at finite temperature, e.g., Langevin dynamics. Let us imagine that the shape of the potential is not known and the temperature is sufficiently low, i.e., k B T≪E act, so that the particle is trapped in the left basin. In such a situation, it is not trivial to decide whether the particle should be pulled toward the left or the right along x to go to the other state. Our key idea is that the local potential experienced by the particle inherently exhibits some degree of anharmonicity, which, in turn, must be reflected in the shape of its probability distribution. By sampling the vicinity of the minimum for a sufficient duration, it is possible to estimate the local probability distribution function and its moments, such as the skewness. The latter, for a one-dimensional distribution, is defined in Equation (1) as:

γ = 1/n * sum_{i=1}^n ((x_i - x_mean)^3) / ((1/n * sum_{i=1}^n (x_i - x_mean)^2)^(3/2)) (1)

where n is the number of samples, x_i represents the i-th value of the random variable, and x_mean is the sample mean. A positive value of γ indicates a left-modal distribution, meaning that the tail of the distribution becomes more pronounced as the random variable x increases. Conversely, a negative γ signifies that the tail is more pronounced as it decreases. In the context of an asymmetric anharmonic potential, a positive skewness implies that shifting the mass along the positive direction of x is likely to guide the system across the lowest energy barrier, where the tail of the distribution is longer, thus leading it into a new metastable state. For this reason, we named our method Loxodynamics, deriving from the Greek word loxós, meaning skewed, to emphasize its core principle of leveraging skewness in probability distributions to guide the sampling process.

In this model system, the initial sampling shows positive skewness, indicating that the particle should move to the right to reach the other basin. To implement this, a harmonic wall is applied, shifted by a constant offset Δx to the right of the mean μ of the distribution sampled in the previous step, and sampling resumes. The new potential, which combines the physical double-well potential and the restraining bias, makes the particle sampling a new landscape and search iteratively for the next metastable state direction based on the skewness of the updated distribution. The lower panel of Fig. 1 shows the evolution of the local equilibrium distribution after applying the restraining biases at each successive iteration. In this simple case, the skew remains positive and, therefore, the direction along which the bias must be applied is always positive along x. In this way, through successive iterations of sample-and-search, the system gradually exits the left basin and reaches the right one, completing the exploration process.

Real chemical systems often require multiple dimensions to capture relevant anharmonicity changes in the potential. A practical approach is to define N_d descriptors d_i(R), functions of atomic coordinates R, which serve as the basis for constructing the CV. The procedure remains similar to the 1D case: after defining these descriptors, the system is briefly sampled at a given temperature, generating an N_d-dimensional distribution. However, extending skewness to higher dimensions is non-trivial due to the complexity of quantifying asymmetry48.

To leverage the benefits of 1D skewness, dimensionality reduction is necessary. This reduction must preserve key system features while maintaining differentiability for CV construction. A tailored bias can then be applied in this reduced space, enabling efficient exploration of the potential energy landscape.

To achieve this, we developed Skewencoder, an Autoencoder-based framework that reduces data to a 1D latent space while incorporating higher-moment statistical information. Autoencoder-based CVs have been widely adopted as an initial means to identify essential features relevant to reaction mechanisms and molecular transformations, thereby serving as a foundation for reaction space exploration. Chen et al.49,50 developed an iterative exploration protocol that integrates umbrella sampling with dynamically updated Autoencoder-based CVs derived from the evolving set of sampled configurations. Their approach employs a grid-based boundary detection method to determine bias placement by identifying outlier configurations potentially indicative of state transitions. In comparison, our framework—while also employing iterative CV-accelerated discovery—prioritizes the exploitation of collective statistical properties within the sampled data, such as skewness, to guide and optimize transition pathways. Ketkaew et al.51,52 introduced an autoencoder-based strategy aimed at constructing data-driven collective variables for enhanced sampling, trained exclusively on reactant-state configurations. Their approach exploits deep autoencoders to extract latent representations from unbiased trajectories and employs these representations within flooding-based enhanced sampling schemes to reconstruct free energy surfaces. While this methodology alleviates the need for explicit product-state information in the definition of collective variables, it remains inherently tied to static CVs and relies on conventional enhanced sampling paradigms, without introducing an explicit directional bias capable of actively driving exploration toward new regions of phase space.

The aforementioned approaches combine CVs with enhanced sampling techniques such as umbrella sampling, metadynamics, or OPES to accelerate phase-space exploration. However, these methods do not encode an intrinsic directionality in the applied bias53. In contrast, Loxodynamics introduces a bias explicitly designed to drive the system toward new states along a CV that is updated on-the-fly during the exploration. Metadynamics and OPES, in particular, are flooding approaches, in which the bias progressively fills local free-energy basins defined along static, pre-defined CVs, thereby enhancing fluctuations within those basins. The core principle of Loxodynamics is fundamentally different: rather than amplifying fluctuations, it deliberately focuses them with the explicit goal of propelling the system uphill in free energy.

By integrating an Autoencoder neural network (NN) with a loss function based on the skewness of projected samples, Skewencoder efficiently utilizes limited early-stage sampling data in a multitask learning framework, as defined in Equation (2).

L_total(x; θ) = L_AE(x; θ) + αL_skew(x0; θ) + β||w||^2 (2)

where L_AE denotes the Autoencoder loss function, and L_skew is a skewness loss term designed to maximize the skewness of the latent space distribution. Here, x and x0 represent the global and local set of input descriptors, respectively, and θ comprises all neural network parameters (i.e., weights and biases across all layers). The term ||w||^2 refers to the L2 regularization applied to the weights, while α and β are hyperparameters that control the relative contributions of the skewness and regularization terms, respectively.

The Autoencoder loss function is defined in Equation (3) using the standard reconstruction error, based on the difference between the inputs and the reconstructed outputs:

L_AE(x; θ) = 1/n * sum_{i=1}^n ||xe_i(x_i; θ) - x*_i||^2 (3)

where n is the number of training examples, x* is the vector of inputs for the i-th training sample, xe_i is the output vector of the Autoencoder, which is a function of the NN parameters θ and the input descriptor x_i and the difference is calculated as the mean square error (MSE) of the samples, expressed in terms of vector norms. The global dataset accumulated across iterations is used for optimizing the Autoencoder.

In contrast, the skewness loss, which is calculated based on the projected CVs in the latent space (depicted by the light-blue path in Fig. 2), is defined in Equation (4) as:

L_skew(x0; θ) = ln(1 + e^(-γ(s(x0; θ))^2)) (4)

where s(x0; θ) denotes the latent representation computed by the Autoencoder NN, γ(s(x0; θ)) denotes the skewness, computed according to Equation (1), and is derived solely from the local data x0 generated by the MD simulation in the most recent iteration (refer to the Supplementary Information (SI) for further details). The use of both global and local datasets for calculating the loss functions defines the multitask nature of our method.

In the early iterations of the sample-and-search process, designing an effective CV is challenging due to the limited data available for training the Autoencoder. To address this, we employed a warm-start training approach, which has been widely successful in convex optimization54. Instead of training the model from scratch at each iteration, the parameters learned from the previous iteration are used as the starting point, and the model is trained to converge from these pre-initialized values (see SI for details).

Case study

2D model potential. As a proof of concept, we tested our method on simple 2D model potentials using Langevin dynamics. For clarity, we present results for one of the three tested potentials, shown in Fig. 3 (panel a). The explicit functional form and results for the other cases are provided in the SI.

Starting in the upper basin, the algorithm progressively sampled the potential. The panel b of Fig. 3 highlights key iterations of the exploration, while the panel c shows the trajectory along the x and y coordinates, illustrating how the system systematically gathers information and transitions toward the second basin.

This potential features two pathways with barriers of different heights. The algorithm iteratively analyzes the data distribution, identifies the direction of maximal anharmonicity—where skewness is highest—and applies harmonic walls to bias the system along the lower-energy route, ensuring a smooth and efficient transition to the next metastable state.

In detail, the blue points in panel b of Figure 3 represent samples from the incremental dataset, which aggregates data from all previous and current iterations to form the global training set for the Autoencoder’s reconstruction loss. In contrast, red points correspond to samples from the current iteration, constituting the local training dataset optimized for the skewness loss. Sampling at each iteration follows these criteria, progressively expanding the explored region and allowing the system to identify new minima on the free energy landscape.

In this example, both pathways lead to the same product. However, in realistic chemical systems, a single reactant can yield multiple products, each requiring distinct transition regions. Our method ensures that the system follows the lowest-barrier pathways, identifying the most probable reaction routes and products.

Moreover, upon algorithm termination, the accumulated samples provide sufficient data to construct a latent space that effectively approximates a meaningful CV (see SI), ensuring an accurate representation of the reaction pathway, even in complex energy landscapes.

Gas phase chemical reactions. To demonstrate the capability of Loxodynamics in capturing chemical reactivity, we selected two prototypical reactions: the SN2 nucleophilic substitution of chloromethane by a fluoride anion and the gas-phase Diels-Alder reaction between ethylene and 1,3-butadiene. These reactions exemplify key aspects of reactivity, including simultaneous bond rearrangement (synchronous bond breaking and formation) and significant differences in the shape of the reactant and product basins as well as asymmetric barrier heights in the forward and reverse directions.

To detect whether the system has relocated into a new chemical metastable state, bond topology changes were assessed by monitoring variations in contacts between neighboring atoms (see “Methods” section for details). In addition, how metastable states are identified in different case studies is also discussed (see SI).

To begin, we investigated the SN2 nucleophilic substitution of chloromethane by a fluoride anion. The reaction is characterized by two key distances: d_F, the distance between the carbon and fluoride ions, and d_Cl, the distance between the carbon and chlorine ions. These distances were used as inputs for our Skewencoder. To ensure consistent training, all distances were normalized to zero mean and unit variance.

Figure 4a,c illustrate the progression of Loxodynamics, tracked via the control variable d = d_Cl - d_F, which provides a simplified yet effective one-dimensional CV to control the reaction progression55. After a few iterations, the Skewencoder identified the most favorable direction, enabling the bias to efficiently drive the system toward the new metastable state.

Due to the asymmetric energy barriers between reactants and products, a larger harmonic constant was applied when steering the reaction in the direction Cl- + CH3F -> CH3Cl + F- to account for the higher barrier. The method consistently determines the optimal pulling direction within the two-dimensional plane defined by key interatomic distances.

This highlights the robustness of the approach: even with a simple basis of chemical bond descriptors, it effectively identifies the reaction pathway with the lowest barrier by targeting the most anharmonic directions on the PES. Statistical asymmetries in locally sampled distributions are leveraged after short sampling trajectories.

Subsequently, we extended the application of our method to the more complex chemical case of the Diels-Alder reaction between 1,3-butadiene and ethylene. For this system, we deliberately chose to include all 15 possible interatomic distances among the heavy (carbon) atoms in the molecule as descriptors (see SI), rather than preselecting only the two most relevant C-C distances involved in the reactive process. This decision, made regardless of whether these distances directly describe relevant bonds in either metastable state, was aimed at testing the robustness and general applicability of the method.

Figure 4b,d illustrate the progression of Loxodynamics for the Diels-Alder reaction. The method successfully identifies key reaction pathways in both forward and reverse directions without relying on predefined mechanistic knowledge.

The forward reaction presents significant challenges due to its broad reactant state with high entropy, where the conformational variability of the separated fragments complicates the identification of an effective biasing direction. As shown in Figure 4b, after a few iterations, the method successfully aligns the fragments, facilitating σ bond formation in the cyclohexene product state. By iteration 8, the system surmounts the energy barrier and transitions to the product state.

The reverse reaction is equally complex but for different reasons. The deep and narrow ring state exhibits low anharmonicity, making it difficult to determine the correct biasing direction without mistakenly exploring alternative anharmonic channels, potentially leading to undesired C-C or C-H rearrangements. Despite these challenges, the method effectively identifies the ring-opening pathway and smoothly guides the system out of the metastable state, as depicted in Figure 4d. This highlights the method’s robustness in traversing complex energy landscapes, enabling accurate and insightful reactivity exploration.

Catalytic ethanol dehydration in acidic chabazite. In this section, we apply our method to explore a highly complex catalytic system under realistic operando conditions. We applied Loxodynamics to a reaction occurring at the Brønsted acid site (BAS) of a zeolite. Here, we focus on the dehydration of ethanol to ethylene at the BAS of acidic chabazite (H-CHA), a widely used zeolite, a complex reaction, and highly relevant for the chemical industry.

Experimental56–58 and computational59,60 studies have shown that one mechanism is dominant, as it has a lower activation barrier: the concerted mechanism. In this pathway, the acidic zeolite protonates the hydroxyl group of the ethanol molecule, facilitating the elimination of water. Simultaneously, a hydrogen atom is abstracted from the ethyl group, regenerating the BAS proton and restoring the catalyst.

However, evidence suggests the existence of an alternative pathway with a slightly higher energy barrier, yet still kinetically feasible. This mechanism, commonly referred to as the stepwise mechanism, involves the formation of an ethoxy-framework complex intermediate. More precisely, initially, the elimination of water leaves a CH3CH2+ cation, which binds to a negatively charged BAS oxygen, forming a relatively stable species. Subsequently, a proton of the methyl group of the ethoxy group is abstracted by a neighboring BAS oxygen, oxidizing the intermediate and yielding ethylene.

To set up our simulation, we started from the adsorption state of ethanol at the BAS. To begin our exploration, our input set comprises 37 descriptors, including interatomic distances among heavy atoms, such as C-O and C-C, but also with light atoms C-H and O-H, with minimal reliance on chemical intuition (see “Methods” and SI for details).

To our surprise, despite the simplicity and generality of these descriptors, Loxodynamics successfully identifies the ethylene product state through two distinct pathways. In one case, it follows a direct route to the product, while in another, it first identifies an intermediate corresponding precisely to the ethoxy-framework complex, which subsequently evolves into ethylene via the expected mechanism.

Figure 5 illustrates two detected reaction pathways from the sample-and-search process, in agreement with the static calculations by Kim et al.59, corresponding to the concerted and stepwise mechanisms, respectively.

As shown in the upper panel, the method successfully identifies the concerted reaction pathway, in which the framework proton is transferred from the BAS to ethanol, leading to simultaneous water elimination and proton abstraction from the ethanol methyl group by an adjacent BAS acidic oxygen.

In contrast, the lower panel illustrates the stepwise mechanism, which first involves the formation of an ethoxy group bound to a BAS oxygen—a crucial intermediate in this pathway. To prevent back-crossing to the original state upon reaching this intermediate state, we restart the Loxodynamics simulation. While the Skewencoder parameters and the original CV are retained, an additional restraining wall is applied along this CV to confine the system to the current intermediate state. The simulation then proceeds with this constraint active. (For further details regarding the restart protocol, refer to the SI.) By restarting another Loxodynamics search from this intermediate, the second step of the stepwise mechanism is observed, consisting of proton abstraction from the methyl group by a vicinal BAS oxygen.

This demonstrates that our method is capable of identifying not only the most likely products but also of agnostically exploring reaction pathways that transit through intermediates with a similar probability of occurrence.

This is a fundamental aspect of studying chemical reactivity in complex catalytic systems, where experimental knowledge is often limited due to the short-lived nature of such intermediates. Nevertheless, these intermediates play a crucial role in determining the overall mechanistic network, which is essential for multiscale modeling of catalytic reactions, such as through microkinetic models.

Catalytic 1-butanol dehydration in acidic chabazite. Finally, we extend our approach to the investigation of another highly complex catalytic system involving a BAS in a zeolite framework. Specifically, Loxodynamics was applied to study the dehydration of 1-butanol to butene at the BAS of acidic chabazite. As observed in the ethanol dehydration mechanism discussed in the section “Catalytic ethanol dehydration in acidic chabazite”, both concerted and stepwise pathways are present, each kinetically accessible but characterized by slightly different energy barriers61. In the concerted mechanism for 1-butanol dehydration, analogous to that of ethanol, protonation of the hydroxyl group by the acidic zeolite proton promotes water elimination, while simultaneous abstraction of a hydrogen atom from the butyl group regenerates the BAS proton and restores catalytic activity. Conversely, in the stepwise pathway, formation of a butoxy intermediate occurs: initial water elimination produces a CH3CH2CH2CH2+ cation that associates with a negatively charged BAS oxygen to yield a relatively stable intermediate. Subsequently, abstraction of a methyl proton from this intermediate by an adjacent BAS oxygen results in oxidation and formation of an alkene double bond.

What radically differentiates the case of 1-butanol from ethanol is the wide variety of possible reaction products, arising from the position of the double bond along the alkyl chain and the associated isomerism. Both experimental studies62 and computational investigations63 have shown that the dehydration of 1-butanol can yield all the different dehydration products and, under conditions of high substrate loading, can also be accompanied by side reactions such as oligomerization, cracking, and aromatization. In the present work, we restrict our attention to single-site adsorption, focusing specifically on the dehydration pathways.

The simulation was initialized from the adsorption state of 1-butanol at the BAS. For exploration, the input descriptor set consisted of 28 features, encompassing interatomic distances among heavy atoms as well as selected light atom pairs, there by minimizing dependence on prior chemical intuition (see “Methods” and SI for further details). Fewer distances were incorporated to direct system evolution toward the intended reaction space.

Figure 6 illustrates the reaction pathways identified through the Loxodynamics sample-and-search procedure. Along the direct dehydration route, the catalyst generates the three canonical products of 1-butanol dehydration: 1-butene, and the cis- and trans-isomers of 2-butene. As anticipated, 1-butene formation may also proceed via stepwise mechanisms. In close analogy to the case of ethanol, this mechanistic motif is preserved for longer alkyl chains. Notably, the stepwise route exhibits high selectivity toward 1-butene, whereas the formation of cis- and trans-2-butene is not observed.

Despite employing a descriptor set of remarkable simplicity and generality, Loxodynamics reliably identifies the butene conformer product. Furthermore, the method uncovers an intermediate corresponding to a butoxy-framework complex, which subsequently evolves into 1-butene through the expected mechanistic pathway, without branching into alternative products.

These products and pathways have previously been reported by John et al.63 for the same reaction in zeolite H-ZSM-5, where static DFT calculations were employed to identify reactants, products, intermediates, and transition states. In their study, however, the initial guesses for all structures relied heavily on chemical intuition and required extensive manual effort, including considerable trial-and-error testing. By contrast, our approach offers a simple yet effective framework to explore complex reaction networks in highly dynamic systems such as zeolite catalysts under operando conditions.

This case study demonstrates that our approach is capable of efficiently exploring energetically favorable conformer product spaces, as well as identifying both the most probable products and reaction pathways involving intermediates with comparable likelihoods of occurrence. This capability is particularly important in investigations of chemical reactivity within complex catalytic systems, where experimental information may be scarce due to the transient nature of intermediates and the presence of numerous conformers with distinct thermochemical properties and implications for product applications.

Methods
Loss function design

Details regarding the Autoencoder architecture and the selection of hyperparameters for Skewencoder are provided in the SI. This section is devoted to the discussion of the formulation and implementation of the skewness-related loss term, L_skew (defined in Equation (4)). Consider the derivative of the additional skewness-related loss function term, L_skew, with respect to the NN parameters {θ}n := {fW_l,j g, f b_l g}n, which include all trainable variables: the weight elements {w_ij^l} from matrix W_l and the bias terms {b_i^l} from vector b_l across all layers l in the Skewencoder. The subscript n denotes that these parameters are optimized based on the n-th batch. The gradient of L_skew, as shown in Equation (5), appears in any gradient-based optimization method, ensuring the integration of skewness considerations into the training process.

∇_θn L_skew(x0; θ) = dL_skew / dγ^2 * ∇_θn γ(s(x0; θ))^2 (5)

We calculated the derivative based on γ^2 instead of γ itself, primarily because, within the scope of Skewencoder, only the magnitude of the skewness is considered, regardless of its sign. Given the form of the skewness loss function in Equation (4), the corresponding derivative with respect to γ^2 can be computed as shown in Equation (6).

dL_skew / dγ^2 = - e^(-γ^2) / (1 + e^(-γ^2)) = - 1 / (1 + e^(γ^2)) (6)

Since the absolute value of this term monotonically decreases as γ^2 increases, it can be regarded as a correction factor for the learning rate, fully based on skewness. When the skewness is small, the derivative value approaches -1/2, allowing for larger search steps during the early training epochs. Conversely, when the skewness is large, the derivative value approaches 0, contributing to the stable convergence of the loss function minimization. In this way, the derivative term in Equation (6) effectively acts as an exponentially decreasing learning rate, enhancing the performance of the standard optimization algorithm. The effectiveness and efficiency of incorporating L_skew into the loss function are quantitatively analyzed using benchmark model potential (see SI).

Biasing wall design

Figure 7 illustrates the procedure for defining the biasing wall in latent space, based on sampling from previous iterations of the Loxodynamics simulation. At each nth iteration, the multitask Autoencoder is retrained using the state-detection restart training protocol, incorporating information from all preceding n − 1 iterations. The resulting latent CV space is constructed, and samples from the local dataset (as defined in Figure 2) at iteration n−1 are projected onto this updated latent space. This projection enables the calculation of statistical moments—specifically, mean, variance, and skewness—for the sample distribution from iteration n−1. These statistical properties are then used to define a harmonic bias wall within this generalized coordinate system, as described in Equation (7).

Apply wall: κ(s_LS - (μ_LS + σ_LS + ϵ))^2 (7)

when s_LS <= sgn(γ_LS) * (μ_LS + σ_LS + ϵ)

where s_LS is the CV in the latent space, κ is the energy constant for the wall, and μ_LS, σ_LS and γ_LS are the sample mean, variance and sample skewness of the local dataset projection in the latent space. For positive skewness, a lower boundary—corresponding to a left-sided half-harmonic potential as depicted in Figure 7—is imposed. In contrast, negative skewness results in the application of an upper boundary. The parameter ϵ > 0 is a custom offset.

As depicted in Figure 7, the inclusion of both σ_LS and ϵ shifts the wall within the latent space, adaptively enhancing the asymmetry of the local distribution and accelerating the search in the vicinity of local minima on the PES. However, it is important to note that the offset terms σ_LS and ϵ also influence the resolution of the exploration, which can be inferred from Figure 7. Values of ϵ <= 1.0 are considered moderate.

Since the latent space is treated as a generalized coordinate with an implicit direction, it is crucial to determine whether the projected distribution is positively or negatively skewed and to set the wall accordingly.

General workflow

In this section, we present an overview of our methodology, summarized in the general workflow outlined in Box 1. The procedure is initialized with a small swarm of independent trajectories, each characterized by a different value of the biasing parameter κ_i. This strategy allows us to explore multiple directions efficiently and identify suitable paths for further exploration. The maximum number of allowed iterations, n_iter, is also specified to regulate computational cost. For each κ_i in the sequence, Loxodynamics begins by performing an unbiased simulation. During each iteration of the subsequent process, Skewencoder is trained using a state-detection-informed restart strategy (see SI for details). To demonstrate the enhanced efficiency of the warm start strategy, we compared Loxodynamics simulations utilizing this approach against exploration training initiated from scratch (see SI). The harmonic wall is configured as illustrated in Figure 7 and Equation (7), and steered MD simulations are executed accordingly.

For every value of κ_i, the Loxodynamics simulation terminates either upon reaching n_iter or upon identification of a new state. Loxodynamics outputs not only newly discovered states but also corresponding MD trajectories and latent space CVs from Skewencoder, which may be utilized for subsequent analysis of minimum energy pathways.

Simulations settings

Model potentials. We performed Langevin dynamics using the pesmd tool in PLUMED 2.9.072,73, defining the potential directly in the input file. All parameters are in Lennard-Jones units, with temperature 1.0, a Langevin thermostat friction coefficient of 10.0, and a timestep of 0.05 for 20,000 steps.

Gas phase reactions. Simulations were performed using CP2K 2023.174, patched with PLUMED 2.9.0, in the NVT ensemble with a 0.5 fs time step at 300 K, controlled by a canonical sampling through velocity rescaling (CSVR) thermostat75. Each sampling iteration consisted of 5000 steps in the biased MD simulation.

Energy and forces were computed using the PM676 semiempirical Hamiltonian with a Self-Consistent Field (SCF) convergence tolerance of 1.0 × 10^-5.

The system was placed within an orthorhombic simulation cell with no periodic boundary conditions, set to 15.0 Å for SN2 and 10.0 Å for Diels-Alder.

To detect new metastable states, bond topology changes were assessed by monitoring atomic contacts, using the function in Equation (8)77. The bond-dependent threshold p (0 < p < 0.5) distinguishes different states:

Bond formed, if c_i,j > 1 - p
Bond broken, if c_i,j < p
Not stable, otherwise

where c_i,j = 1 / (1 + (d_i,j / d_0)^m) (8)

where d_0 represents the typical bond length between atoms of species i and j. The use of machine-learning-based classifiers that could offer a more general or automated strategy.

While the implemented state detection method is straightforward and effective for the presented case studies, complex systems may necessitate more robust approaches. State-of-the-art machine learning classifiers, particularly those designed for transition state detection78, offer a generalized and automated strategy. These methods could provide the fast and precise state determination required for the warm-start training of Loxodynamics. Alternatively, for simpler organic fragments or structures that can be chemically isolated from their environment, SMILES-based fingerprinting79 serves as a valuable tool, allowing changes in bond topology to be rapidly identified via string analysis.

Catalytic ethanol dehydration in acidic chabazite. Simulations were conducted with CP2K 2023.1 patched with PLUMED 2.9.0 in the NVT ensemble, using a 0.5 fs timestep at 473 K, regulated by the CSVR thermostat. Each iteration comprised 5000 steps in both biased and unbiased MD simulations.

Energy and forces were computed with the semi-empirical GFN-XTB80 method, using an SCF convergence tolerance of 1.0 × 10^-5. The simulation box was the relaxed primitive rhombohedral unit cell of natural H-CHA, with lattice constants of 9.514 Å and angles of 94.07°.

Catalytic 1-butanol dehydration in acidic chabazite. Simulations were conducted with CP2K 2025.1 patched with PLUMED 2.9.0 in the NVT ensemble, using a 0.5 fs timestep at 473 K, regulated by the CSVR thermostat. Each iteration comprised 5000 steps in both biased and unbiased MD simulations.

Energy and forces were computed with the semi-empirical GFN-XTB method, using an SCF convergence tolerance of 1.0 × 10^-5. The simulation box was the relaxed primitive rhombohedral unit cell of natural H-CHA, with lattice constants of 9.514 Å and angles of 94.07°.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript presents Loxodynamics, a novel enhanced sampling scheme meant to explore chemical pathways in a\ndata-driven way. At the core of the method is a deep autoencoder architecture (Skewencoder) featuring a bespoke loss\nfunction which uses the skewness of the latent space distribution to identify chemical reaction coordinates (CV). This CV is\niteratively updated and used to bias the dynamics via a static half wall potential.\nThe problem of identifying chemical reaction pathways is extremely relevant and any development in this area is of interest\nto many researchers. The proposed method introduces some interesting new approaches, such as the clever combination of\ndeep learning and statistical analysis of the anharmonicity to guide the sampling. I believe it deserves publication, especially\nafter the following minor points have been addressed.\n1. Since it is important for the method to properly work, it should be mentioned in the main paper that one has to have a way\nto identify the metastable states. At least in the SI, please provide an example of how the metastable states were identified in\nthe systems presented.\n2. It would be interesting to report, in the main text or in the SI, the trajectory of the learned CV, i.e. of the latent space. This\nwould give a more intuitive idea of the anharmonicities the method is set up to follow. How does this CV compare with the\nones shown in the plots? How much does this CV change at each iteration?\n3. My understanding is that the skewness gamma of the latent space can switch sign between iterations, spontaneously\nbreaking the symmetry of the loss L_skew. This leads to the annoyance of having to keep track of the sign of gamma when\ndefining the biasing wall. But, instead of that, one could simply change the definition of the CV to ` s’_LS = sign(gamma_ls) *\ns_LS `, this would simplify eq. 8, fig. S1, and some of the code, since one would always need a lower wall.\n4. Could you briefly comment on the choice of eq. 4 for the skewness loss? Why not a simple gaussian? Why not an\nasymmetric function that would always push the skewness of the latent space to be of the same sign?\n5. Please be more explicit in the mathematical formulation on the loss function, eq. 2, and report it in full, e.g in the methods\nsection. This would help the reader to understand. The loss L should be a function of the parameters of the NN, theta, and of\nthe descriptors x. Gamma is computed on the latent space y, thus gamma(y) is also a function of the parameters theta.\n6. “A larger harmonic constant was applied when steering the reaction in the direction [...] to account for the higher barrier”\nhow do you choose the harmonic constant? Is it difficult?\n7. It’s not clear how Loxodynamics finds multiple paths. Trial and error? When should I stop? Does this mean that the first\npath I find could be a higher FES pathway, and not the most common transition state?", "reviewer_2": "The manuscript by Zhang and Piccini presents an original and promising approach to reaction discovery based on the\nskewness of probability distribution using an autoencoder framework. The idea behind the method—Loxodynamics—is\nclearly explained and introduces a creative use of probability distribution asymmetry to guide sampling. However, several\nmethodological and validation aspects require clarification and improvement before the work can be recommended for\npublication.\nMajor Comments\n>> Related Work and Literature.\nThe paper lacks a thorough discussion of prior literature on reaction discovery using machine learning collective variables,\nespecially autoencoder-based methods. Relevant prior work includes:\n- Chen & Ferguson (https://doi.org/10.1063/1.5023804) demonstrated the use of autoencoders for interleaved CV discovery\nand enhanced sampling, iteratively expanding the sampled space.\n- Ketkaew & Luber (https://doi.org/10.1021/acs.jcim.2c00883) applied autoencoders to generate blind CVs for chemical\nreactions starting only from reactants\n- Raucci et al. (https://doi.org/10.1021/acs.jpclett.1c03993) proposed a method for blind reaction discovery using general-\npurpose CVs with OPES explore.\n>> Method Validation\n- Comparison with standard autoencoders: what advantage does the “Skewencoder” offer over a standard autoencoder with\nthe same architecture? A comparative experiment (e.g., on the toy model) would help clarify its contribution.\n- Minimum Free Energy Pathways: How can it be guaranteed that the discovered paths correspond to minimum free energy\npaths? Except for the first toy model, the others only have one pathway. The last example has two pathways, and the fact\nthat it can also find a higher pathway seems to contradict this aspect. Also, it seems that this will depend on the\nhyperparameters (e.g., \\kappa, \\epsilon)? Is prior knowledge of energy barriers required to set them?\n- The discussion section claims a wide range of potential applications, but the validation is limited to a few relatively simple\ncases (except for the last one). A more realistic assessment of current limitations, along with suggestions for how to\ntransition from exploration to detailed mechanistic modeling, would strengthen the manuscript. For instance: after discovery,\ncan the CVs produced by Skewencoder be used in conventional biased simulations to compute free energy surfaces?\n>> Methodological Details\n- It is unclear whether \\kappa is chosen by the user or adapted during the iterations (as suggested. In the SN2 example, it is\nsaid: “Due to the asymmetric energy barriers… a larger harmonic constant was applied … to account for the higher barrier.”\nDoes setting this parameter require prior barrier knowledge?\n- The rationale behind the inclusion of \\sigma in the wall position (i.e., \\mu + \\sigma + \\epsilon) is not explained in the main\ntext. Figure 1 should also reflect this logic explicitly. Also, the role of \\epsilon seems important for determining the degree of\nexploration, but is not discussed. A discussion of how to choose these parameters and what effect they have is needed.\n- The claim that the reaction discovery is “unbiased” is misleading. Since the skewness of the biased distribution depends\non the applied parameters of the harmonic bias and no reweighting is performed, the method introduces a directional\nsampling artifact.\n- A complex procedure to decide whether to warm-start the autoencoder or not is put forward, but not justified. A comparison\nwith an exploration training only from scratch on the toy model would clarify its impact. Also, this procedure requires a\ncriterion to determine whether a new metastable state has been reached, which is not described for the examples reported.\nThis contradicts the claim of minimal user intervention and should be detailed more clearly.\n- Section A of the methods could benefit from a few revisions. I think that the workings of neural network layers are now part\nof common knowledge and so can be moved to the SI. On the other hand, it is not clear if the details about the L_{skew} are\njust conceptual considerations (as it seems to me looking at the code) or technical details that need to be implemented to\nimprove the performance of the method. Whichever the case, I think this should be made clear to the reader.\n- In general, it would be much better to include in the methods section a summary of the method, including the loss functions\n(which are only defined in the introduction) and the workflow of the iterative procedure (which is again discussed in the\nintroduction and summarized only in the SI).\nBelow there are other comments specific to the applications presented.\n>> Toy Model\n- The evolution of the skewness and the mean of the CV across iterations is not shown. Including this would help illustrate\nhow the method progresses and adapts the biasing direction based on the measured skewness?\n- How does the exploration change as a function of the chosen hyperparameters, \\kappa and \\epsilon? Showing\nsupplementary examples would help clarify how to set them.\n- Although the potential includes two distinct paths, only one is sampled. Can the method recover both, depending on the\nhyperparameters used? This would also help clarify how this is obtained in the chabazite case.\n>> SN2 Reaction\n- The input CV space essentially consists of d_{Cl} and d_F, which already capture the transition without needing any\ncombination of them. In my opinion, this system should be moved to the SI, as it does not significantly validate the method.\n>> Chabazite\n- How is the second pathway found? In the SI it is said that different elastic constants are used to reflect barrier differences.\nSince barrier information is usually unknown, how would one proceed in general?\n- In the stepwise mechanism, how is backward motion prevented? Are additional walls required? Would this become limiting\nfor more complex catalytic cycles involving many intermediate steps?\n>> Minor Comments\n- Many figure captions are not self-contained. For example, the meaning of color schemes and symbols (e.g., Figure 1\nsymbols, Figure 3 colorbar, color of scatter plots, Figure 4 color of scatter plots, Figure S1 symbols, Figure S2 colorbars,\nFigure S4 color of scatter plots, colorbars ) is only defined in the main text and should be reported in the captions.\n- The provided code is well-structured and detailed, which ensures reproducibility. However, from the SI and the code, it\nappears to rely extensively on the mlcolvar library, which is not cited in the main text.\n- When quantities are discussed in the main text, refer to the specific equation number to help guide the reader (e.g.,\nL_{skew} in methods sec. A).", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature\nCommunications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career\nResearchers who co-review manuscripts.", "reviewer_4": "In this manuscript, Zhang and Piccini propose a novel approach to explore reactive paths using an encoder that includes the\nskewness of the probability distribution observed in a (set of) collective variables into its loss function. I find the manuscript to\nbe clearly written and the method is well-explained. The authors, moreover, support their method with a handful of examples\nranging from toy potentials to heterogeneous catalysis. While I would most certainly recommend publication in a more\nspecialized journal, I believe that, in their present form, the results might not be sufficiently groundbreaking or of broad\nenough applicability given the multidisciplinary audience of Nature Communications. My main points of concern, in this\nregard, are the following:\n1. The authors state “Moreover, upon algorithm termination, the accumulated samples provide sufficient data to construct a\nlatent space that effectively approximates a meaningful CV”. While it looks like that is the case with the simple toy systems,\nit would be interesting to see if the conclusion holds when going to the more realistic reactions (e.g., plotting the isolines on\nthe d6 vs d4 graph in Figure 4). I would find a reaction finder that also provides an (ok-ish) CV to sample the reaction with\nenhanced sampling quite interesting.\n2. Concerning the last case study, I find it presently unclear what the difference was between the simulations that led to the\nexploration of the concerted vs stepwise mechanisms (aside from the wall kappa, as reported in the SI). Was sampling the\ntwo paths in the two simulations a coincidence? How different are they in energy according to the semiempirical method\nused by the authors? Also, in the stepwise mechanism, how likely it is to actually resample the backward reaction from the\nsurface species to ethene instead of going forward?\n3. I find that some comparison with literature is missing. For example, how does the method proposed by the authors\ncompare with OPES-based exploration methods (e.g., https://pubs.acs.org/doi/10.1021/acs.jpclett.1c03993) in terms of cost\nand ability of finding new reaction products?\n4. The reported case studies, while comprehensive, are all pretty simplistic. To remain in zeolite catalysis, I guess that using\na larger adsorbed molecule (with, e.g., double bonds) could open up quite a lot more reactive paths (cyclization, cracking,\nand so on). Going back to comment 2, it would be interesting to see which paths is the method able to discover.\nIf the authors can reply to these comments, the manuscript might become suitable for Nature Communications. I additionally\nhave the following minor comment:\n1. I am not fully sure I understand all the passages behind warm-start training explained in the SI. Am I correct in stating that\nwarm start is used unless the system falls into a new basin of attraction, in which case everything restarts from scratch? In\nthe final part, the authors mention the case in which the system falls back to the original minimum. I do not immediately see\nhow training from scratch would suddenly push the system towards another barrier."}, "author_response": "Response to the reviewer’s comments to manuscript\nNCOMMS-25-28852-T\n“Exploring Chemistry and Catalysis by Biasing Skewed\nDistributions via Deep Learning”\nZhikunZhang,GiovanniMariaPiccini\nOctober2,2025\nWewouldliketosincerelythankthereviewersfortheirthorough,thoughtful,andconstructivecriticism. The\ncommentswereceivedwerebothrigorousandfair,andwegreatlyappreciatedtheopportunitytocarefully\naddresseachofthem. Followingyoursuggestionshasbeenanintellectuallyrewardingprocess,andweare\nconfidentthatthemanuscripthasimprovedsubstantiallyinclarity,scope,andscientificimpactasaresult.\nWe are very glad to have incorporated all your requests, and we hope you will find the revised version\nsignificantlystrengthenedthankstoyourvaluableinput.\nReviewer #1\nRemarkstotheAuthor:\nThe manuscript presents Loxodynamics, a novel enhanced sampling scheme meant to explore chemical\npathwaysinadata-drivenway. Atthecoreofthemethodisadeepautoencoderarchitecture(Skewencoder)\nfeaturingabespokelossfunctionwhichusestheskewnessofthelatentspacedistributiontoidentifychemical\nreactioncoordinates(CV).ThisCVisiterativelyupdatedandusedtobiasthedynamicsviaastatichalfwall\npotential.\nTheproblemofidentifyingchemicalreactionpathwaysisextremelyrelevantandanydevelopmentinthis\nareaisofinteresttomanyresearchers. Theproposedmethodintroducessomeinterestingnewapproaches,\nsuchastheclevercombinationofdeeplearningandstatisticalanalysisoftheanharmonicitytoguidethe\nsampling. Ibelieveitdeservespublication,especiallyafterthefollowingminorpointshavebeenaddressed.\n1. Sinceitisimportantforthemethodtoproperlywork,itshouldbementionedinthemainpaperthat\nonehastohaveawaytoidentifythemetastablestates. AtleastintheSI,pleaseprovideanexampleof\nhowthemetastablestateswereidentifiedinthesystemspresented.\n#Response#:\nsubsection\"Metastablestatesidentification\"(Sec. S1. C,Page5,Line78,seeSI)detailingthemethodologyusedto\nidentifymetastablestatesinboththereactionsystemandthemodelsystemmentionedinourstudy. Additionally,we\nhaveexplicitlystatedthecriterioninthemaintext(Page11,Line759-760,seemaintext).\n2. Itwouldbeinterestingtoreport,inthemaintextorintheSI,thetrajectoryofthelearnedCV,i.e. ofthe\nlatentspace. Thiswouldgiveamoreintuitiveideaoftheanharmonicitiesthemethodissetuptofollow.\nHowdoesthisCVcomparewiththeonesshownintheplots? HowmuchdoesthisCVchangeateach\niteration?\n#Response#:\nWethankthereviewerforthisinsightfulcomment. Assuggested,wehaveaddedCVcontourplotsforrepresentative\nsimulationiterationsinthemodelsystemstotheSI(Figure. S4,Page8,seeSI),illustratinghowthelatentspaceCV\nevolvesovertime. Furthermore,wereportthattheLoxodynamicsbiasesenhancesamplingtowardsaddlepoints,thereby\ncapturinglocalanharmonicfeaturesmoreeffectively.\n3. Myunderstandingisthattheskewnessgammaofthelatentspacecanswitchsignbetweeniterations,\nspontaneouslybreakingthesymmetryofthelossL . Thisleadstotheannoyanceofhavingtokeep\nskew\ntrackofthesignofgammawhendefiningthebiasingwall. But,insteadofthat,onecouldsimplychange\nthedefinitionoftheCVtos′ =sign(γ )∗s ,thiswouldsimplifyeq. 8,fig. S1,andsomeofthecode,\nLS ls LS\nsinceonewouldalwaysneedalowerwall.\n#Response#:\nWedeeplyappreciatethiscommentfromthereviewerandweembracedthevaluablesuggestion. Inthisregard,the\nsimplifiedformatproposedbythereviewerismathematicallyequivalenttotheformulapresentedinEquation8ofthe\noriginalmanuscript. Tofurtheranalyzetheseaspects,weconductedcomparativeexperimentsonthreemodelPESto\nsubstantiatethisassertion,substitutingtheinitialwallapplicationwithanewlysimplifiedversionwhileconsistently\nimposingalowerboundarycondition. Theresultsareentirelyreproducibleusingthisrevisedharmonicwallapplication\nformat. Resultscanbefoundineithergithubcodeorcosinedatabaseofthecorrespondingmanuscript. Consequently,we\nhaveamendedthemaintext’sformulaconcerningwallapplicationaccordingtotherecommendedform(Equation7,Page\n10,Line687-688). Nonetheless,weintendtoretaintheoriginalformulationwithinthealgorithmdescriptionsectionof\ntheSI.Theterms\"lowerwall\"and\"upperwall\"possessgreaterphysicalsignificance,enhancingcomprehensionofhow\nskewnesscorrespondstolocalanharmonicityandhowLoxodynamicsbiasisappliedaccordingly(FigureS1,Page4,see\nSI).Thisapproachmayfacilitatereadersinrecognizingcoherenceandimprovetheirunderstanding.\n4. Couldyoubrieflycommentonthechoiceofeq. 4fortheskewnessloss? Whynotasimplegaussian?\nWhynotanasymmetricfunctionthatwouldalwayspushtheskewnessofthelatentspacetobeofthe\nsamesign?\n#Response#:\nThepointraisedbythereviewerisveryimportantandwearewillingtoclarifyitinourresponsehereaswellasby\naddingnewinformationintherevisedSI(videinfra).\nspacesuchthatprojectionsof\"localdataset\"sampleswithinthisspaceexhibitrelativelymaximalskewness. Therefore,\ntheskewlosstermshouldserveasametricforthemagnitude(i.e.,absolutevalue)ofsampleskewnessprojectedonto\nthelatentspace,excludingitssign. Uponcompletingoptimizationbasedonthisabsolutevalue,thesignofoptimized\nprojectedsampleskewnessrevealstheslowestascentdirectionwithinthelatentspace. AsnotedinPoint3byReviewer\n#1,withcertainmathematicalmanipulations,onecanflipthissign,implyinghowsystemprogressestowardsareaction\nchannelisreflectedentirelybythisabsolutevalueofsampleskewness. However,employinganasymmetricfunction\nthatconsistentlyforceslatentspaceskewnesstomaintainauniformsignaffectslossfunctionoptimizationandmay\nresultinunforeseenperformancedegradationduringexploration.\nNonetheless,varioussymmetricformsofskewness-relatedlossfunctionscantheoreticallyserveasvalidalternativesto\nthecurrentapplicationwithintheSkewencoderframework. Newtestsandanalyseshavebeenconductedondifferent\nmodelPESforcomparison,asdetailedintherevisedSI(SectionS2.A.1,Page8-9). Weevaluatedtheperformanceof\nSkewencoderusingthreedistinctformulationsoftheskewnesslossterm: 1)theoriginalimplementation;2)aGaussian-\nlike form, represented as exp(−γ2); and 3) a simplified quadratic form of projected sample skewness, expressed as\n−γ2. Considering criteria primarily focused on computational efficiency and training resource consumption, our\noriginalimplementationoftheskewnesslossfunctiondemonstratedsuperiorperformanceamongalltestedformats.\nThiscomparisonsubstantiatesourassertionthattheinitiallydesignedlossfunctionoffersadvantagesinfacilitating\nmoreefficienttrainingandsamplingprocesses.\n5. Pleasebemoreexplicitinthemathematicalformulationonthelossfunction,eq. 2,andreportitinfull,\ne.ginthemethodssection. Thiswouldhelpthereadertounderstand. ThelossLshouldbeafunctionof\ntheparametersoftheNN,theta,andofthedescriptorsx. γiscomputedonthelatentspacey,thusγ(y)is\nalsoafunctionoftheparameterstheta.\n#Response#:\nWeappreciatethereviewerforhighlightingthisissue. Inresponse,wehaveprovidedamorecomprehensiveexplanation\nofthelossfunction. Specifically,theautoencoderreconstructionlossL isdependentontheglobalinputdataset\nAutoencoder\nxanditsassociatedneuralnetworkparametersθ,whereasL iscontingentuponthelocalinputdatasetx′and\nskewness\ncorrespondingneuralnetworkparametersθ. Similarly,thelatentspace,akintoL ,isafunctionofthelocalinput\nskewness\ndatasetx′anditsrespectiveNNparametersθ. Asthelocalinputdatasetx′isinvariablyasubsetoftheglobaldataset\nx,thesetwosetswerenotdifferentiatedinournotationofthetotallossfunction. Consequently,wehaverevisedthe\nnotationaccordinglyinEquations2(Page4,Line236-237),3(Page4,Line252-253),4(Page4,Line263-264),and5\n(Page10,Line644-645)withinthemaintext.\n6. \"Alargerharmonicconstantwasappliedwhensteeringthereactioninthedirection[...] toaccountfor\nthehigherbarrier\"howdoyouchoosetheharmonicconstant? Isitdifficult?\n#Response#:\nThereviewerraisedanimportantpointregardingthesimulationsetup. Whilethisaspectisindeedrelevant,itdoesnot\nposesignificantdifficultiesinpractice. WerealizedthattheoriginalflowchartinFigureS1didnotaccuratelyreflectthe\nactualprotocolused.\nsmallswarmofparallelexplorations,eachwithadifferentκvalue,increasingwithinaspecificrange. Thisstrategy\nensuresmoreefficientsamplingofthereactionnetworkwhilekeepingcomputationalcostandtimeundercontrol.\nWehavenowupdatedFigureS1(Page4,seeSI)toreflectthisparallelexplorationscheme,whichbetterrepresentsthe\nmethodologywetypicallyemploy.\n7. It’snotclearhowLoxodynamicsfindsmultiplepaths. Trialanderror? WhenshouldIstop? Doesthis\nmeanthatthefirstpathIfindcouldbeahigherFESpathway,andnotthemostcommontransitionstate?\n#Response#:\nWeagreewiththereviewerthat,initscurrentform,Loxodynamicsidentifiesmultiplereactionpathwaysthrougha\nprocessthatpartiallyreliesontrialanderror. Atitscore, thisprocessisstronglyinfluencedbythechoiceoftheκ\nparameter,whichdefinesthestrengthofthebiasingwallandthusguidesthesystemtowardthelowest-lyingfreeenergy\npathwaysonthesurface. Nonetheless,theprocedureishighlyautomatedand,asdescribedinthereponsetocomment6\n\"thisstrategyensuresmoreefficientsamplingofthereactionnetworkwhilekeepingcomputationalcostandtimeunder\ncontrol\".\nToillustratethispoint,wehaveincludedaseriesofcomputationalexperimentsusingthetoymodelPESpresentedin\nthemaintext,whichfeaturetwoaccessiblereactionchannels(FigureS6,Page11,seeSI).Theseexampleshelpclarify\nhowvaryingκvaluescanaffecttheexplorationoutcome.\nLookingforward,weareactivelyworkingonintegratingareinforcementlearning(RL)frameworkintoLoxodynamics.\nInthisapproach,anagentwilliterativelylearntodiscriminateandselectthemostpromisingdirectionsamongaswarm\nofindependenttrajectories. Theobjectiveistoincreasethesuccessrateofdiscoveringlow-barriertransitionpathswhile\nmaintainingdiversitybydiscouragingconvergenceontothesamepathway.\nAlthoughstillinitsearlystages, thisRL-drivenextensionrepresentsakeydirectionforthefuturedevelopmentof\nLoxodynamicsandwillbethesubjectofforthcomingpublications.\nRemarksoncodeavailability:\nThecodeisclearlywrittenandreasonablydocumented. Itcontainsallthenecessarytoolstoreproducethe\nresultsofthepaperorusetheproposedmethodondifferentsystems,includingaJupyternotebooktutorial.\nReviewer #2\nRemarkstotheAuthor:\nThemanuscriptbyZhangandPiccinipresentsanoriginalandpromisingapproachtoreactiondiscovery\nbasedontheskewnessofprobabilitydistributionusinganautoencoderframework. Theideabehindthe\nmethod—Loxodynamics—is clearly explained and introduces a creative use of probability distribution\n4\nasymmetrytoguidesampling. However,severalmethodologicalandvalidationaspectsrequireclarification\nandimprovementbeforetheworkcanberecommendedforpublication.\nMajorComments\n>>RelatedWorkandLiterature.\nThepaperlacksathoroughdiscussionofpriorliteratureonreactiondiscoveryusingmachinelearning\ncollectivevariables,especiallyautoencoder-basedmethods. Relevantpriorworkincludes: Weappreciate\nthereviewer’scarefulattentioninhighlightingfundamentalliteraturethatshouldbecited. Below,we\noutlinetheadditionswehaveincorporatedintothemanuscripttoaddressthesesuggestions.\n-Chen&Fergusonhttps://doi.org/10.1063/1.5023804demonstratedtheuseofautoencodersfor\ninterleavedCVdiscoveryandenhancedsampling,iterativelyexpandingthesampledspace.\n#Response#:\nWeaddedthesuggestedreferencetoChenandFergusonJ.Chem. Phys. onewealreadycitedonthesametopic(Page3,\nLine208-219,seemaintext).\n-Ketkaew&Luberhttps://doi.org/10.1021/acs.jcim.2c00883appliedautoencoderstogener-\nateblindCVsforchemicalreactionsstartingonlyfromreactants.\n#Response#:\nWeaddedthesuggestedreference(Page3,Line219-231,seemaintext).\n-Rauccietal. https://doi.org/10.1021/acs.jpclett.1c03993proposedamethodforblindre-\nactiondiscoveryusinggeneral-purposeCVswithOPESexplore.\n#Response#:\nWeaddedthesuggestedreference(Page2,Line90-91,seemaintext).\nNonetheless, we believe that a further clarification on the distinction between Loxodynamics and other schemes is\nwarranted. TheaforementionedapproachesemployCVstogetherwithenhancedsamplingmethodssuchasumbrella\nsampling,metadynamics,orOPEStoacceleratesampling. Thesemethods,however,donotencodeaspecificdirection-\nality. Incontrast,Loxodynamicsintroducesabiasexplicitlydesignedtodrivethesystemtowardanewstatealong\ntheCV.MetadynamicsandOPES,inparticular,are\"flooding\"approaches,inwhichthebiasgraduallyfillsthelocal\nfree-energybasintoenhancefluctuationswithinit. ThephilosophyofLoxodynamicsisfundamentallydifferent: rather\nthanamplifyingfluctuations,itreshapesthemwiththespecificpurposeofpropellingthesystemuphill.\n>>MethodValidation\n-Comparisonwithstandardautoencoders: whatadvantagedoesthe\"Skewencoder\"offeroverastandard\n5\nautoencoderwiththesamearchitecture? Acomparativeexperiment(e.g.,onthetoymodel)wouldhelp\nclarifyitscontribution.\n#Response#:\nWe thank the reviewer for raising this important point. Indeed, in our initial attempts, we employed a standard\nAutoencodertoachieveourgoal. Thisinvolvedprojectingtheconfigurationsontothelatentspaceandevaluatingthe\nskewnessoftheresultingdistributionwithoutexplicitlyoptimizingforit. Inprinciple,suchanapproachcouldwork.\nAccordingly,weperformedcomparativeexperimentsinvolvingboththecurrentSkewencoderconfigurationandapure\nAutoencoderformat(SectionS2.A.1,Page10-11,seeSI),alongsideevaluationsofvariousformulationsoftheskewness\nlossterm(refertoReviewer#1’sPoint4).\nHowever,weconsistentlyobservedthat,especiallyintheearlystagesofexplorationwhendataarescarce,theskewness\nsignal tends to be weak and significantly affected by noise. As a consequence, the direction selected by the biasing\nwalls—based on this unreliable signal—amounts to an effectively random search (see Figure S5, Page 10, see SI).\nThebiasthenfailstoconsistentlysteerthesystemalongmeaningfulreactionpaths, greatlyreducingthemethod’s\neffectiveness.\nTheselimitationsmotivatedthedesignoftheSkewencoderarchitecture,wherethelatentspaceisexplicitlyoptimizedto\nenhancetheskewnesssignal,therebyprovidingamorerobustanddirectionaldrivingforceforexploration.\n-MinimumFreeEnergyPathways: Howcanitbeguaranteedthatthediscoveredpathscorrespondto\nminimumfreeenergypaths? Exceptforthefirsttoymodel,theothersonlyhaveonepathway. Thelast\nexamplehastwopathways,andthefactthatitcanalsofindahigherpathwayseemstocontradictthis\naspect. Also,itseemsthatthiswilldependonthehyperparameters(e.g.,κ,ϵ)? Ispriorknowledgeof\nenergybarriersrequiredtosetthem?\n#Response#:\nWeappreciatethereviewer’scommentandwegladlyintroducedfurtherinformationinthesubmittedrevisedversion.\nFirstly,wemustclarifythatκisnotahyperparameteroftheSkewencoderbutratheravariablefortheexploration\nprocess. AsindicatedintherevisedFigureS1inSI(Page4,seeSI)(alsoreferencedinReviewer#1’sPoint6),noprior\nknowledgeofκisnecessary. InthenewlyupdatedversionofSectionIV.C(Page11,Line712-731,seemaintext),we\nspecifythatinsteadofsequentiallyincreasingκaftereachexplorationcampaign,ourmethoddeploysasmallswarm\nofparallelexplorations,eachwithdistinctκvalues. Thisapproachfacilitatesmoreefficientsamplingofthereaction\nnetworkwhilemaintainingcontrolovercomputationalcostandtime. Thus,κisinitializedasaparameterforagnostic\nexplorationofLoxodynamicstothoroughlyandefficientlysearchthechemicalreactionspace.\nConcerning the influence of κ, we invite the reviewer to refer to Point 7 from Reviewer #1, where we detailed our\nevaluationofexplorationperformancebyinitializingvariousκvalueswithintheLoxodynamicssimulationappliedto\nthedouble-wellmodelPESillustratedinthemaintext. Briefly,ahigherκtendstoidentifyreactionpathwayswith\nelevatedreactionbarriers,whereasasuitablylowerκisinclinedtodiscoverpathwayswithreducedreactionbarriers,\ncloselyapproximatingtheanticipatedMEP.Thus,selectingκcansignificantlyaffecttheexplorationofdistinctreaction\n6\npathways.\nIntheoriginalmanuscript,wehadalreadydiscussedtheroleofbothσandϵindeterminingthepositionoftheharmonic\nwall within the latent space (Page 10, Line 699-705, see main text). In the revised version, we have expanded the\nmethodologysectionandincludedanewfigure(Figure7,Page10,seemaintext). Thisfigureillustratesthatϵactsasa\ntunableoffset,effectivelycontrollingthedegreeof\"bruteforce\"appliedduringLoxodynamicssimulationstoexplore\nnewregionsofthelandscape.\nImportantly,byadjustingthenewly-addedwallbasedonprojectedsamplevariance,thiswallcanexploreunknown\nchemicalspacesregardlessofϵ’svalue. Providedthatϵisgreaterthan0,itwillnotconstraintherangeofextended\nchemicalspaceexploration.Therefore,Loxodynamics’performanceisnotparticularlysensitivetoϵ’sselection,especially\nwhenitissmallandpositive. Sinceallsamplingdataareprojectedintothelatentspaceandnormalized—i.e.,explored\n(known) space latent CVs are normalized to [-1, 1]—selecting ϵ at 0.1 represents a balanced choice that enhances\nsimulationefficiencywhilepreservingcontinuousandsmoothexplorationprocesses. Thechosenϵremainsconsistent\nacrossalltestcases,ensuringitsneutralityconcerningourmethod. Inotherwords,controllinginsensitivevariables\nenablesustoconcentrateonanalyzingκ’seffectonLoxodynamicssimulationperformance.\nAdditionally,toascertainwhethertheMDtrajectorygeneratedfromLoxodynamicssimulationapproximatestheMEP\noftheidentifiedreactionpathway,weutilizedtheClimbing-ImageNudgedElasticBand(CINEB)methodtooptimize\nourMDtrajectorysnippetdepictingreactionsinbothconcertedandstepwisemechanisms(SectionS2.C.3,Page18,\nseeSI).Insummary,thegeometriesoftheMDtrajectoryfromLoxodynamicssimulationcloselyresembleoptimized\nconfigurationalgeometriesonNEBpaths,suggestingthattheMDtrajectoryisalreadyneartheMEPofdiscovered\nreactionsandindicatingthatLoxodynamicscanguidesystemsapproximatelytowardscalculatedMEPdirectionsfor\ncertainreactions. FurthersimulationdetailsareavailableinSectionS2.C.3(Page18,seeSI).\n-Thediscussionsectionclaimsawiderangeofpotentialapplications,butthevalidationislimitedtoa\nfewrelativelysimplecases(exceptforthelastone). Amorerealisticassessmentofcurrentlimitations,\nalongwithsuggestionsforhowtotransitionfromexplorationtodetailedmechanisticmodeling,would\nstrengthenthemanuscript. Forinstance: afterdiscovery,cantheCVsproducedbySkewencoderbeused\ninconventionalbiasedsimulationstocomputefreeenergysurfaces?\n#Response#:\nWeappreciatethereviewerforthisinsightfulcomment. InthenewlyaddedSectionS2.B.3(Page15-16,seeSI),we\npresentresultsfromwell-temperedMetadynamicssimulationsemployingthelearnedSkewencoderCVinbothS 2and\nN\nDiels-Alderreactions. TheseSkewencoderCVsarederivedfromiterationsofLoxodynamicssimulationswherestate\ntransitionsoccurred. ByproperlyconfiguringMetadynamicsparameters,FigureS11(Page16,seeSI)demonstrates\nthattheSkewencoderCVcanfacilitatetransitionsbetweentwostatesinMetadynamicssimulationsofbothreactions.\nTheSkewencoderCValsoeffectivelydifferentiatesreactantandproductstates,correspondingtothestateinformation\nobtainablefromheuristicCVs. FigureS12(Page16,seeSI)illustratesthatsamplesfromwell-temperedMetadynamics\nusing Skewencoder CVs can reconstruct the FES for specific reactions. The reconstructed FES aligns with those\ngeneratedbyMetadynamicssimulationsemployingheuristicCVs. Consequently,SkewencoderCVscanbefurther\nappliedinconventionalbiasedsimulations,providedtheyencapsulateconfigurationalinformationwithinbothreactant\nandproductspaces.\nhasbeenexploredinanagnosticmannerusingLoxodynamics. Pleasecheckouranswertocomment4toReviewer#4.\n>>MethodologicalDetails\n-Itisunclearwhetherκischosenbytheuseroradaptedduringtheiterations(assuggested. IntheSN2\nexample,itissaid: \"Duetotheasymmetricenergybarriers... alargerharmonicconstantwasapplied...\ntoaccountforthehigherbarrier.\"Doessettingthisparameterrequirepriorbarrierknowledge?\n#Response#:\nThereviewerraisedanimportantpointregardingthesimulationsetup. Whilethisaspectisindeedrelevant,itdoesnot\nposesignificantdifficultiesinpractice. WerealizedthattheoriginalflowchartinFigureS1didnotaccuratelyreflectthe\nactualprotocolused.\nsmallswarmofparallelexplorations,eachwithadifferentκvalue,increasingwithinaspecificrange. Thisstrategy\nensuresmoreefficientsamplingofthereactionnetworkwhilekeepingcomputationalcostandtimeundercontrol.\nWehavenowupdatedFigureS1(Page4,seeSI)toreflectthisparallelexplorationscheme,whichbetterrepresentsthe\nmethodologywetypicallyemploy.\nIngeneral,noprecisepriorknowledgeofthebarrierisrequired.\n-Therationalebehindtheinclusionofσinthewallposition(i.e.,µ+σ+ϵ)isnotexplainedinthemain\ntext. Figure1shouldalsoreflectthislogicexplicitly. Also,theroleofϵseemsimportantfordetermining\nthedegreeofexploration,butisnotdiscussed. Adiscussionofhowtochoosetheseparametersandwhat\neffecttheyhaveisneeded.\n#Response#:\nWeappreciatethereviewerforhighlightingthisissue.Inthenewlyrevisedmanuscript,weprovideadetailedexplanation\ninthemaintextregardingtheimpactofutilizingprojectedsamplevarianceσandacustomoffsetϵ,asillustratedin\nFigure7(Page10,seemaintext),ratherthanmerelyportrayingtheminFigure1asbefore.\nInthecurrentrevisedversion,Figure1(Page3,seemaintext)primarilyillustratesthattherelationshipbetweensample\nskewnessandlocalanharmonicityarisesfrombiasedsimulationsonanartificial1-DPES.Duringeachiterationof\nthe1-Dsimulation,abiasharmonicwallisintroducedandconsistentlyshiftedwithaconstantoffset∆xfromtheleft\ntotherightalongthex-axis. Weobservedthatasthesystemtransitionedfromoneminimumtoanother,thesample\ndistributionbecameincreasinglyasymmetric,reflectingitsgradualevolutiontowardsstatetransition. Inthecurrent\nmanuscriptversion,wehaveclarifiedthatthisobservationfromthe1-Dsimulationservedasinspirationforgoing\nfromasimpleidealcasetomorerealisticscenarios. Consequently,toaddresschallengesinextendingthisknowledgeto\nmultidimensionalscenarios,wedevelopedtheSkewencoderNNarchitectureandLoxodynamicsworkflowasinnovative\ntoolsforchemicalspaceexploration.\nThe impact of the custom offset ϵ has been addressed in the response to the Method Validation section of Major\nCommentsfromReviewer#2. Similarly,weemployprojectedsamplevariancetodeterminetheharmonicwall’sposition.\nouralgorithmincursnoadditionalcomputationalcost. Furthermore,samplevariancereflectsthedistributionwidth\nwithinlatentspace. Byusingthistermasanoffsetfortheharmonicwall,previouslyexploredregionsinlatentspace\narebypassedinnewsimulationiterations,reducingoverlapandenhancingexplorationefficiency. Inpractice,onemay\nreplaceprojectedsamplevarianceσorevenσcombinedwithϵwithaconstantoffset(asillustratedinFigure1,page\n3,seemaintext). ThissubstitutionyieldsvaryingeffectsacrossdifferentPES.Byintroducingσandamoderateϵ,\nwedevisedamechanismthatadaptivelyadjustsLoxodynamics’exploration\"speed\"basedonlocaldatasetswithout\nincurringextracomputationalexpense.\n-Theclaimthatthereactiondiscoveryis“unbiased”ismisleading. Sincetheskewnessofthebiased\ndistributiondependsontheappliedparametersoftheharmonicbiasandnoreweightingisperformed,\nthemethodintroducesadirectionalsamplingartifact.\n#Response#:\nWe appreciate the reviewer for highlighting this issue. We would like to clarify that we do not claim the reaction\ndiscoveryprocesstobeunbiased,butratheragnostic. Althoughitistruethatnoreweightingisperformed,such\nanactionwouldprimarilyreproducetheequilibriumdistributionoftheinitialstate,providinglimitedinsightinto\ntransitionpathways. Therefore,acertaindegreeofdirectionalitymustbeintroduced. Importantly,thisdirectionalbias\nisnotimposedaprioribytheuserbutemergesduringtheearlystagesofthesimulationwhenthesystemexperiences\nminimalornobias. Asexplorationprogresses,themethodselectsspecificdirectionsforsamplingreactionpathways.\nThisselectionisnotmanuallyadjustedbutdrivenbyaniterativesamplingprotocolthatdynamicallyadaptstothe\nunderlyingfreeenergylandscape. Forclarity,wehaveremoved\"unbiased\"fromthesentence: Byfollowingpathways\nwiththelowestenergybarriers,thisapproachenablesefficientandunbiasedexplorationoflikelyreactionchannels\n(Page2,Line101,seemaintext).\n-Acomplexproceduretodecidewhethertowarm-starttheautoencoderornotisputforward,butnot\njustified. Acomparisonwithanexplorationtrainingonlyfromscratchonthetoymodelwouldclarifyits\nimpact. Also,thisprocedurerequiresacriteriontodeterminewhetheranewmetastablestatehasbeen\nreached,whichisnotdescribedfortheexamplesreported. Thiscontradictstheclaimofminimaluser\ninterventionandshouldbedetailedmoreclearly.\n#Response#:\nWeappreciatethereviewerforraisingthisinsightfulissue.\nIntherevisedversionofourmanuscript,wehaveincludedSectionS1.C(Page5,seeSI)toelucidatethecriterionfor\ndistinguishingbetweenmetastablestatesandnon-metastablestates. Generally,thiscriterioncanbeadaptedtoany\nconnectivitymatrix-basedmethodcapableofdeterminingatomicconnectivitythroughcustomizedswitchingfunctions.\nThepresentimplementationisanin-housesimplifiedversiondesignedforefficiency.\nAdditionally,acomparativeexperimentwasconductedbetweentheLoxodynamicssimulationutilizingthewarmstart\ntrainingstrategyandexplorationtraininginitiatedfromscratchonthetoymodelPES.Theresultsarepresentedin\nSectionS2.A.3(Page11,seeSI).Insummary,thewarmstarttrainingstrategyenhancesexplorationefficiencyby\nloweringtrainingandcomputationalcosts,therebyfulfillingtheintendedobjectiveofthisstrategywithinLoxodynamics\n9\nsimulation.\n-SectionAofthemethodscouldbenefitfromafewrevisions. Ithinkthattheworkingsofneuralnetwork\nlayersarenowpartofcommonknowledgeandsocanbemovedtotheSI.\n#Response#:\nWethankthereviewerforthisvaluablesuggestion. ThiscontenthasbeenrelocatedtoSectionS1.B.2(Page3,Line\n43-51,inSI),integratedwiththediscussionontheNNstructureofSkewencoder,inthenewlyrevisedversionofthe\nmanuscript.\nOn the other hand, it is not clear if the details about the L are just conceptual considerations (as\nskew\nit seems to me looking at the code) or technical details that need to be implemented to improve the\nperformanceofthemethod. Whicheverthecase,Ithinkthisshouldbemadecleartothereader.\n#Response#:\nWe appreciate the reviewer for identifying this issue. The current format of L encompasses both conceptual\nskew\nconsiderationsandtechnicalimplementation,enhancingthemethod’sperformance. RefertoPoint4fromReviewer#1\nandthediscussioninthenewlyaddedSectionS2.A.1(Page8-9,inSI)regardingthismatter. Insummary,thecurrent\nimplementationofL hasbeenshowntoimproveexplorationefficiencywhencomparedwithotherformsofL .\nskew skew\nAdditionally,SkewencoderutilizingL inapurelyGaussian-likeformsuccessfullyidentifiedalternativestatesand\nskew\nexhibitedonlyslightlyreducedperformanceduringLoxodynamicssimulationonthetoyPES,suggestingthatL\nskew\ncanbeappliedflexiblyandisthereforealsoaconceptualconsideration.\n-Ingeneral,itwouldbemuchbettertoincludeinthemethodssectionasummaryofthemethod,including\nthelossfunctions(whichareonlydefinedintheintroduction)andtheworkflowoftheiterativeprocedure\n(whichisagaindiscussedintheintroductionandsummarizedonlyintheSI).\n#Response#:\nWeappreciatethereviewerforthishelpfulsuggestion. Wehavereorganizedthemethodologysectioninthemaintext\naccordingly. Toeliminateredundancy,wecontinuetofocusondiscussingthedesignoftheskewnessloss-relatedterm\ninSectionIV.A(Page9-10,seemaintext),buthaveremovedtheintroductiontoageneralNNstructure. Subsequently,\nwepartiallyrevisedSectionVI.B(Page10,seemaintext). WeincludedFigure7(Page10,seemaintext),whichmore\nclearlyillustrateshowtheharmonicwallisappliedbasedonsamplemean,variance,andskewnessvalues. Additionally,\nweemployedasimplifiedequationformindicatinghowtheharmonicwallisappliedinlatentspace,assuggestedby\nReviewer#1inPoint3. Asaconclusiontothemethodologicalsection,wepresentasimplifiedversionoftheentire\nalgorithminSectionIV.C(Page11,seemaintext),correspondingtodetaileddiscussionsinSectionS1oftheSI.\nCommentsspecifictotheapplicationspresented.\n>>ToyModel\n-TheevolutionoftheskewnessandthemeanoftheCVacrossiterationsisnotshown. Includingthis\nwouldhelpillustratehowthemethodprogressesandadaptsthebiasingdirectionbasedonthemeasured\nskewness?\n10\n#Response#:\nWeappreciatethereviewerforhighlightingthisissue. Wewishtoclarifythatthelatentspaceundergoeschangesat\neachiteration,renderingcomparisonsofdistributionmomentsacrossdifferentvariablesimpractical.\nConcerningCVevolutionasrequestedbyReviewer#1,pleaserefertoPoint1fromReviewer#1whereweillustratethe\nevolutionofCVsintoymodels.\n- How does the exploration change as a function of the chosen hyperparameters, κ and ϵ? Showing\nsupplementaryexampleswouldhelpclarifyhowtosetthem.\n#Response#:\nWeappreciatethereviewerforthiscomment. Wehavepreviouslyaddressedtheimpactofϵinearlierquestions(e.g.,\nMethodValidationsectionofMajorCommentsfromReviewer#2). Regardingtheeffectofκ,asdiscussedinprevious\nfeedback,briefly,ahigherκtendstoidentifyreactionpathwayswithelevatedreactionbarriers,whereasasuitablylower\nκisinclinedtodiscoverpathwayswithreducedreactionbarriers,closelyapproximatingtheanticipatedMEP.Thus,\nselectingκcansignificantlyinfluencetheexplorationofdistinctreactionpathways. WehaveaddedSectionS2.A.2\n(Page9-11,seeSI)toassesstheeffectofκusingatoymodelinLoxodynamics.\nConcerningSkewencoderhyperparameters,weassertthatwithatypicaldeepautoencoderneuralnetworkstructure,\nSkewencoder can effectively learn features from latent space and drive Loxodynamics. We have already compared\nSkewencoderperformanceusingbothdeeperandshallowerneuralnetworkstructuresinSectionS2.B.2(Page13-15,see\nSI).Bothconfigurationsoftheneuralnetworksuccessfullydrivethesystemtoidentifystatetransitions.Therefore,given\ntherobustfittingcapabilityofneuralnetworks,ageneralAutoencoderstructurecanbeemployedwithinSkewencoder\nwithoutextensivehyperparameteroptimizationefforts. Thishighlightsourmethod’sflexibilityandscalability.\n- Although the potential includes two distinct paths, only one is sampled. Can the method recover\nboth,dependingonthehyperparametersused? Thiswouldalsohelpclarifyhowthisisobtainedinthe\nchabazitecase.\n#Response#:\nWeappreciatethereviewerforhighlightingthisissue. WekindlyrequestthereviewertorefertoPoint7fromReviewer\n#1andSectionS2.A.2(Page9-11,seeSI).\n>>SN2Reaction\n- The input CV space essentially consists of d and d , which already capture the transition without\nCl F\nneedinganycombinationofthem. Inmyopinion,thissystemshouldbemovedtotheSI,asitdoesnot\nsignificantlyvalidatethemethod.\n#Response#:\nWeappreciatethereviewerforthissuggestion. However,wewouldliketoemphasizethatalthoughtheS 2reaction\nN\nisarelativelysimpleexample,itservesasarealistictestcase. Giventhattheinputdimensionsofthepreviousmodel\nPESareall2D,andtheS 2reactioninvolvestwoSkewencoderinputdescriptors,itcanbeconsideredanappropriate\nN\nintermediarybetweenmorecomplexreal-worldexamplesandsimplermodelPES,whichareeasiertotest. Consequently,\nwehavedecidedtoretainthissectioninthemaintext,asitcontinuestovalidatethemethod.\n>>Chabazite\n-Howisthesecondpathwayfound? IntheSIitissaidthatdifferentelasticconstantsareusedtoreflect\nbarrierdifferences. Sincebarrierinformationisusuallyunknown,howwouldoneproceedingeneral?\n#Response#:\nThisiscertainlyafundamentalpoint. AsdiscussedinourresponsetoPoint6fromReviewer#1andillustratedin\nFigureS1(Page4,seeSI),ourmethodinitiatesasmallswarmofparallelexplorations,eachwithadistinctκvaluethat\nincreaseswithinaspecifiedrange. ThisapproachallowsLoxodynamicssimulationtoidentifycorrespondingreaction\npathwayswithdifferentκvalueswithoutrequiringbarrierinformation. IntheSI,wementionthatvariousconstants\nareemployedtoreflectbarrierdifferences,indicatingthatdifferentκvaluesconfiguringharmonicwallscorrespond\ntobarrierheightdifferencescalculatedusingDFTmethods(asreportedinhttps://pubs.acs.org/doi/10.\n1021/jp513024z). Wedonotimplythatκdirectlycorrespondstoenergyvalues;rather,smallerκvaluesarelikely\ntouncoverreactionpathwayswithlowerbarriers,whereaslargerκvaluestendtorevealpathwayswithhigherbarriers.\nThisobservationalignswithourfindingsinthedoublepathwaymodelPESpresentedinSectionS2.A.2(Page9-11,see\nSI).\n-Inthestepwisemechanism,howisbackwardmotionprevented? Areadditionalwallsrequired? Would\nthisbecomeliitingformorecomplexcatalyticcyclesinvolvingmanyintermediatesteps?\n#Response#:\nWeappreciatethiscommentverymuchaswerealizethispointwasnotclearintheoriginalsubmission. Thus,we\nupdatedaccordinglySectionS2.C.2(Page18,seeSI),discussinghowwerestarttheLoxodynamicssimulationtoexplore\nareactionpathwaythatincludesmultiplereactionwells,e.g. thestepwisemechanismintheethanoldehydrationcase. In\nbrief,FigureS13(Page17,seeSI)showsthatoncewereachedtheintermediatestate,werestarttheLoxodynamcis,keep\nalltheparametersrelatedtotheSkewencoder,butweretaintheoldCValongwhichweapplyanadditionalrestraining\nwall that guide the system to identify the current intermediate state. Then in the newly restarted Loxodynamics\nsimulation,an“unbiased”simulationwillbeconducted,withtherestrainingwallkept. ThennewLoxodynamicsbias\nwillbeaddedinthefollowingiteration,withtherestrainingwallstillexisting. Tocompensatetheadditionalpotential\nbiascausedbytherestrainingwall,κvaluemustbeincreasedtoexplorethesecondpotentialwell,thoughthebarrier\nheight of the first reaction and the second reaction in the multiple-well reaction path are similar. We reported this\nobservationinSectionS2.C.2(Page18,seeSI)andS2.D(Page18-20,seeSI),correspondingtodifferentdehydration\nreactionsontheacidicChabasite. Basedonourobservation,weassumethatthismightbealimitationtoareaction\npathwaywithmorewells. However,inthetestcaseswehaveconductedLoxodynamicssimulationsonsofarwedonot\nobservesuchalimitation.\n>>MinorComments\n-Manyfigurecaptionsarenotself-contained. Forexample,themeaningofcolorschemesandsymbols\n(e.g.,Figure1symbols,Figure3colorbar,colorofscatterplots,Figure4colorofscatterplots,FigureS1\nsymbols,FigureS2colorbars,FigureS4colorofscatterplots,colorbars)isonlydefinedinthemaintext\nandshouldbereportedinthecaptions.\n#Response#:\nWeappreciatethereviewerforidentifyingthisissue. Below,wedetailthecontentaddedtothecaptionsoftherelevant\nfigurestoaddressthisconcern. Generally,weemphasizethattheprecisevaluesonthecolorbar(Figure3andFigure\nS2)forboththemodelPESandprojectedlatentspace(FigureS4)arenotcritical. Theclosedformathasalreadybeen\nprovidedforthePES.AsLoxodynamicssimulationsprogress,theprojectedCVspaceevolvescontinuously,withexact\nvaluesinprojectedlatentspaceCVbeingupdatediteratively. Consequently,relativerelationshipsaremoresignificant\nthanexactvaluesforclearlyillustratingCVevolution.\nInthecaptionofFigure1(Page3,seemaintext),weadded:\n-Theinitialwallislocatedatx=µ,correspondingtothemeanoftheunbiasedsampledistribution(indicatedbyred\ndotsintheupperpanel). Ineachsubsequentiteration,theharmonicwallistranslatedbyafixedoffset∆xrelativetoits\npositioninthepreviousiteration. Thespringconstantκforallwalls(depictedasblue,green,andyellowhalf-parabolic\ncurvesinthelowerpanel)remainsconstantthroughoutalliterations.\nInthecaptionofFigure3(Page5,seemaintext),weadded:\n-Thedarkerregionsignifiesanareawithlowerpotentialvalues,whereasthelighterregiondenotesanareawith\nhigherpotentialvalues.\n- The central and lower panels display the evolution of the x (in dark green) and y (in light green) coordinates,\nrespectively,overtime(inLennard-Jonesunits). Thebluepointsrepresentsamplesfromtheincrementaldataset,while\nredpointscorrespondtosamplesfromthecurrentiteration.\nInthecaptionofFigure4(Page6,seemaintext),weadded:\n-Inthe2Dmaps,asinFigure3,bluepointsdenotesamplesfromtheglobaldataset,whereasredpointsindicate\nsamplesfromthelocaldatasetcorrespondingtothecurrentiteration. Relevantstructuresextractedfromthetrajectory\nblocksofeachsuccessiveiterationareshownintheuppersubpanels.\nInthecaptionofFigureS1(Page4,seeSI),weadded:\n- In this flowchart, {κ } represents a sequence of wall spring constants arranged in an arithmetic progression\ni\nwithacommondifferenceofδκ. Theparametern denotesthemaximumnumberofiterationspermittedforthe\niter\nLoxodynamicssimulation. Thequantityγ referstotheprojectedsampleskewnessinlatentspace. ThevariableL\nls state\nindicatesthelistofidentifiedmetastablestates,whileS servesastheidentifierforthecurrentstate.\np\nInthecaptionofFigureS3(Page7,seeSI),weadded:\n-Thedarkerregionsignifiesanareawithlowerpotentialvalues,whereasthelighterregiondenotesanareawith\nhigherpotentialvalues.\nInthecaptionofFigureS4(Page8,seeSI),weadded:\n13\n-Ineachleftsubcolumn,Thebluepointsrepresentsamplesfromtheincrementaldataset,whileredpointscorrespond\ntosamplesfromthecurrentiteration.\n-LighterregionscorrespondtohighervaluesinthelatentCVspace,whereasdarkerregionsrepresentlowervalues.\n-Theprovidedcodeiswell-structuredanddetailed,whichensuresreproducibility. However,fromtheSI\nandthecode,itappearstorelyextensivelyonthemlcolvarlibrary,whichisnotcitedinthemaintext.\n#Response#:\nWe appreciate the reviewer for this suggestion. The relevant citation has been incorporated in the code and data\navailabilitysection(Page12,seemaintext).\n-Whenquantitiesarediscussedinthemaintext,refertothespecificequationnumbertohelpguidethe\nreader(e.g.,L inmethodssec. A).\nskew\n#Response#:\nWeappreciatethereviewerforthissuggestion. Wehaveaddedthefollowinginthemaintext(Page9-11,line631-633):\nThissectionisdevotedtothediscussionoftheformulationandimplementationoftheskewness-relatedlossterm,L\nskew\n(definedinEquation4),toaddressthisissue.\nRemarksoncodeavailability:\nTheprovidedcodeiswell-structuredanddetailed,whichensuresreproducibility. However,fromtheSIand\nthecode,itappearstorelyextensivelyonthemlcolvarlibrary,whichisnotcitedinthemaintext.\nReviewer #3\nRemarkstotheAuthor:\nIco-reviewedthismanuscriptwithoneofthereviewerswhoprovidedthelistedreports. Thisispartofthe\nNatureCommunicationsinitiativetofacilitatetraininginpeerreviewandtoprovideappropriaterecognition\nforEarlyCareerResearcherswhoco-reviewmanuscripts.\nReviewer #4\nRemarkstotheAuthor:\nInthismanuscript,ZhangandPicciniproposeanovelapproachtoexplorereactivepathsusinganencoder\nthatincludestheskewnessoftheprobabilitydistributionobservedina(setof)collectivevariablesintoits\nlossfunction. Ifindthemanuscripttobeclearlywrittenandthemethodiswell-explained. Theauthors,\nmoreover,supporttheirmethodwithahandfulofexamplesrangingfromtoypotentialstoheterogeneous\ncatalysis. WhileIwouldmostcertainlyrecommendpublicationinamorespecializedjournal,Ibelievethat,\nintheirpresentform,theresultsmightnotbesufficientlygroundbreakingorofbroadenoughapplicability\n14\ngiventhemultidisciplinaryaudienceofNatureCommunications. Mymainpointsofconcern,inthis\nregard,arethefollowing:\n1. Theauthorsstate“Moreover,uponalgorithmtermination,theaccumulatedsamplesprovidesufficient\ndatatoconstructalatentspacethateffectivelyapproximatesameaningfulCV”. Whileitlookslikethatis\nthecasewiththesimpletoysystems,itwouldbeinterestingtoseeiftheconclusionholdswhengoingto\nthemorerealisticreactions(e.g.,plottingtheisolinesonthed vsd graphinFigure4). Iwouldfinda\n6 4\nreactionfinderthatalsoprovidesan(ok-ish)CVtosamplethereactionwithenhancedsamplingquite\ninteresting.\n#Response#:\nWe appreciate the reviewer for highlighting this issue. As pointed out also in our response to Reviewer #2, in the\nnewlyaddedSectionS2.B.3(Page15-16,seeSI),wepresentresultsfromwell-temperedMetadynamicssimulations\nemployingthelearnedSkewencoderCVinbothS 2andDiels-Alderreactions. TheseSkewencoderCVsarederived\nN\nfromiterationsofLoxodynamicssimulationswherestatetransitionsoccurred. ByproperlyconfiguringMetadynamics\nparameters,FigureS11(Page16,seeSI)demonstratesthattheSkewencoderCVcanfacilitatetransitionsbetween\ntwostatesinMetadynamicssimulationsofbothreactions. TheSkewencoderCValsoeffectivelydifferentiatesreactant\nandproductstates,correspondingtothestateinformationobtainablefromheuristicCVs. FigureS12(Page16,seeSI)\nillustratesthatsamplesfromwell-temperedMetadynamicsusingSkewencoderCVscanreconstructtheFESforspecific\nreactions. ThereconstructedFESalignswiththosegeneratedbyMetadynamicssimulationsemployingheuristicCVs.\nConsequently,SkewencoderCVscanbefurtherappliedinconventionalbiasedsimulations,providedtheyencapsulate\nconfigurationalinformationwithinbothreactantandproductspaces.\n2. Concerningthelastcasestudy,Ifinditpresentlyunclearwhatthedifferencewasbetweenthesimula-\ntionsthatledtotheexplorationoftheconcertedvsstepwisemechanisms(asidefromthewallkappa,as\nreportedintheSI).Wassamplingthetwopathsinthetwosimulationsacoincidence?\n#Response#:\nWeappreciatethereviewerforthiscommentwhichhasalsobeenraisedbyotherReviewers.\nThesamplingoftwopathsisnotcoincidentalbutrathercontingentuponκasaparameterinitializedwithinaswarm.\nAsdiscussedinourresponsetoPoint6fromReviewer#1andillustratedinFigureS1(Page4,seeSI),ourmethod\ninitiatesasmallswarmofparallelexplorations,eachwithadistinctκvaluethatincreaseswithinaspecifiedrange. This\napproachallowsLoxodynamicssimulationtoidentifycorrespondingreactionpathwayswithdifferentκvalueswithout\nrequiringbarrierinformation. IntheSI,wementionthatvariousconstantsareemployedtoreflectbarrierdifferences,\nindicatingthatdifferentκvaluesconfiguringharmonicwallscorrespondtobarrierheightdifferencescalculatedusing\nDFTmethods(asreportedinhttps://pubs.acs.org/doi/10.1021/jp513024z). Wedonotimplythat\nκdirectlycorrespondstoenergyvalues;rather,smallerκvaluesarelikely\ntouncoverreactionpathwayswithlowerbarriers,whereaslargerκvaluestendtorevealpathwayswithhigherbarriers.\nThisobservationalignswithourfindingsinthedoublepathwaymodelPESpresentedinSectionS2.A.2(Page9-11,see\nSI).\nHowdifferentaretheyinenergyaccordingtothesemiempiricalmethodusedbytheauthors?\n15\n#Response#:\nWe appreciate the reviewer for highlighting this insightful issue. To evaluate how closely the trajectory obtained\nfrom Loxodynamics simulation approximates the MEP of a specific reaction, particularly ethanol dehydration at\nacidic Chabazite, we employed the Climbing-Image Nudged Elastic Band (CI-NEB) method to optimize our MD\ntrajectorysnippetdepictingreactionsinbothconcertedandstepwisemechanisms. WeintroducedarefinedCI-NEB\nworkflowinSectionS2.C.3(Page18,seeSI),newlyadded. ThisworkflowbeginswithconductingCI-NEBcalculations\nusingarelativelylongMDtrajectorysnippetasinput,followedbyselectingrepresentativeimagesfromtheinitial\ncalculationresultstoserveasinitialimagesforasecondCI-NEBcalculation. ThisapproachallowstherefinedCI-\nNEB workflow to utilize information from the Loxodynamics simulation trajectory while reducing computational\ncosts. TheNEBenergiescalculatedthroughthismethodarecomparabletoDFTcalculationresults(asreportedin\nhttps://pubs.acs.org/doi/10.1021/jp513024z). Furthermore, the geometries of the MD trajectory\ncloselyresembleoptimizedconfigurationalgeometriesonNEBpaths,indicatingthattheMDtrajectoryusedasinput\nisalreadyneartheMEPofdiscoveredreactionsandsuggestingthatLoxodynamicscanguidesystemsapproximately\ntowardscalculatedMEPdirectionsforcertainreactions. FurthersimulationdetailscanbefoundinSectionS2.C.3\n(Page18,seeSI)\nAlso,inthestepwisemechanism,howlikelyitistoactuallyresamplethebackwardreactionfromthe\nsurfacespeciestoetheneinsteadofgoingforward?\n#Response#:\nWeappreciatethereviewerforhighlightingthisissue.\nAs pointed out also in our response to Reviewer #2, we realize this point was not clear in the original submission.\nThus,weupdatedaccordinglySectionS2.C.2(Page18,seeSI),discussinghowwerestarttheLoxodynamicssimulation\nto explore a reaction pathway that includes multiple reaction wells, e.g. the stepwise mechanism in the ethanol\ndehydrationcase. Inbrief,FigureS13(Page17,seeSI)showsthatoncewereachedtheintermediatestate,werestartthe\nLoxodynamcis,keepalltheparametersrelatedtotheSkewencoder,butweretaintheoldCValongwhichweapplyan\nadditionalrestrainingwallthatguidethesystemtoidentifythecurrentintermediatestate. Theninthenewlyrestarted\nLoxodynamics simulation, an “unbiased” simulation will be conducted, with the restraining wall kept. Then new\nLoxodynamicsbiaswillbeaddedinthefollowingiteration,withtherestrainingwallstillexisting. Tocompensatethe\nadditionalpotentialbiascausedbytherestrainingwall,κvaluemustbeincreasedtoexplorethesecondpotentialwell,\nthoughthebarrierheightofthefirstreactionandthesecondreactioninthemultiple-wellreactionpatharesimilar. We\nreportedthisobservationinSectionS2.C.2(Page18,seeSI)andS2.D(Page18-20,seeSI),correspondingtodifferent\ndehydrationreactionsontheacidicChabasite. Basedonourobservation,weassumethatthismightbealimitationtoa\nreactionpathwaywithmorewells. However,inthetestcaseswehaveconductedLoxodynamicssimulationsonsofar\nwedonotobservesuchalimitation.\n3. Ifindthatsomecomparisonwithliteratureismissing. Forexample,howdoesthemethodproposed\nbytheauthorscomparewithOPES-basedexplorationmethods(e.g.,https://pubs.acs.org/doi/10.\n1021/acs.jpclett.1c03993)intermsofcostandabilityoffindingnewreactionproducts?\n#Response#:\nWeappreciatethereviewerforhighlightingthispoint. Wehavenowincludedasentenceinthemaintext(Page2,Line\n16\n89-90,seemaintext)explicitlycomparingLoxodynamicstoOPES_E:\"Ingeneral,designingafullyagnosticCVfor\nreactivityexplorationwithhighsamplingefficiencyisimpracticalwithoutpriorknowledgeoftheexpectedproducts.\"\nWhileOPES_Ecanidentifynewmetastablestates,itstilldependsonaninitiallydefinedCV,necessaryforapplying\nbias. Additionally, OPES_E—similar to other flooding techniques like metadynamics or variationally enhanced\nsampling—operates by gradually filling free energy wells along the predefined CV, enabling transitions to other\nmetastablestates.\nConversely,ourmethodresemblesdirectionalpullingapproachessuchassteeredMDorratchet-and-pawlMD,where\nbiasisappliedinaspecificdirection. TheprimaryadvantageofLoxodynamicsisthatboththeCVandpullingdirection\narenotpredeterminedbytheuserbutaredynamicallyandautomaticallyupdatedduringexploration. Thisfacilitates\nmoreagnosticandadaptivesamplingofreactionpathways.\n4. Thereportedcasestudies,whilecomprehensive,areallprettysimplistic. Toremaininzeolitecatalysis,\nIguessthatusingalargeradsorbedmolecule(with,e.g.,doublebonds)couldopenupquitealotmore\nreactivepaths(cyclization,cracking,andsoon). Goingbacktocomment2,itwouldbeinterestingtosee\nwhichpathsisthemethodabletodiscover.\n#Response#:\nWeappreciatethereviewerfortheinsightfulsuggestion. Therefore,weembracedthesuggestionandstartedamuch\nmorecomplexcase,withalargerreactionnetwork. Moreprecisely,inthenewlyrevisedmanuscript,weintroducea\nmorecomplexexampleinzeolitecatalysis,namely1-ButanolDehydrationinacidicchabazite.\nThe1-Butanolcasesharessimilarcharacteristicswithethanoldehydration,consistingoftworeactionmechanisms—concerted\nand stepwise pathways. A key distinction between 1-butanol and ethanol is the diverse array of potential reaction\nproducts,influencedbydoublebondpositioningalongthealkylchainandassociatedisomerism. Literatureindicates\nthat1-butanoldehydrationcanyieldvariousproductsandmaybeaccompaniedbysidereactionslikeoligomerization,\ncracking,andaromatizationunderhighsubstrateloadingconditions. Intherevisedmanuscript,wefocusonsingle-site\nadsorptionspecificallyexaminingdehydrationpathways. AdditionaldetailsareprovidedinSectionII.D(Page8-9,see\nmaintext).\nFigure6(Page8,seemaintext)depictsreactionpathwaysidentifiedviaLoxodynamicssimulations. Alongthedirect\ndehydration route, the catalystproducesthreecanonical products of1-butanol dehydration: 1-buteneand cis- and\ntrans-isomersof2-butene. Asexpected,1-buteneformationmayalsoproceedthroughstepwisemechanismsakinto\nthoseinethanoldehydrationpathways. FurtherdetailedanalysisisavailableinSectionII.D(Page8-9,see\nmaintext).\nTheκvaluesreportedinS2.D(Page18-20,seeSI)correspondtomagnituderelationshipsamongdifferentproducts,\nparticularlyconformers.\nThiscasestudydemonstratesourmethod’scapabilitytoefficientlyexplorealargercatalyticreactionnetworkwhile\nidentifyingbothlikelyproductsandpathwaysinvolvingintermediateswithcomparableoccurrencelikelihoods. We\nbelievethiscasestudysufficientlyillustratesLoxodynamics’potentialapplicationforexploringanextendedrealmof\nreactionspaces.\nIftheauthorscanreplytothesecomments,themanuscriptmightbecomesuitableforNatureCommuni-\n17\ncations.\nIadditionallyhavethefollowingminorcomment:\n1. IamnotfullysureIunderstandallthepassagesbehindwarm-starttrainingexplainedintheSI.AmI\ncorrectinstatingthatwarmstartisusedunlessthesystemfallsintoanewbasinofattraction,inwhich\ncaseeverythingrestartsfromscratch? Inthefinalpart,theauthorsmentionthecaseinwhichthesystem\nfallsbacktotheoriginalminimum. Idonotimmediatelyseehowtrainingfromscratchwouldsuddenly\npushthesystemtowardsanotherbarrier.\n#Response#:\nWeappreciatethereviewerforhighlightingthisissue. Wewouldliketoclarifythatthewarmstarttrainingstrategyis\nemployedunless: 1)thesystementersanewbasinofattraction,or2)thesystemrevertstotheoriginalreactantspace\nwithoutidentifyinganotherstate. Inthelatterscenario,thesystemshouldtransitionfrom\"unstable\"statesbacktoits\noriginalstate. Definitionsregardingwhetherastateis(meta)stableorunstableareprovidedinSectionIV.D.2(Page11,\nseemaintext)andthenewlyaddedSectionS1.C(Page5,seeSI).\nInthescenariomentionedbythereviewerinSectionS2.B.2(Page13-15,seeSI),thesystemwasabruptlyrevertedto\nitsoriginalstatebecausetheκvaluewassufficientlylargetodrivethesystembackandforth,inabsenceofadditional\nrestrainingwallstopreventthisbehavior.ThisobservationalignswiththeshapeoftheFESobtainedfromMetadynamics\nsimulations using either heuristic CV or Skewencoder CV for Diels-Alder reaction (Figure S12, Page 16, see SI).\nFigure S9 (Page 14, see SI) illustrates that a large κ forcefully pushes the system from a wide but higher basin to\na lower and narrower basin, compared with trajectories generated by Loxodynamics simulation featuring deeper\nSkewencoderstructures(Figure4, Page6, seemaintext). Despitetraininginthenewstatestartingentirelyfrom\nscratch,Skewencodercanstillcaptureinformationfromalreadysampledconfigurationalspace.Moreover,Loxodynamics\nharmonicwallbiasisadaptivelypositionedusingsamplemeanandvariancevalues(Figure7,Page10,seemaintext),\npreventingredundantexplorationofknownconfigurationalspace. Undertheseconditions,thesystemcanswiftlyexit\nshallowbasinsandreturntooriginalstates. Insummary,hereweaimtodemonstratethatSkewencoder’sstructureis\nflexibleandcapableoftrainingrepresentativelatentspaceCVs,underscoringourworkflow’sconceptualrobustness. We\nacknowledgepotentialmisconceptionscausedbythisexample;henceweanalyzewarm-starttrainingstrategyeffectson\nmodelPESinSectionS2.A.3(Page11-12,seeSI).\n18"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors addressed all my concerns and the new version of the paper is in excellent shape, especially the\nmuch improved SI. I believe it makes a significant contribution to the field and I recommend it for publication.", "reviewer_2": "See attached file", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature\nCommunications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career\nResearchers who co-review manuscripts.", "reviewer_4": "I appreciate the effort that the authors have undertaken to reply to all the reviewers' comments. I find the additional\nexamples and information included in the manuscript and SI to drastically improve the reader’s understanding of the method as well as\nits relevance. Therefore, I would be happy to recommend the manuscript for publication in Nature Communications.\nI only have one final remark for the authors: in section S3, a comparison between the free energy profiles obtained from a\nchemically-inspired CV and the CV from Loxodynamics is reported (Figure S12). However, the free energy profile is CV\ndependent and comparing two profiles directly is not per se meaningful (see, e.g., https://doi.org/10.1016/j.jcat.2020.04.015\nor https://doi.org/10.1063/5.0083423). I would find it interesting to evaluate the actual CV-independent free energy difference\nbetween the two metastable states from the simulations results, to showcase that the Loxodynamics CV is (likely) fine and\ndiscrepancies in the free energy profile with chemically-inspired CVs are meaningless."}, "author_response": "Response to the reviewer’s comments to manuscript\nNCOMMS-25-28852A\n“Exploring Chemistry and Catalysis by Biasing Skewed\nDistributions via Deep Learning”\nZhikunZhang,GiovanniMariaPiccini\nDecember22,2025\nWesincerelythankthereviewersfortheirpositiveassessmentofourrevisedmanuscriptandtheiracknowl-\nedgement of the improvements made in the previous round. We appreciate the additional suggestions\nprovided,whichhavebeenveryhelpfulinfurtherpolishingthetextandclarifyingthefinaldetails. Wehave\naddressedtheseremainingpointsbelowandbelievethemanuscriptisnowfullyreadyforpublication. All\nadditionalsuggestionsfromthereviewersaremarkedinred,whileourcorrespondingresponsesareinblue\nitalics.\nReviewer #1\nRemarkstotheAuthor:\nTheauthorsaddressedallmyconcernsandthenewversionofthepaperisinexcellentshape,especiallythe\nmuchimprovedSI.IbelieveitmakesasignificantcontributiontothefieldandIrecommenditforpublication.\nReviewer #2\nRemarkstotheAuthor:\nIappreciatetheAuthors’effortsinaddressingthereviewers’questions,whichhaveclarifiedseveralimportant\naspectsofthework. Thatsaid,Ibelievethatsomeminoraspectscouldstillberefinedtofurtherimprovethe\nmanuscript.\nInparticular,inseveralcasestheclarificationsprovidedintheresponseletterhavenotbeenincorporated\nintothemanuscript,orhavebeenincludedonlyintheSupportingInformationwithoutbeingmentioned\ninthemaintext. Thismakesitdifficultforreaderstoaccessorappreciatethesedetails. Sincesomeof\ntheseissueswereraisedindependentlybymultiplereviewers,addingthecorrespondingexplanationsor\ncross-referencesinthemanuscriptwouldhelpmaketheworkclearerandmorerobust.\n1\n#Response#:\nWethankthereviewerforthesegeneralcomments. Wehaveincorporatedthecorrespondingexplanationsfromthe\noriginalresponseletterintotherevisedmanuscript. Foradetaileddiscussion,pleaserefertoourresponsestothespecific\nremarksbelow.\nFurthermore,theinclusionofthemorecomplexexample(1-butanoldehydration)isawelcomeimprove-\nmentandrepresentsavaluableextensionofthestudy. However,thisexamplelackssufficientsupporting\ndata:Inadditiontotheschematicofthe“discovered”species,somequantitativeevidence(e.g.,time-series\ndataorscatterplots)shouldbereportedintheSI.\n#Response#:\nWethankthereviewerforthiscomment. WehaveincludedFiguresS15andS16intheSupportingInformation(SI),\nwhichillustratetheevolutionofthesystemusingrepresentativebonddistances,analogoustotheapproachinFigure5.\nThecorrespondingtrajectoriesareavailableinthedatasetpublishedonZenodo.\nRegardingSectionD.2,itwouldalsobeusefultobrieflymentionalternativeorcomplementaryapproaches\nforstateidentification,assuggestedbyotherreviewers;forinstance,theuseofmachine-learning-based\nclassifiersthatcouldofferamoregeneralorautomatedstrategy.\n#Response#:\nWe thank the reviewer for this insightful suggestion. We have added the following content to Section D.2 of the\nmanuscripttoprovideacompleteelaborationofthestatedetectionmethod.\nWhiletheimplementedstatedetectionmethodisstraightforwardandeffectiveforthepresentedcasestudies,complex\nsystems may necessitate more robust approaches. State-of-the-art machine learning classifiers, particularly those\ndesignedfortransitionstatedetection,offerageneralizedandautomatedstrategy. Thesemethodscouldprovidethe\nfastandprecisestatedeterminationrequiredforthewarm-starttrainingofLoxodynamics. Alternatively,forsimpler\norganicfragmentsorstructuresthatcanbechemicallyisolatedfromtheirenvironment,SMILES-basedfingerprinting\nservesasavaluabletool,allowingchangesinbondtopologytoberapidlyidentifiedviastringanalysis.\nAdditionalPointsRelatedtotheOriginalRemarks\nMajorComments\n>>RelatedWorkandLiterature.\nThepaperlacksathoroughdiscussionofpriorliteratureonreactiondiscoveryusingmachinelearning\ncollectivevariables,especiallyautoencoder-basedmethods. Relevantpriorworkincludes: Weappreciate\nthereviewer’scarefulattentioninhighlightingfundamentalliteraturethatshouldbecited. Below,we\noutlinetheadditionswehaveincorporatedintothemanuscripttoaddressthesesuggestions.\n-Chen&Fergusonhttps://doi.org/10.1063/1.5023804demonstratedtheuseofautoencodersfor\ninterleavedCVdiscoveryandenhancedsampling,iterativelyexpandingthesampledspace.\n#Response#:\nWeaddedthesuggestedreferencetoChenandFergusonJ.Chem. Phys. onewealreadycitedonthesametopic(Page3,\nLine208-219,seemaintext).\n-Ketkaew&Luberhttps://doi.org/10.1021/acs.jcim.2c00883appliedautoencoderstogener-\nateblindCVsforchemicalreactionsstartingonlyfromreactants.\n#Response#:\nWeaddedthesuggestedreference(Page3,Line219-231,seemaintext).\nReviewer: Thewaythisreferenceismentionedisnotcontextualizedcorrectly;itseemstoevenimplythat\nitalreadysolvedtheproblemaddressedinthepresentworkbytheAuthors,whichcontrastswiththe\nclarificationgivenbelow.\n#Response#:\nWeappreciatethereviewerforhighlightingthisissue. Wehavereformulatedthediscussioninthemanuscripttoclarify\nthattheprimaryadvantageofLoxodynamicsisitsabilitytoactivelydriveexplorationtowardnewregionsofphase\nspaceviaanexplicitdirectionalbias. Notably,ourmethodachievesthiswithoutrelyingonconventionalenhanced\nsamplingparadigmsorstaticCVsderivedfromextensivepriorsimulations.\n\"Ketkaewetal. introducedanautoencoder-basedstrategyaimedatconstructingdata-drivencollectivevariablesfor\nenhancedsampling,trainedexclusivelyonreactant-stateconfigurations. Theirapproachexploitsdeepautoencoders\ntoextractlatentrepresentationsfromunbiasedtrajectoriesandemploystheserepresentationswithinflooding-based\nenhancedsamplingschemestoreconstructfreeenergysurfaces. Whilethismethodologyalleviatestheneedforexplicit\nproduct-stateinformationinthedefinitionofcollectivevariables,itremainsinherentlytiedtostaticCVsandrelieson\nconventionalenhancedsamplingparadigms,withoutintroducinganexplicitdirectionalbiascapableofactivelydriving\nexplorationtowardnewregionsofphasespace.\"\n-Rauccietal. https://doi.org/10.1021/acs.jpclett.1c03993proposedamethodforblindre-\nactiondiscoveryusinggeneral-purposeCVswithOPESexplore.\n#Response#:\nWeaddedthesuggestedreference(Page2,Line90-91,seemaintext).\nissuetheAuthorsdiscuss,itshouldbeatleastmentionedinthepartaboutsimilarapproaches;together\nwiththeclarificationbelow.\n#Response#:\nWe appreciate the reviewer for raising this point. Raucci et al. proposed an exploration strategy based on Social\nPermutationInvariantCoordinates(SPRINT)combinedwiththeOPESmethod. However,theirsingle-dimensional\nCV—thelargesteigenvalueoftheconnectivitymatrix—dependsstrictlyontheinitialgeometry. Consequently,the\nsearchdirectioniseffectivelypredeterminedatthestart. Incontrast,Loxodynamicsexploreschemicalspacedynamically;\nbyvaryingthewallconstantκ,differentproductspacescanbeaccessed. Furthermore,duetoitsdistinctsampling\npattern,Loxodynamicsdemonstratessuperiorefficiency,drivingthesystemtonewstateswithfivetotentimesless\nsimulationtime.\nRegardingthereactioncoordinatedefinition,Rauccietal. employDeepLDAtorefinetheCVdescribingthereactant\nandproductonlyafteranewstateisdetected. Conversely,LoxodynamicsyieldsaCVduringtheexplorationphase\nitself,whichisoftenimmediatelysuitableforsubsequentenhancedsamplingmethodslikeMetadynamics.\nWeinitiallyomittedadetaileddiscussionofthispaperbecauseitisnotanautoencoder-basedmethod. Instead,itfalls\nunderthecategoryofmethodsutilizingpredefinedCVswithconventionalenhancedsampling—approacheswhich,\naswenoted,oftenlacktheflexibilitytosamplediversedirections. However,toprovideacompletecontext,wehave\nincorporatedthefollowingexplanationintothemaintextandaddedthecitationforthereader’sreference.\nNonetheless, we believe that a further clarification on the distinction between Loxodynamics and other schemes is\nwarranted. TheaforementionedapproachesemployCVstogetherwithenhancedsamplingmethodssuchasumbrella\nsampling,metadynamics,orOPEStoacceleratesampling. Thesemethods,however,donotencodeaspecificdirection-\nality. Incontrast,Loxodynamicsintroducesabiasexplicitlydesignedtodrivethesystemtowardanewstatealong\ntheCV.MetadynamicsandOPES,inparticular,are\"flooding\"approaches,inwhichthebiasgraduallyfillsthelocal\nfree-energybasintoenhancefluctuationswithinit. ThephilosophyofLoxodynamicsisfundamentallydifferent: rather\nthanamplifyingfluctuations,itreshapesthemwiththespecificpurposeofpropellingthesystemuphill.\nReviewer: This answer clarifies the difference with respect to the other methods. Adding it to the\nmanuscriptwouldbehighlybeneficial.\n#Response#:\nWethankthereviewerforthissuggestion. Asnotedinthepreviouspoint,wehaveincorporatedtherelevantdiscussion\nintothemaintext.\n>>MethodValidation\n-Comparisonwithstandardautoencoders: whatadvantagedoesthe\"Skewencoder\"offeroverastandard\n4\nautoencoderwiththesamearchitecture? Acomparativeexperiment(e.g.,onthetoymodel)wouldhelp\nclarifyitscontribution.\n#Response#:\nWe thank the reviewer for raising this important point. Indeed, in our initial attempts, we employed a standard\nAutoencodertoachieveourgoal. Thisinvolvedprojectingtheconfigurationsontothelatentspaceandevaluatingthe\nskewnessoftheresultingdistributionwithoutexplicitlyoptimizingforit. Inprinciple,suchanapproachcouldwork.\nAccordingly,weperformedcomparativeexperimentsinvolvingboththecurrentSkewencoderconfigurationandapure\nAutoencoderformat(SectionS2.A.1,Page10-11,seeSI),alongsideevaluationsofvariousformulationsoftheskewness\nlossterm(refertoReviewer#1’sPoint4).\nHowever,weconsistentlyobservedthat,especiallyintheearlystagesofexplorationwhendataarescarce,theskewness\nsignal tends to be weak and significantly affected by noise. As a consequence, the direction selected by the biasing\nwalls—based on this unreliable signal—amounts to an effectively random search (see Figure S5, Page 10, see SI).\nThebiasthenfailstoconsistentlysteerthesystemalongmeaningfulreactionpaths, greatlyreducingthemethod’s\neffectiveness.\nTheselimitationsmotivatedthedesignoftheSkewencoderarchitecture,wherethelatentspaceisexplicitlyoptimizedto\nenhancetheskewnesssignal,therebyprovidingamorerobustanddirectionaldrivingforceforexploration.\nReviewer: Thanksforthisadditionalinvestigation. Ibelievethattheseresultsshouldbealsomentioned\ninthemaintext,sincetheyclarifytheneedtousetheproposedarchitectureasopposedtoaplainautoen-\ncoder.\n#Response#:\nWethankthereviewerforthissuggestion. WehaveaddedthefollowingsentencestoSectionVI.Aofthemaintextto\ndirectthereadertothecorrespondingdiscussionintheSI.\n\"TheeffectivenessandefficiencyofincorporatingL intothelossfunctionarequantitativelyanalyzedusingbench-\nskew\nmarkmodelpotential(seeSI).\"\n-MinimumFreeEnergyPathways: Howcanitbeguaranteedthatthediscoveredpathscorrespondto\nminimumfreeenergypaths? Exceptforthefirsttoymodel,theothersonlyhaveonepathway. Thelast\nexamplehastwopathways,andthefactthatitcanalsofindahigherpathwayseemstocontradictthis\naspect. Also,itseemsthatthiswilldependonthehyperparameters(e.g.,κ,ϵ)? Ispriorknowledgeof\nenergybarriersrequiredtosetthem?\n#Response#:\nWeappreciatethereviewer’scommentandwegladlyintroducedfurtherinformationinthesubmittedrevisedversion.\nFirstly,wemustclarifythatκisnotahyperparameteroftheSkewencoderbutratheravariablefortheexploration\nprocess. AsindicatedintherevisedFigureS1inSI(Page4,seeSI)(alsoreferencedinReviewer#1’sPoint6),noprior\nknowledgeofκisnecessary. InthenewlyupdatedversionofSectionIV.C(Page11,Line712-731,seemaintext),we\nspecifythatinsteadofsequentiallyincreasingκaftereachexplorationcampaign,ourmethoddeploysasmallswarm\nofparallelexplorations,eachwithdistinctκvalues. Thisapproachfacilitatesmoreefficientsamplingofthereaction\nnetworkwhilemaintainingcontrolovercomputationalcostandtime. Thus,κisinitializedasaparameterforagnostic\nexplorationofLoxodynamicstothoroughlyandefficientlysearchthechemicalreactionspace.\nConcerning the influence of κ, we invite the reviewer to refer to Point 7 from Reviewer #1, where we detailed our\nevaluationofexplorationperformancebyinitializingvariousκvalueswithintheLoxodynamicssimulationappliedto\nthedouble-wellmodelPESillustratedinthemaintext. Briefly,ahigherκtendstoidentifyreactionpathwayswith\nelevatedreactionbarriers,whereasasuitablylowerκisinclinedtodiscoverpathwayswithreducedreactionbarriers,\ncloselyapproximatingtheanticipatedMEP.Thus,selectingκcansignificantlyaffecttheexplorationofdistinctreaction\n6\npathways.\nIntheoriginalmanuscript,wehadalreadydiscussedtheroleofbothσandϵindeterminingthepositionoftheharmonic\nwall within the latent space (Page 10, Line 699-705, see main text). In the revised version, we have expanded the\nmethodologysectionandincludedanewfigure(Figure7,Page10,seemaintext). Thisfigureillustratesthatϵactsasa\ntunableoffset,effectivelycontrollingthedegreeof\"bruteforce\"appliedduringLoxodynamicssimulationstoexplore\nnewregionsofthelandscape.\nImportantly,byadjustingthenewly-addedwallbasedonprojectedsamplevariance,thiswallcanexploreunknown\nchemicalspacesregardlessofϵ’svalue. Providedthatϵisgreaterthan0,itwillnotconstraintherangeofextended\nchemicalspaceexploration.Therefore,Loxodynamics’performanceisnotparticularlysensitivetoϵ’sselection,especially\nwhenitissmallandpositive. Sinceallsamplingdataareprojectedintothelatentspaceandnormalized—i.e.,explored\n(known) space latent CVs are normalized to [-1, 1]—selecting ϵ at 0.1 represents a balanced choice that enhances\nsimulationefficiencywhilepreservingcontinuousandsmoothexplorationprocesses. Thechosenϵremainsconsistent\nacrossalltestcases,ensuringitsneutralityconcerningourmethod. Inotherwords,controllinginsensitivevariables\nenablesustoconcentrateonanalyzingκ’seffectonLoxodynamicssimulationperformance.\nAdditionally,toascertainwhethertheMDtrajectorygeneratedfromLoxodynamicssimulationapproximatestheMEP\noftheidentifiedreactionpathway,weutilizedtheClimbing-ImageNudgedElasticBand(CINEB)methodtooptimize\nourMDtrajectorysnippetdepictingreactionsinbothconcertedandstepwisemechanisms(SectionS2.C.3,Page18,\nseeSI).Insummary,thegeometriesoftheMDtrajectoryfromLoxodynamicssimulationcloselyresembleoptimized\nconfigurationalgeometriesonNEBpaths,suggestingthattheMDtrajectoryisalreadyneartheMEPofdiscovered\nreactionsandindicatingthatLoxodynamicscanguidesystemsapproximatelytowardscalculatedMEPdirectionsfor\ncertainreactions. FurthersimulationdetailsareavailableinSectionS2.C.3(Page18,seeSI).\n-Thediscussionsectionclaimsawiderangeofpotentialapplications,butthevalidationislimitedtoa\nfewrelativelysimplecases(exceptforthelastone). Amorerealisticassessmentofcurrentlimitations,\nalongwithsuggestionsforhowtotransitionfromexplorationtodetailedmechanisticmodeling,would\nstrengthenthemanuscript. Forinstance: afterdiscovery,cantheCVsproducedbySkewencoderbeused\ninconventionalbiasedsimulationstocomputefreeenergysurfaces?\n#Response#:\nWeappreciatethereviewerforthisinsightfulcomment. InthenewlyaddedSectionS2.B.3(Page15-16,seeSI),we\npresentresultsfromwell-temperedMetadynamicssimulationsemployingthelearnedSkewencoderCVinbothS 2and\nN\nDiels-Alderreactions. TheseSkewencoderCVsarederivedfromiterationsofLoxodynamicssimulationswherestate\ntransitionsoccurred. ByproperlyconfiguringMetadynamicsparameters,FigureS11(Page16,seeSI)demonstrates\nthattheSkewencoderCVcanfacilitatetransitionsbetweentwostatesinMetadynamicssimulationsofbothreactions.\nTheSkewencoderCValsoeffectivelydifferentiatesreactantandproductstates,correspondingtothestateinformation\nobtainablefromheuristicCVs. FigureS12(Page16,seeSI)illustratesthatsamplesfromwell-temperedMetadynamics\nusing Skewencoder CVs can reconstruct the FES for specific reactions. The reconstructed FES aligns with those\ngeneratedbyMetadynamicssimulationsemployingheuristicCVs. Consequently,SkewencoderCVscanbefurther\nappliedinconventionalbiasedsimulations,providedtheyencapsulateconfigurationalinformationwithinbothreactant\nandproductspaces.\nhasbeenexploredinanagnosticmannerusingLoxodynamics. Pleasecheckouranswertocomment4toReviewer#4.\nReviewer: MyearliercommentreferredspecificallytotheDiscussionsection,whichstillappearsun-\nchanged. Evenwiththeadditionalresultsnowincluded,thereremainsanoticeablegapbetweenthe\nscopeofthepresenteddemonstrationsandtheextensivelistofpotentialapplicationsdescribed:\n“... heterogeneouscatalysis,reactivesurfacerearrangements,electrocatalysis,enzymaticcatalysis,and\nasymmetriccatalysis. Furthermore,itmayofferinsightsintofundamentalscientificquestionssuchaspre-\nbioticchemistry—includingtheUrey–Millerreactions—aswellasenvironmentallyandtechnologically\nrelevantprocesses,includingatmosphericandcombustionchemistry,andpyrolysis. Likewise,itcan\nbeappliedeventonon-chemicalactivatedeventssuchasligand–proteinunbindingandpolymorphic\ntransitionsincrystals,wherenewphasesarenotknownapriori.”\nIwouldsuggestmakingthissectionmorebalancedandrealistic. Itwouldalsobehelpfultoexplicitly\nacknowledgethecurrentmethodologicalconstraints,forexample,theneedofsystem-specificrestraining\nwallsorstatedetections,andhowthesecouldbeimprovedinfuturework.\n#Response#:\nWethankthereviewerforthissuggestion. Accordingly,wehaverevisedthecorrespondingpartofthediscussionas\nfollows:\n\"Thisframeworkhasthepotentialtodrivesignificantadvancesinareaswheretheexplorationofactivatedeventsis\ncentral,particularlyincondensed-phasephenomena. Theseincludesurface-mediatedprocessessuchasheterogeneous\ncatalysis—thermal,asexemplifiedbythecasespresentedhere,electrochemical,andpotentiallyphotocatalytic—aswell\nasselectivity-drivenchallengesinasymmetric,enzymatic,andsupramolecularcatalysis. Furthermore,themethodmay\nprovideinsightsintoothercomplexsystemswithunknownpathways,suchasthosefoundinprebiotic,atmospheric,\nandcombustionchemistry. Conceptually, itisalsoapplicabletonon-chemicaltransformationslikeligand–protein\n7\nunbinding,proteinfolding,orcrystalnucleationandpolymorphism,offeringatooltoprobetransitionswherethefinal\nstateisnotknownapriori.\"\nWealsoexplicitlyacknowledgethenecessityofsystem-specificwallsandstatedetectionmethodsintheDiscussion\nsection,asfollows:\n\"Weacknowledgethatthecurrentimplementationreliesonsystem-specificsetups, includingmanuallydefinedre-\nstrainingwallsandstate-detectioncriteria,whichstillrequirelimitedbutnecessarychemicalintuitionandhuman\nintervention. Futureworkwillfocusonfurtherautomatingtheseaspects,forinstancethroughthedevelopmentof\nan autonomous agent capable of managing the exploration phase, launching multiple exploration runs, designing\nandadaptingbiasingpotentialsonthefly,andautomaticallydetectingnewstateswhilestoringandorganizingthe\nassociatedinformation.\"\n>>MethodologicalDetails\n-Itisunclearwhetherκischosenbytheuseroradaptedduringtheiterations(assuggested. IntheSN2\nexample,itissaid: \"Duetotheasymmetricenergybarriers... alargerharmonicconstantwasapplied...\ntoaccountforthehigherbarrier.\"Doessettingthisparameterrequirepriorbarrierknowledge?\n#Response#:\nWethankthereviewerforraisingthispoint. Inpractice,weinitiateasmallswarmofindependenttrajectories(typically\n5to10dependingonthecomplexityoftheproblem),eachcharacterizedbyadifferentvalueofthebiasingparameterκ.\nThisstrategyallowsustoexploremultipledirectionsefficientlyandidentifysuitablepathsforfurtherexploration.\nWeacknowledgethatthisdetailwasnotclearlyexplainedinouroriginalsubmission. Toaddressthis,wehaveupdated\nboththemainmanuscript(SectionIV.C,Page11,seemaintext)andtheSupportingInformation(FigureS1,Pages4,\nseeSI),wherewenowexplicitlydescribethisinitialswarmapproach. Additionally,werevisedtheflowchartinFigure\nS1toreflectthisprocedure,clearlyindicatingtheuseofaninitialsetoftesttrajectoriesthatguidethesubsequentpath\nexploration.\nIngeneral,noprecisepriorknowledgeofthebarrierisrequired.\nReviewer: While this answer is clear, I feel that this point is still not properly explained in the main\ntext. Inthegeneralworkflowsection,itisdescribedas: “Theprocedureisinitializedwithasequenceκ ,\ni\nwhereeachvalueisdistributedatuniformintervals...”Iwouldsuggestexplainingmoreindetaillelikein\ntheanswerprovidedhereinstead(“Weinitiateasmallswarmofindependenttrajectories... forfurther\nexploration.”) \n#Response#:\nWethankthereviewerforthisobservation. WehaveamendedtherelevanttextinSectionVI.Ctoaddressthispoint,as\nshownbelow:\n\"Theprocedureisinitializedwithasmallswarmofindependenttrajectories,eachcharacterizedbyadifferentvalueof\nthebiasingparameterκ . Thisstrategyallowsustoexploremultipledirectionsefficientlyandidentifysuitablepathsfor\ni\nfurtherexploration.\"\nForthedetailsoftheLoxodynamicsworkflow,wehaverevisedthecorrespondingtextintheSIasfollows:\n\"The procedure is initialized with a sequence {κ }, where each value is distributed at uniform intervals such that\ni\nκ =κ +δκ,withδκdenotingaconstantincrement.\"\ni+1 i\n-Therationalebehindtheinclusionofσinthewallposition(i.e.,µ+σ+ϵ)isnotexplainedinthemain\ntext. Figure1shouldalsoreflectthislogicexplicitly. Also,theroleofϵseemsimportantfordetermining\nthedegreeofexploration,butisnotdiscussed. Adiscussionofhowtochoosetheseparametersandwhat\neffecttheyhaveisneeded.\n#Response#:\nWeappreciatethereviewerforhighlightingthisissue.Inthenewlyrevisedmanuscript,weprovideadetailedexplanation\ninthemaintextregardingtheimpactofutilizingprojectedsamplevarianceσandacustomoffsetϵ,asillustratedin\nFigure7(Page10,seemaintext),ratherthanmerelyportrayingtheminFigure1asbefore.\nInthecurrentrevisedversion,Figure1(Page3,seemaintext)primarilyillustratesthattherelationshipbetweensample\nskewnessandlocalanharmonicityarisesfrombiasedsimulationsonanartificial1-DPES.Duringeachiterationof\nthe1-Dsimulation,abiasharmonicwallisintroducedandconsistentlyshiftedwithaconstantoffset∆xfromtheleft\ntotherightalongthex-axis. Weobservedthatasthesystemtransitionedfromoneminimumtoanother,thesample\ndistributionbecameincreasinglyasymmetric,reflectingitsgradualevolutiontowardsstatetransition. Inthecurrent\nmanuscriptversion,wehaveclarifiedthatthisobservationfromthe1-Dsimulationservedasinspirationforgoing\nfromasimpleidealcasetomorerealisticscenarios. Consequently,toaddresschallengesinextendingthisknowledgeto\nmultidimensionalscenarios,wedevelopedtheSkewencoderNNarchitectureandLoxodynamicsworkflowasinnovative\ntoolsforchemicalspaceexploration.\nThe impact of the custom offset ϵ has been addressed in the response to the Method Validation section of Major\nCommentsfromReviewer#2. Similarly,weemployprojectedsamplevariancetodeterminetheharmonicwall’sposition.\nouralgorithmincursnoadditionalcomputationalcost. Furthermore,samplevariancereflectsthedistributionwidth\nwithinlatentspace. Byusingthistermasanoffsetfortheharmonicwall,previouslyexploredregionsinlatentspace\narebypassedinnewsimulationiterations,reducingoverlapandenhancingexplorationefficiency. Inpractice,onemay\nreplaceprojectedsamplevarianceσorevenσcombinedwithϵwithaconstantoffset(asillustratedinFigure1,page\n3,seemaintext). ThissubstitutionyieldsvaryingeffectsacrossdifferentPES.Byintroducingσandamoderateϵ,\nwedevisedamechanismthatadaptivelyadjustsLoxodynamics’exploration\"speed\"basedonlocaldatasetswithout\nincurringextracomputationalexpense.\nReviewer: Iappreciatetheadditions,whichIbelievearehelpfultobetterunderstandthemethodology.\nHowever,Ithinkthatafewwordsonepsilonwouldstillbeuseful,evenjusttogiveanideaofwhich\nvalues should be used based on your experience if, as it seems, it is not a very critical parameter. In\naddition,Fig.7isnotexplainedproperlybythecaption,inwhichmanyelementsarenotexplained(colors,\narrow,symbols,abbreviations..). Inaddition,itisnotcleartomewhythenotationfortheshiftofthe\npotentialwallisdifferenthere(σandϵ)withrespecttothenotationusedinFig.1(∆x).\n#Response#:\nWethankthereviewerforthisfurthercomment. Wehaveexpandedthediscussionofϵinthemaintext,specificallyby\nprovidingarangeofvaluessuitableforexploration. WebelievethisfurtherelucidatestheroleofϵinLoxodynamics\nsimulations. Therevisionisasfollows:\n\"Valuesofϵ≤1.0areconsideredmoderate.\"\nWeapologizefortheconfusionregardingthesymbolchanges. Wehaverevisedthecaptionsofbothfigurestoclarifythe\nnotation. Specifically,wenotethat∆xinFigure1representsaconstantshiftwithinthebiasedsamplingschemethat\ninspiredLoxodynamics. Incontrast,σrepresentsastatisticalpropertyofsamplesprojectedontothelatentspace,while\nϵisauser-definedparameter. Ifoneweretodrawananalogy,thecombinedterm(σ+ϵ)inLoxodynamicsactsastheef-\nfectiveshiftinthelatentspace—analogousto∆x—withthekeydistinctionthat(σ+ϵ)isdynamicratherthanconstant.\n-Theclaimthatthereactiondiscoveryis“unbiased”ismisleading. Sincetheskewnessofthebiased\ndistributiondependsontheappliedparametersoftheharmonicbiasandnoreweightingisperformed,\nthemethodintroducesadirectionalsamplingartifact.\n#Response#:\nWe appreciate the reviewer for highlighting this issue. We would like to clarify that we do not claim the reaction\ndiscoveryprocesstobeunbiased,butratheragnostic. Althoughitistruethatnoreweightingisperformed,such\nanactionwouldprimarilyreproducetheequilibriumdistributionoftheinitialstate,providinglimitedinsightinto\ntransitionpathways. Therefore,acertaindegreeofdirectionalitymustbeintroduced. Importantly,thisdirectionalbias\nisnotimposedaprioribytheuserbutemergesduringtheearlystagesofthesimulationwhenthesystemexperiences\nminimalornobias. Asexplorationprogresses,themethodselectsspecificdirectionsforsamplingreactionpathways.\nThisselectionisnotmanuallyadjustedbutdrivenbyaniterativesamplingprotocolthatdynamicallyadaptstothe\nunderlyingfreeenergylandscape. Forclarity,wehaveremoved\"unbiased\"fromthesentence: Byfollowingpathways\nwiththelowestenergybarriers,thisapproachenablesefficientandunbiasedexplorationoflikelyreactionchannels\n(Page2,Line101,seemaintext).\n-Acomplexproceduretodecidewhethertowarm-starttheautoencoderornotisputforward,butnot\njustified. Acomparisonwithanexplorationtrainingonlyfromscratchonthetoymodelwouldclarifyits\nimpact. Also,thisprocedurerequiresacriteriontodeterminewhetheranewmetastablestatehasbeen\nreached,whichisnotdescribedfortheexamplesreported. Thiscontradictstheclaimofminimaluser\ninterventionandshouldbedetailedmoreclearly.\n#Response#:\nWeappreciatethereviewerforraisingthisinsightfulissue.\nIntherevisedversionofourmanuscript,wehaveincludedSectionS1.C(Page5,seeSI)toelucidatethecriterionfor\ndistinguishingbetweenmetastablestatesandnon-metastablestates. Generally,thiscriterioncanbeadaptedtoany\nconnectivitymatrix-basedmethodcapableofdeterminingatomicconnectivitythroughcustomizedswitchingfunctions.\nThepresentimplementationisanin-housesimplifiedversiondesignedforefficiency.\nAdditionally,acomparativeexperimentwasconductedbetweentheLoxodynamicssimulationutilizingthewarmstart\ntrainingstrategyandexplorationtraininginitiatedfromscratchonthetoymodelPES.Theresultsarepresentedin\nSectionS2.A.3(Page11,seeSI).Insummary,thewarmstarttrainingstrategyenhancesexplorationefficiencyby\nloweringtrainingandcomputationalcosts,therebyfulfillingtheintendedobjectiveofthisstrategywithinLoxodynamics\n10\nsimulation.\nReviewer: While I appreciate the additional material provided in the Supporting Information, I still\ndonotunderstandwhytheseelementsarenotatleastmentionedinthemaintext. Theyappearhighly\nrelevant,particularlyinlightoftheresultsshowninFigureS7. Abriefreferenceorsummaryinthemain\ntextwouldhelpreadersbettergrasptheirsignificance.\n#Response#:\nWethankthissuggestion. WehaveaddedabriefsummarytoSectionVI.Cofthemaintextasfollows:\n\"Todemonstratetheenhancedefficiencyofthewarmstartstrategy,wecomparedLoxodynamicssimulationsutilizing\nthisapproachagainstexplorationtraininginitiatedfromscratch(seeSI).\"\n-SectionAofthemethodscouldbenefitfromafewrevisions. Ithinkthattheworkingsofneuralnetwork\nlayersarenowpartofcommonknowledgeandsocanbemovedtotheSI.\n#Response#:\nWethankthereviewerforthisvaluablesuggestion. ThiscontenthasbeenrelocatedtoSectionS1.B.2(Page3,Line\n43-51,inSI),integratedwiththediscussionontheNNstructureofSkewencoder,inthenewlyrevisedversionofthe\nmanuscript.\nOn the other hand, it is not clear if the details about the L are just conceptual considerations (as\nskew\nit seems to me looking at the code) or technical details that need to be implemented to improve the\nperformanceofthemethod. Whicheverthecase,Ithinkthisshouldbemadecleartothereader.\n#Response#:\nWe appreciate the reviewer for identifying this issue. The current format of L encompasses both conceptual\nskew\nconsiderationsandtechnicalimplementation,enhancingthemethod’sperformance. RefertoPoint4fromReviewer#1\nandthediscussioninthenewlyaddedSectionS2.A.1(Page8-9,inSI)regardingthismatter. Insummary,thecurrent\nimplementationofL hasbeenshowntoimproveexplorationefficiencywhencomparedwithotherformsofL .\nskew skew\nAdditionally,SkewencoderutilizingL inapurelyGaussian-likeformsuccessfullyidentifiedalternativestatesand\nskew\nexhibitedonlyslightlyreducedperformanceduringLoxodynamicssimulationonthetoyPES,suggestingthatL\nskew\ncanbeappliedflexiblyandisthereforealsoaconceptualconsideration.\n-Ingeneral,itwouldbemuchbettertoincludeinthemethodssectionasummaryofthemethod,including\nthelossfunctions(whichareonlydefinedintheintroduction)andtheworkflowoftheiterativeprocedure\n(whichisagaindiscussedintheintroductionandsummarizedonlyintheSI).\n#Response#:\nWeappreciatethereviewerforthishelpfulsuggestion. Wehavereorganizedthemethodologysectioninthemaintext\naccordingly. Toeliminateredundancy,wecontinuetofocusondiscussingthedesignoftheskewnessloss-relatedterm\ninSectionIV.A(Page9-10,seemaintext),buthaveremovedtheintroductiontoageneralNNstructure. Subsequently,\nwepartiallyrevisedSectionVI.B(Page10,seemaintext). WeincludedFigure7(Page10,seemaintext),whichmore\nclearlyillustrateshowtheharmonicwallisappliedbasedonsamplemean,variance,andskewnessvalues. Additionally,\nweemployedasimplifiedequationformindicatinghowtheharmonicwallisappliedinlatentspace,assuggestedby\nReviewer#1inPoint3. Asaconclusiontothemethodologicalsection,wepresentasimplifiedversionoftheentire\nalgorithminSectionIV.C(Page11,seemaintext),correspondingtodetaileddiscussionsinSectionS1oftheSI.\nCommentsspecifictotheapplicationspresented.\n>>ToyModel\n-TheevolutionoftheskewnessandthemeanoftheCVacrossiterationsisnotshown. Includingthis\nwouldhelpillustratehowthemethodprogressesandadaptsthebiasingdirectionbasedonthemeasured\nskewness?\n11\n#Response#:\nWeappreciatethereviewerforhighlightingthisissue. Wewishtoclarifythatthelatentspaceundergoeschangesat\neachiteration,renderingcomparisonsofdistributionmomentsacrossdifferentvariablesimpractical.\nConcerningCVevolutionasrequestedbyReviewer#1,pleaserefertoPoint1fromReviewer#1whereweillustratethe\nevolutionofCVsintoymodels.\n- How does the exploration change as a function of the chosen hyperparameters, κ and ϵ? Showing\nsupplementaryexampleswouldhelpclarifyhowtosetthem.\n#Response#:\nWeappreciatethereviewerforthiscomment. Wehavepreviouslyaddressedtheimpactofϵinearlierquestions(e.g.,\nMethodValidationsectionofMajorCommentsfromReviewer#2). Regardingtheeffectofκ,asdiscussedinprevious\nfeedback,briefly,ahigherκtendstoidentifyreactionpathwayswithelevatedreactionbarriers,whereasasuitablylower\nκisinclinedtodiscoverpathwayswithreducedreactionbarriers,closelyapproximatingtheanticipatedMEP.Thus,\nselectingκcansignificantlyinfluencetheexplorationofdistinctreactionpathways. WehaveaddedSectionS2.A.2\n(Page9-11,seeSI)toassesstheeffectofκusingatoymodelinLoxodynamics.\nConcerningSkewencoderhyperparameters,weassertthatwithatypicaldeepautoencoderneuralnetworkstructure,\nSkewencoder can effectively learn features from latent space and drive Loxodynamics. We have already compared\nSkewencoderperformanceusingbothdeeperandshallowerneuralnetworkstructuresinSectionS2.B.2(Page13-15,see\nSI).Bothconfigurationsoftheneuralnetworksuccessfullydrivethesystemtoidentifystatetransitions.Therefore,given\ntherobustfittingcapabilityofneuralnetworks,ageneralAutoencoderstructurecanbeemployedwithinSkewencoder\nwithoutextensivehyperparameteroptimizationefforts. Thishighlightsourmethod’sflexibilityandscalability.\n- Although the potential includes two distinct paths, only one is sampled. Can the method recover\nboth,dependingonthehyperparametersused? Thiswouldalsohelpclarifyhowthisisobtainedinthe\nchabazitecase.\n#Response#:\nWeappreciatethereviewerforhighlightingthisissue. WekindlyrequestthereviewertorefertoPoint7fromReviewer\n#1andSectionS2.A.2(Page9-11,seeSI).\n>>SN2Reaction\n- The input CV space essentially consists of d and d , which already capture the transition without\nCl F\nneedinganycombinationofthem. Inmyopinion,thissystemshouldbemovedtotheSI,asitdoesnot\nsignificantlyvalidatethemethod.\n#Response#:\nWeappreciatethereviewerforthissuggestion. However,wewouldliketoemphasizethatalthoughtheS 2reaction\nN\nisarelativelysimpleexample,itservesasarealistictestcase. Giventhattheinputdimensionsofthepreviousmodel\nPESareall2D,andtheS 2reactioninvolvestwoSkewencoderinputdescriptors,itcanbeconsideredanappropriate\nN\nintermediarybetweenmorecomplexreal-worldexamplesandsimplermodelPES,whichareeasiertotest. Consequently,\nwehavedecidedtoretainthissectioninthemaintext,asitcontinuestovalidatethemethod.\n>>Chabazite\n-Howisthesecondpathwayfound? IntheSIitissaidthatdifferentelasticconstantsareusedtoreflect\nbarrierdifferences. Sincebarrierinformationisusuallyunknown,howwouldoneproceedingeneral?\n#Response#:\nThisiscertainlyafundamentalpoint. AsdiscussedinourresponsetoPoint6fromReviewer#1andillustratedin\nFigureS1(Page4,seeSI),ourmethodinitiatesasmallswarmofparallelexplorations,eachwithadistinctκvaluethat\nincreaseswithinaspecifiedrange. ThisapproachallowsLoxodynamicssimulationtoidentifycorrespondingreaction\npathwayswithdifferentκvalueswithoutrequiringbarrierinformation. IntheSI,wementionthatvariousconstants\nareemployedtoreflectbarrierdifferences,indicatingthatdifferentκvaluesconfiguringharmonicwallscorrespond\ntobarrierheightdifferencescalculatedusingDFTmethods(asreportedinhttps://pubs.acs.org/doi/10.\n12\n1021/jp513024z). Wedonotimplythatκdirectlycorrespondstoenergyvalues;rather,smallerκvaluesarelikely\ntouncoverreactionpathwayswithlowerbarriers,whereaslargerκvaluestendtorevealpathwayswithhigherbarriers.\nThisobservationalignswithourfindingsinthedoublepathwaymodelPESpresentedinSectionS2.A.2(Page9-11,see\nSI).\n-Inthestepwisemechanism,howisbackwardmotionprevented? Areadditionalwallsrequired? Would\nthisbecomelimitingformorecomplexcatalyticcyclesinvolvingmanyintermediatesteps?\n#Response#:\nWeappreciatethiscommentverymuchaswerealizethispointwasnotclearintheoriginalsubmission. Thus,we\nupdatedaccordinglySectionS2.C.2(Page18,seeSI),discussinghowwerestarttheLoxodynamicssimulationtoexplore\nareactionpathwaythatincludesmultiplereactionwells,e.g. thestepwisemechanismintheethanoldehydrationcase. In\nbrief,FigureS13(Page17,seeSI)showsthatoncewereachedtheintermediatestate,werestarttheLoxodynamcis,keep\nalltheparametersrelatedtotheSkewencoder,butweretaintheoldCValongwhichweapplyanadditionalrestraining\nwall that guide the system to identify the current intermediate state. Then in the newly restarted Loxodynamics\nsimulation,an“unbiased”simulationwillbeconducted,withtherestrainingwallkept. ThennewLoxodynamicsbias\nwillbeaddedinthefollowingiteration,withtherestrainingwallstillexisting. Tocompensatetheadditionalpotential\nbiascausedbytherestrainingwall,κvaluemustbeincreasedtoexplorethesecondpotentialwell,thoughthebarrier\nheight of the first reaction and the second reaction in the multiple-well reaction path are similar. We reported this\nobservationinSectionS2.C.2(Page18,seeSI)andS2.D(Page18-20,seeSI),correspondingtodifferentdehydration\nreactionsontheacidicChabasite. Basedonourobservation,weassumethatthismightbealimitationtoareaction\npathwaywithmorewells. However,inthetestcaseswehaveconductedLoxodynamicssimulationsonsofarwedonot\nobservesuchalimitation.\nReviewer: Alsointhiscase,theuseofawallfromthepreviousiterationisnotdescribedwellinthetext,\ndespitesoundingimportantinpreventingback-crossing. Thedetailedexplanationgivenhereshould\nthussomehowbeincorporatedinthemaintext. IalsobelievethatFig. S13shouldatleastbementioned\ninthemaintextcontextually,ifnotaddeddirectlytothemaintext.\n#Response#:\nWethankthereviewerforthissuggestion.Wehaveaddedabriefdiscussiontothemaintextregardinghowback-crossing\nispreventedafterrestartingaLoxodynamicssimulation. Therevisedpassageisasfollows:\n\"Topreventback-crossingtotheoriginalstateuponreachingthisintermediatestate, werestarttheLoxodynamics\nsimulation. WhiletheSkewencoderparametersandtheoriginalCVareretained, anadditionalrestrainingwallis\nappliedalongthisCVtoconfinethesystemtothecurrentintermediatestate. Thesimulationthenproceedswiththis\nconstraintactive. (Forfurtherdetailsregardingtherestartprotocol,refertotheSI.)\"\n>>MinorComments\n-Manyfigurecaptionsarenotself-contained. Forexample,themeaningofcolorschemesandsymbols\n(e.g.,Figure1symbols,Figure3colorbar,colorofscatterplots,Figure4colorofscatterplots,FigureS1\nsymbols,FigureS2colorbars,FigureS4colorofscatterplots,colorbars)isonlydefinedinthemaintext\nandshouldbereportedinthecaptions.\n#Response#:\nWeappreciatethereviewerforidentifyingthisissue. Below,wedetailthecontentaddedtothecaptionsoftherelevant\nfigurestoaddressthisconcern. Generally,weemphasizethattheprecisevaluesonthecolorbar(Figure3andFigure\nS2)forboththemodelPESandprojectedlatentspace(FigureS4)arenotcritical. Theclosedformathasalreadybeen\nprovidedforthePES.AsLoxodynamicssimulationsprogress,theprojectedCVspaceevolvescontinuously,withexact\nvaluesinprojectedlatentspaceCVbeingupdatediteratively. Consequently,relativerelationshipsaremoresignificant\nthanexactvaluesforclearlyillustratingCVevolution.\nInthecaptionofFigure1(Page3,seemaintext),weadded:\n-Theinitialwallislocatedatx=µ,correspondingtothemeanoftheunbiasedsampledistribution(indicatedbyred\ndotsintheupperpanel). Ineachsubsequentiteration,theharmonicwallistranslatedbyafixedoffset∆xrelativetoits\npositioninthepreviousiteration. Thespringconstantκforallwalls(depictedasblue,green,andyellowhalf-parabolic\ncurvesinthelowerpanel)remainsconstantthroughoutalliterations.\nReviewer: Asmentionedinapreviouscomment,theAuthorsshouldensurethatthenotationusedhereis\nconsistentwiththeoneusedinFig. 7(∆xvsϵandσ).\n#Response#:\nWethankthereviewerforraisingthispoint. Wehaverevisedthecaptionsofbothfigurestoclarifyourexplanation.\nInthecaptionofFigure3(Page5,seemaintext),weadded:\n-Thedarkerregionsignifiesanareawithlowerpotentialvalues,whereasthelighterregiondenotesanareawith\nhigherpotentialvalues.\nReviewer: Iwouldatleastindicatethespacingoftheisolinesinthecaption(sameforFigureS3below)to\nprovideaquantitativeinformation.\n#Response#:\nWethankthereviewerforthissuggestion. WehaveaddedacolorbartoFigure3,whichalsoalignswiththeeditor’s\nrequirements. ForFigureS3,wehaveexplicitlyindicatedthespacingoftheisolinesinthefigurecaption.\n- The central and lower panels display the evolution of the x (in dark green) and y (in light green) coordinates,\nrespectively,overtime(inLennard-Jonesunits). Thebluepointsrepresentsamplesfromtheincrementaldataset,while\nredpointscorrespondtosamplesfromthecurrentiteration.\nInthecaptionofFigure4(Page6,seemaintext),weadded:\n-Inthe2Dmaps,asinFigure3,bluepointsdenotesamplesfromtheglobaldataset,whereasredpointsindicate\nsamplesfromthelocaldatasetcorrespondingtothecurrentiteration. Relevantstructuresextractedfromthetrajectory\nblocksofeachsuccessiveiterationareshownintheuppersubpanels.\nInthecaptionofFigureS1(Page4,seeSI),weadded:\n- In this flowchart, {κ } represents a sequence of wall spring constants arranged in an arithmetic progression\ni\nwithacommondifferenceofδκ. Theparametern denotesthemaximumnumberofiterationspermittedforthe\niter\nLoxodynamicssimulation. Thequantityγ referstotheprojectedsampleskewnessinlatentspace. ThevariableL\nls state\nindicatesthelistofidentifiedmetastablestates,whileS servesastheidentifierforthecurrentstate.\np\nInthecaptionofFigureS3(Page7,seeSI),weadded:\n-Thedarkerregionsignifiesanareawithlowerpotentialvalues,whereasthelighterregiondenotesanareawith\nhigherpotentialvalues.\nInthecaptionofFigureS4(Page8,seeSI),weadded:\n-Ineachleftsubcolumn,Thebluepointsrepresentsamplesfromtheincrementaldataset,whileredpointscorrespond\ntosamplesfromthecurrentiteration.\n-LighterregionscorrespondtohighervaluesinthelatentCVspace,whereasdarkerregionsrepresentlowervalues.\n-Theprovidedcodeiswell-structuredanddetailed,whichensuresreproducibility. However,fromtheSI\nandthecode,itappearstorelyextensivelyonthemlcolvarlibrary,whichisnotcitedinthemaintext.\n#Response#:\nWe appreciate the reviewer for this suggestion. The relevant citation has been incorporated in the code and data\navailabilitysection(Page12,seemaintext).\n-Whenquantitiesarediscussedinthemaintext,refertothespecificequationnumbertohelpguidethe\nreader(e.g.,L inmethodssec. A).\nskew\n#Response#:\nWeappreciatethereviewerforthissuggestion. Wehaveaddedthefollowinginthemaintext(Page9-11,line631-633):\nThissectionisdevotedtothediscussionoftheformulationandimplementationoftheskewness-relatedlossterm,L\nskew\n(definedinEquation4),toaddressthisissue.\nRemarksoncodeavailability:\nTheprovidedcodeiswell-structuredanddetailed,whichensuresreproducibility. However,fromtheSIand\nthecode,itappearstorelyextensivelyonthemlcolvarlibrary,whichisnotcitedinthemaintext.\nReviewer #3\nRemarkstotheAuthor:\nIco-reviewedthismanuscriptwithoneofthereviewerswhoprovidedthelistedreports. Thisispartofthe\nNatureCommunicationsinitiativetofacilitatetraininginpeerreviewandtoprovideappropriaterecognition\nforEarlyCareerResearcherswhoco-reviewmanuscripts.\nReviewer #4\nRemarkstotheAuthor:\nIappreciatetheeffortthattheauthorshaveundertakentoreplytoallthereviewers’comments. Ifind\ntheadditionalexamplesandinformationincludedinthemanuscriptandSItodrasticallyimprovethe\nreader’sunderstandingofthemethodaswellasitsrelevance. Therefore,Iwouldbehappytorecommend\nthemanuscriptforpublicationinNatureCommunications. Ionlyhaveonefinalremarkfortheauthors:\ninsectionS3,acomparisonbetweenthefreeenergyprofilesobtainedfromachemically-inspiredCVand\ntheCVfromLoxodynamicsisreported(FigureS12).However,thefreeenergyprofileisCVdependentand\ncomparingtwoprofilesdirectlyisnotpersemeaningful(see,e.g.,https://doi.org/10.1016/j.jcat.2020.04.015\norhttps://doi.org/10.1063/5.0083423). IwouldfinditinterestingtoevaluatetheactualCV-independentfree\nenergydifferencebetweenthetwometastablestatesfromthesimulationsresults,toshowcasethatthe\nLoxodynamicsCVis(likely)fineanddiscrepanciesinthefreeenergyprofilewithchemically-inspired\nCVsaremeaningless.\n#Response#:\nWethankthereviewerforhighlightingthispoint. WecalculatedtheCV-independentfreeenergydifferencebetween\nthereactantsandproductsfortheSN2andDiels-Alderreactionsbyintegratingtheprobabilitydensitieswithinthe\nrespectivebasins. ThefreeenergydifferencesobtainedfromLoxodynamicssamplingmatchedcloselywiththoseobtained\nusingchemically-inspiredCVsinthesetwotestcases. FurtherdetailshavebeenaddedtotheSI,asrevisedbelow:\n\"Thefreeenergydifferencesbetweenreactantsandproductswerecalculatedbyintegratingtheprobabilitydensities\nwithintherespectivebasins. Thedividingboundarybetweenbasinswassetto−0.28forthechemically-inspiredCVand\n0.0fortheSkewencoderCV.Thechemically-inspiredCVsimulationyieldedafreeenergydifferenceof143.177kJ/mol,\nwhile the Skewencoder CV simulation resulted in 143.204 kJ/mol, demonstrating excellent agreement between the\ntwo simulation results. ... Still, we calculated the free energy difference by integrating over basins defined by a\nthresholdof−2.2(chemically-inspiredCV)and−0.7(SkewencoderCV).Theresultingvalueswere106.223kJ/moland\n102.342kJ/mol,respectively,indicatingreasonableagreementbetweenthetwulationresults.\""}]

Decision letter