Thermodynamic and Kinetic Analysis of Molecular Conformational Dynamics in a Riemannian Framework

TL;DR

This paper introduces a new geometric framework using Riemannian geometry to improve the accuracy of molecular simulations by addressing issues with coordinate transformations.

Contribution

A novel Riemannian framework is introduced to ensure invariance of key thermodynamic and kinetic quantities under coordinate transformations.

Findings

The framework resolves noninvariance issues in conventional PMF definitions through Riemannian geometry.

A generalized Riemannian diffusion model enables rigorous estimation of kinetic properties in CV spaces.

Bayesian integration with the framework improves the reliability of free energy and transition rate calculations.

Abstract

We have formulated
a Riemannian framework for describing the geometry
of collective variable (CV) spaces in molecular simulations and demonstrate
its applicability through both toy model potentials and a biomolecular
example. This formalism addresses significant theoretical challenges
arising from the inherent nonlinearity of CV transformations, ensuring
critical quantities such as the potential of mean force (PMF), minimum
free energy path (MFEP), and rate constant remain invariant under
coordinate transformations. Our framework identifies and addresses
the noninvariance issues of conventional PMF definitions, clearly
illustrating their limitations through simple illustrative examples.
To overcome these issues, we introduce invariant definitions of PMF
and MFEP using Riemannian geometry. Moreover, we propose a generalized
Riemannian diffusion model applicable to diffusive dynamics within
CV spaces, allowing rigorous estimation of kinetic properties. Through
this model, we derive practical numerical methods for determining
the PMF, diffusion constant, metric tensor, and transition rates from
unbiased simulations conducted along identified transition pathways.
By integrating Bayesian approaches with the Riemannian framework,
our method provides a statistically robust technique for accurately
calculating free energy landscapes and transition kinetics, thereby
enhancing the reliability and interpretability of biomolecular simulations.

Figures

Toy Model: (A) A one-dimensional toy model described by potential energy V(x) = x 2. (B) Several CVs (ξ­(x)) are used as examples of smooth coordinate transformations ( x→ξζ ). (C) The conventional PMF along the CVs defined in (B). (D) The projection of conventional PMF (shown in C) onto the x space. (E) The Riemannian PMF along the CVs defined in (B). (F) The projection of Riemannian PMF (shown in E) onto the x space (all PMFs are exactly the same when projected onto the x space, irrespective of the CV used).

Comparison of thermodynamic and kinetic quantities for a 2-D toy model, analyzed using conventional and Riemannian frameworks. (A) Muller-Brown potential energy surface. The gray line represents the true MFPT connecting the two primary energy minima. (B) Diffusion profiles along the arbitrarily parametrized path r in the log-polar (X, Y) spcae. The black line shows the constant diffusion coefficient D used in the Riemannian framework. The blue line shows the position-dependent diffusion coefficient D(r) = D/h(r). (C) The ground truth free energy profile, G(s), calculated along the true arc length s of the MFEP in the Cartesian (x, y) space. (D) Free energy profiles calculated along the same arbitrarily parametrized path r.

Biomolecular example (alanine dipeptide): Comparison of thermodynamic and kinetic quantities for the C7 eq → C7 ax transition, analyzed using conventional and Riemannian frameworks. (A) Two-dimensional free-energy surface in the (ϕ, ψ) space; the gray curve shows a transition pathway estimated as MFEP, discretized into 20 equally spaced bins. (B) PMF profiles along the path coordinate (r = 0 and r = 1 correspond to C7 eq and C7 ax , respectively), computed using the conventional definition and the Riemannian definition. (C) Diffusion profiles inferred along the same discretized pathway: the conventional position-dependent diffusion D(r) and the Riemannian effective diffusion D/h(r). (D) The corresponding one-dimensional metric h(r) estimated along the pathway.

Full text

Introduction

1

Biomolecular simulations have made enormous progress in recent years. The advent of molecular dynamics (MD) in the study of biomolecular processes ?−? ? has given researchers new molecular level insights into previously unobservable phenomena. MD has overcome ?,? the common limitation of experimental techniques which forces the researchers to choose between high-resolution static (such as X-ray crystallography) or low-resolution dynamic (such as single-molecule fluorescence resonance energy transfer spectroscopy) pictures of biomolecular systems. MD does however have a few drawbacks. A key issue is the “time scale gap” in that MD simulations typically have a shorter time scale than many relevant biological phenomena. A related issue is the metastability; the system gets trapped in local free energy minima, preventing the system from evenly visiting the entire free energy landscape. Various enhanced sampling techniques and path-finding algorithms have been developed over the last few decades in order to overcome these limitations. ?−? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? These techniques have proven quite successful for simple toy models (e.g., dialanine peptide ?,?,?−? ? ? ? ) and, over the last few decades, have been successfully applied to a wide range of realistic biomolecular processes, including protein folding and binding, ?,? conformational transitions, ?−? ? and ligand binding. ?,? Despite this significant progress, their application to increasingly complex systems remains a considerable challenge.

Perhaps the most obvious line of attack for solving the sampling problems is simply to have stronger computing power and more efficient or specialized computer hardware. In the past few decades, extraordinary advances have been made in this regard. Current state of the art architectures such as very large ?−? ? ? or MD specialized ?−? ? supercomputers and GPU-enabled computing ?−? ? ? ? have capabilities that allow for simulations with large number of particles to run over long periods of time. The size of supercomputers have also given rise to algorithms developed for speeding up the calculations within brute force MD such as particle mesh Ewald? and dynamic load balancing.? The focus of this article, however, is the enhanced sampling techniques that are specifically based on enhancing the sampling within a statistical mechanical framework, in particular those methods that, in some form or another, use CV within their formalism.

In order to obtain relevant thermodynamic and kinetic? properties of a statistical mechanical system, one must integrate over high dimensional spaces, which requires large sets of independent and identically distributed samples. In order to manage the dimensionality of biomolecular systems, the assumption is often made that the vast majority of the high-dimensional space is practically “empty” (or occupied by microstates of negligible probability) and the occupied space can be approximated by a lower-dimensional manifold which contains all relevant conformations from stable states to important transition states. Ideally, any elementary reaction, can be thought of as a transition between two stable states, characterized by a transition pathway ?,? (or committor function?).

The dimensionality reduction may be employed explicitly or implicitly in a sampling scheme. For instance, a path-optimization technique ?,?,?−? ? ? ? ? ? ? ? can be thought of as a dimensionality reduction technique, if the optimized transition path can be approximated as a thin transition tube, where the areas outside the tube are associated with low probabilities. In this case, one may parametrize this path and potentially deviation from this path ?,?,? and use them as CVs. Other methods attempt to identify the intrinsic manifold by using statistical learning methods such as principal component analysis,? time-lagged independent component analysis,? isomap,? diffusion map, ?,? and autoencoders. ?−? ? Often these techniques are used to analyze MD trajectories, ?,?,?,? or to build Markov state models ?,? or are combined with enhanced sampling as in metadynamics ?,? or adaptive biasing force,? or integrated within adaptive sampling frameworks.? In this manuscript, we focus on the geometric structure induced by a chosen set of CVs, rather than on how those CVs are discovered.

Whether a set of CVs is defined in a systematic manner as described above or it is defined intuitively, it is a convenient way of reducing the dimensionality of the configuration space in both free energy calculation methods and path-finding algorithms. Various algorithms have been developed to estimate free energies ?,?,?,?,?,?−? ? ? ? ? or find transition pathways ?,?,?,?−? ?,? in CV spaces. What is mostly missing is a robust theoretical framework that allows for a rigorous treatment of the issues one needs to deal with when working with CVs. For instance, the CVs used in collective-variable-based enhanced sampling methods are often nonlinear transformations of atomic coordinates. This complicates their application as noted previously by Johnson and Hummer,? who show the conventional MFEP obtained from various path-finding algorithms are CV specific and are not invariant under nonlinear coordinate transformations. Such difficulties have often been ignored in the past in the majority of the applications of the collective-variable-based simulations; however, there has been attempts in addressing them as in the aforementioned work? or within the framework of Transition Path Theory. ?,? We also previously introduced a Riemannian framework for the rigorous treatment of the collective-variable-based spaces within the context of enhanced sampling and path-finding algorithms,? where the distance, free energy, and MFEP are redefined in a way that are all invariant quantities and do not change under coordinate transformations. Our previous work focused on the thermodynamic characterization of protein dynamics within a Riemannian diffusion model.? We did not previously discuss how to derive rate constants that, like the Riemannian PMF, are invariant under coordinate transformation. Here we extend the formalism in two different directions. First, in Section, we show that the Riemannian formulation is more general than its diffusion model manifestation and applies to a broad class of CV-based enhanced sampling and analysis methods. In Section, we extend the Riemannian diffusion model by focusing on the kinetics and rate calculation techniques, which we previously did not elaborate on.? In Section, we describe how to estimate transition rates in practice, and in Section we use a 2-D toy model to illustrate how the correct geometric treatment of CV spaces affects both thermodynamic and kinetic properties. While the 1-D and 2-D examples clarify the underlying concepts, we also apply the framework to alanine dipeptide in Section to demonstrate its feasibility in a realistic biomolecular system.

From Euclidean to Riemannian

CV Space

2

Euclidean geometry was developed by the Greek Euclid ca. 300 BC. Early in the 19th century, mathematicians such as Gauss,? Schweikart,? Bolyai? and Lobachevsky ?,? began to formulate non-Euclidean geometries. With the work of Bernhard Riemann, specifically his lecture “On the Hypotheses which lie at the Bases of Geometry”,? published first in 1873, geometry started exploring new and more diverse applications. Riemann’s work on generalizing the differential surfaces of $[eqn]$ ^3^ led to progress in many fields of science. Further progress upon his ideas allowed for the formulation of Einstein’s General Theory of Relativity? and progress in group theory.?

Riemannian geometry provides a robust mathematical framework to develop a formalism for the geometry of CV spaces, that are often defined to reduce the dimensionality of the atomic models of macromolecular systems. For instance, consider a transmembrane protein whose transmembrane helices rotate under certain conditions to allow opening or closing of a gate and transporting materials across the membrane. An intuitive CV for such a system would be the orientation of the transmembrane helices that can be determined using principal axes of the helices or their orientation quaternions. ?,? In order to work in these spaces, we must first have a formalism that allows us to answer common questions which are posed in a typical collective-variable-based simulation. Examples of such questions include: what is the distance between two points in the CV space? How is the potential of the mean force (PMF) defined at a given point in the CV space? How does the system diffuses along a transition pathway? How can we estimate the rate of a transition along a transition pathway? The questions have previously been answered within a Euclidean framework but with our Riemannian treatment of the CV space, some of these concepts and quantities need to be revisited. We will begin to answer these questions one at a time.

Imagine a system containing N atomic coordinates described by position vector ** x ** under a potential energy surface V(** x **). In order to reduce the dimensionality of the system and quantify important functional states, a coarser space is desirable to be defined such that ξ: $[eqn]$ ^n^ → $[eqn]$ ^m^, where n = 3N (N being the number of particles) and ξ is a multidimensional CV. The PMF is typically defined as

$$[eqn]$$

where Ã(ζ) is the PMF at any given point ζ in the CV space and δ(.) is the conventional Dirac delta function.

The PMF as defined above is sometimes considered as the effective potential energy of the reduced system (i.e, a system which is described by the CV and not the atomic coordinates). For multidimensional CV spaces, the MFEP or other related 1D pathways are often used for extensive sampling and characterization in place of the entire multidimensional space. In other words, the multidimensional CV space can be reduced itself to a 1D curve defined in the CV space.

The above approach is quite common and provides a powerful tool for characterizing the energetics of large biomolecular systems. However, we argue here that the above definitions of PMF and MFEP are not well-suited for the purposes that they are defined for. The main problem with these quantities is that they are not invariant under coordinate transformation. Unlike our previous work,? the discussions in this manuscript do not make any assumption regarding the diffusivity of the effective dynamics in the CV space and are more general. In the next Section, we will discuss the implications of the diffusivity and the Riemannian diffusion model.

Let us consider a simple toy model to clearly illustrate the noninvariance of conventional PMF. Our toy model is a 1D system in thermal equilibrium with a heat reservoir of temperature T, where β = (k _ B _ T)^−1^ = 1. The system is governed by the potential energy V(x). The PMF along any CV ξ­(x) is defined as

$$[eqn]$$

Defining η­(ξ) as the inverse function of ξ­(x) such that x = η­(ξ), and assuming dη/dξ > 0 everywhere, we can show:

$$[eqn]$$

In other words, if $[eqn]$ , we have

$$[eqn]$$

Projecting the CV back to the x space, we have

$$[eqn]$$

where we have used ξ′(x) = 1/η′(ζ) since η is the inverse function of ξ.

Relations (?) and (?) clearly show that the shape of PMF is not determined by the potential energy of the system only but it also depends on the derivative of the CV. One may use any free energy calculation method (such as umbrella sampling,? metadynamics, ?,? or ABF?) to calculate the PMF; however, without taking into account the second term in Relation or ?, the PMF’s shape or the shape of its projection onto the x space, does not represent the underlying potential energy for such a 1D system.

As an example, let us assume V(x) = x ^2^. The shape of PMF may or may not be similar to this potential energy, depending on the definition of CV. Figure shows V(x) along with the PMF as a function of several CVs (up to an additive constant) that were specifically designed to result in various PMFs. If ξ­(x) = x, the PMF would be the same as V(x). However, if ξ­(x) = erf(x), the PMF becomes flat since:

$$[eqn]$$

Using Relation, one may solve the following nonlinear differential equation to generate any arbitrary PMF function Ã(ζ):

$$[eqn]$$

or by some rearrangement:

$$[eqn]$$

Once the above differential equation is solved for η­(ζ), one can easily find its inverse function ξ­(x). We used this procedure using numerical methods to find functions α­(x) and β­(x) that result in Ã(ζ) = −ζ^2^ and Ã(ζ) = sin­(10ζ), respectively. Figure illustrates how the PMF, as defined conventionally based on these CVs, could qualitatively behave differently from what we expect intuitively from the underlying energetics (here, the potential energy) of this simple system. This is true when PMF is plotted against ζ (which is generally expected as ζ and x are not the same) and more importantly when it is projected back onto the x space.

Toy Model: (A) A one-dimensional toy model described by potential energy V(x) = x 2. (B) Several CVs (ξ­(x)) are used as examples of smooth coordinate transformations ( x→ξζ ). (C) The conventional PMF along the CVs defined in (B). (D) The projection of conventional PMF (shown in C) onto the x space. (E) The Riemannian PMF along the CVs defined in (B). (F) The projection of Riemannian PMF (shown in E) onto the x space (all PMFs are exactly the same when projected onto the x space, irrespective of the CV used).

The above example clearly shows that the PMF could be quite misleading if interpreted as a typical potential energy surface, where the minima are considered as locally stable states and the maxima are interpreted as transition states. However, since the CV function ξ­(x) is known, log­(ξ′(x)) can also be calculated and subtracted from the PMF to result in V(x). To do this, first one needs to project Ã(ζ) onto the x space by using ζ = ξ­(x) before subtracting log­(ξ′(x)) (Relation). Alternatively, one may add log­(η′(ζ)) to Ã(ζ) to get V(η­(ζ)) (Relation) and then project that onto x to get V(x). We define an invariant PMF as

$$[eqn]$$

Note that A(ζ) is a function of ζ similar to conventional PMF but A(ξ­(x)) = V(x) for any ξ­(x). We argue A(ζ) is conceptually more useful than Ã(ζ) since its connection to V(x) is more straightforward and its minima and maxima correctly represent the minima and maxima of V(x). Figure illustrates how the same CVs used for conventional PMF calculations can also be used to calculate the invariant PMF A(ζ). The resulting PMFs all qualitatively look similar but once projected back onto the x space, they all result in an identical function, i.e., V(x).

The above definition of invariant PMF can be reformulated as below

$$[eqn]$$

Assuming x = η­(ζ) is the only solution to ζ = ξ­(x) (i.e., assuming ξ is an invertible function of x), we can show δ­(ζ – ξ­(x)) = δ­(x – η­(ζ)) /|ξ′(x)| (which is the result of identity δ­(y) = ∑_ i _ δ­(x – x _ i _) /|y′(x _ i _)|, where x _ i _’s are the roots of y(x)). Therefore, using Relation, we have

$$[eqn]$$

After illustrating the noninvariance nature of conventional PMF and showing that the problem can be resolved by modifying the definition of PMF, we can now generalize the definition of PMF to multidimensional spaces. Relation is very easily generalizable for a full transformation from n-dimensional ** x ** to m-dimensional Ξ(** x **), where we also no longer assume β = 1:

$$[eqn]$$

in which J Ξ is the determinant of Jacobian ** J ** Ξ of transformation Ξ, i.e., $[eqn]$ where i = 1,···, m and j = 1,···, n, ** Z ** is any given point in the CV space, η(** Z **) is the inverse map that returns the unique ** x ** coordinates whose CV value is ** Z . Unfortunately, the full transformation is generally not desirable and η( Z **) is not available.

Of course, J Ξ, is defined only when n = m. When a CV ξ(** x **) is lower-dimensional than ** x , one can assume ξ( x ) is part of a full transformation. First, we assume ξ­( x ) is one-dimensional to simplify the discussion. We can keep the definition of ξ­( x ) general, yet choose Ξ( x ) to be the full transformation, such that it contains ξ­( x ) and a (N – 1) -dimensional vector orthogonal to ξ­( x ), denoted by ϕ( x ). Together, Ξ( x **) = (ξ, ϕ), and orthogonality implies ** ∇ ξ· ∇ **ϕ_ i _ = 0 for all i. Here ** ∇ ** is the gradient in x. Now we can write

$$[eqn]$$

Since the determinant of the product of matrices is equal to the product of their determinants, we can also write

$$[eqn]$$

where $[eqn]$ is the determinant of matrix ** J ** Ξ ^2^ (i.e., $[eqn]$ ). Subsequently, one can write

$$[eqn]$$

Since ξ­(** x ) is orthogonal to all ϕ_ i _( x **), one can easily show that J Ξ = J ξ J ϕ, where J ξ ^2^ = ** ∇ ξ. ∇ **ξ and J ϕ ^2^ is the determinant of a (N – 1) × (N –

1. matrix containing elements ** ∇ ϕ_ i _. ∇ **ϕ_ j _. We have

$$[eqn]$$

where ** Z ** = (ζ, θ). The N-dimensional invariant PMF can be used to define the 1D invariant PMF A(ζ) as

$$[eqn]$$

Note that J ϕ ^–1^ appears in the integral since we have change the variable to ϕ. Inserting into the above relation gives

$$[eqn]$$

A(ζ) is invariant since (1) A Ξ(** Z **) is invariant and (2) A(ζ) does not depend on the choice of ϕ as long as ϕ is orthogonal to ξ.

Now we define what we refer to as the metric g, using the conditional ensemble average of J ξ:

$$[eqn]$$

With the above definition of metric, one can now write

$$[eqn]$$

By introducing a fixed reference constant with the same units as g, g 0 ≡ g(ζ_ ref _), we can safely take a logarithm:

$$[eqn]$$

The above definition of 1D invariant PMF A(ζ) can be generalized to any arbitrary number d of dimensions:

$$[eqn]$$

The metric inverse g ^–1^ is defined to satisfy the relationship g ^–1/2^ = ⟨J ξ⟩_ ζ=ξ(** x ). This can be achieved by constructing elements of metric inverse g ^–1^ as g _ ij _ ^–1^ = ⟨(J ξ ^2^) ij ⟩ ζ=ξ( x ), where matrix J ξ ^2^ is composed of components ** ∇ **ξ i _· ∇ **ξ_ j _.

The above definition of invariant PMF differs from the conventional definition of PMF as it involves the metric tensor g. The metric tensor is a well-known quantity in differential geometry and is needed to do invariant measurements in non-Euclidean spaces. The PMF is only an example of a quantity that can be redefined to be invariant with the help of the metric tensor. The noninvariance of the conventional PMF and any other geometric quantity in the CV space is known, although it is often overlooked. The PMF can be made invariant by adding an additional term, which depends on the derivatives of CVs, as derived above. The difference between the invariant PMF and conventional PMF is related to the Fixman potential as also discussed elsewhere. ?,? Note that the Fixman potential was originally developed to relate the PMF associated with a constrained dynamics to that of an unconstrained one. ?,? Closely related “geometric/covariant” PMFs have been advocated in recent work. ?−? ? ? These insert the same metric/Jacobian factor, either to define a Transition State Theory-consistent geometric PMF,? to obtain a CV-form-independent activation free energy via the mass-metric,? or to derive a unique covariant PMF from stochastic dynamics with g = D ^–1^,? where D is the diffusion tensor.

A conceptually more straightforward approach to ensure the invariance of not only the PMF but also any other quantity of interest is to treat the CV space as a Riemannian space. In that view the metric tensor replaces the Euclidean dot product, so all distances, gradients and integrals are performed with g. For instance, the ordinary volume element d ^ N ^ξ is replaced by the invariant one $[eqn]$ ; likewise, the invariant gradient is $[eqn]$ , where ** e ^**_ i ’s are unit vectors along ξ i ’s. g _ ij _ = ⟨∂ k ξ i ∂ k ξ j _⟩ measures how much the underlying Cartesian space stretches when mapped into CV space. Riemannian geometry provides the conceptual framework and mathematical tools to treat the nonlinear behavior of curved spaces that are smooth but have potentially different curvatures at different points of space. The rigorous treatment of the CV space fits well within the Riemannian geometry framework. In this framework, the Riemannian PMF is defined exactly the same way as the conventional PMF is defined but the Riemannian Dirac delta function replaces the conventional Dirac delta function:

$$[eqn]$$

in which δ_ ζ (ξ) is the Riemannian Dirac delta function. Comparing Relation to Relation, one can make a connection between the Riemannian and conventional Dirac delta functions: δ­(ζ – ξ(** x ))|J ξ( x **)| = δ ζ _(ξ(** x **)).

We note that A(ζ) is the same invariant object (up to notation and a constant) as the “geometric/gauge-invariant” PMF. ?,?,? Our approach is not a different correction but a different framework: we treat CV space as a Riemannian manifold, so the profile and the calculus on it (volume element, gradients) are invariant by construction, rather than applying an a-posteriori fix to a Euclidean profile. We also note that the difference between the conventional and Riemannian PMF, $[eqn]$ , does not affect free-energy differences obtained from basin integrals when the same CVs cleanly separate the metastable states as discussed elsewhere.? By contrast, pointwise quantities read off a conventional PMF well depths, barrier heights, curvatures-depend on the parametrization of the CV (see Example 1, FigureC–D). The practical magnitude of this dependency is determined by the variation of the metric term; for simple CVs such as Euclidean distances, the metric is often nearly constant, resulting in a uniform shift that preserves relative barrier heights. However, for nonlinear reparameterizations or curvilinear CVs, the position-dependent metric can introduce significant corrections. Using the Riemannian PMF (Rels. ?–?) removes this CV dependency and collapses the profiles across CV choices (See Example 1, FigureE–F).

With the same argument made above, one can show the conventional MFEP is not invariant under coordinate transformation as shown by Johnson and Hummer.? Our Riemannian treatment of the CV space, however, provides a robust framework for defining the MFEP in an invariant way.? The MFEP in a Euclidean space is defined as any path whose tangent is parallel to the free energy gradient, i.e., dξ/ds∥∇_ ξ _ Ã(ξ) with ∇_ ξ _ A = ∑_ i _ (∂Ã/∂ ξ)** e ~**_ i _. As à itself is not invariant, that MFEP is not invariant under coordinate transformation, which questions its importance as a meaningful quantity. The Riemannian framework allows us to define the MFEP simply based on the Riemannian/invariant gradient of Riemannian/invariant PMF. In other words, not only both PMF and its gradient are well-defined, invariant quantities within the Riemannian framework, the Riemannian MFEP that is simply a path parallel to the gradient of PMF, $[eqn]$ , is also well-defined and invariant under coordinate transformation. With the invariant definitions of PMF and MFEP established, we now require an appropriate dynamic model consistent with these definitions. The Riemannian diffusion model, detailed next, addresses precisely this need, providing a robust theoretical basis for kinetic characterization.

Riemannian Diffusion and

Transition Rate Estimation

3

In the previous Section, we did not make any assumptions regarding the dynamics. However, quantities such as PMF and particularly MFEP are difficult to interpret if the reduced system follows a nondiffusive dynamics in the CV space. The intuitive interpretation of minima and saddle points of PMF representing the stable and transition states and the MFEP representing the most probable pathway relies on the diffusive nature of the effective dynamics. Such a condition is only satisfied if at all with specific choices of CVs. Here, the focus of our discussion is not on how to find such CVs. However, if we can successfully identify a set of CVs such that the projected motion of the system on this reduced CV space is diffusive, we can write

$$[eqn]$$

where D is the diffusion constant and b _ i _ = −Dg ^ jk ^Γ_ jk _ ^ i ^, in which Γ_ jk _ ^ i ^’s are Christoffel symbols and Einstein summation convention is used. Note that D is not assumed to be position-dependent here. Instead the position dependence is absorbed in metric tensor ** g . d W ** is a Riemannian Wiener process, where ⟨W ^ i ^(t)⟩ = 0 and ⟨ Ẇ ^ i ^(t 1)Ẇ ^ j ^(t 2)⟩ = g ^ ij ^δ­(t 1 – t 2), where Ẇ ^ i ^ is the time derivative of W ^ i ^.

Rewriting this equation as Fokker–Planck or Smoluchowski equations,? we have

$$[eqn]$$

which contains the Laplace-Beltrami operator, $[eqn]$ , and u(ζ, t) is the probability of finding the system at ζ at time t with a boundary condition of u(ζ, 0) = u 0(ζ). This summarizes our Riemannian diffusion model, which was previously introduced in ref? In the following we derive a new relation that allows the estimation of the PMF, metric, diffusion constant, and transition rate from unbiased simulations performed along an approximate MFEP.

We have previously discussed how one can find a Riemannian MFEP using the Riemannian implementation of string method with swarms of trajectories (SMwST).? Let us assume we have found such a pathway, ξ(s), parametrized by its arclength s. At any given point along this path, ê _ s _ is the unit vector parallel to ξ(s). An N – 1-dimensional submanifold (Σ_ s _) can be defined perpendicular to ê _ s . Let us assume local coordinates ζ = (s, κ) describes any point on Σ s _. One can write

$$[eqn]$$

where ∇Σ s _ _ denotes the ∇operator in the submanifold Σ_ s _.

Let us now define U(s, t)≡ ∫Σ s _ _ u(κ, t)** dΩ ** _ κ , where Σ s _ is an arbitrary portion of the subspace of κ around its origin. On the MFEP, let us also define a univariate PMF G(s) = A(ξ­(s)). Continuing, we can integrate over individual terms in Relation, where the LHS of the equation becomes $[eqn]$ and the RHS terms can be approximated as follows. Assuming u(ζ, t) is much larger within a relatively narrow “tube” around the MFEP, we can integrate over the cross section of this tube and the areas around the tube. In other words, we choose Σ_ s _ to be a portion of the κ space that covers the cross section of transition tube that falls within this space as well as some low-probability areas around it. The first term of the RHS of Relation will thus be

$$[eqn]$$

where we assume inside a cross section Σ_ s _ the gradient ∂_ s _ A varies negligibly, so we treat it as a constant equal to its value on the MFEP center line A(ξ­(s)) = G(s). The next term becomes $[eqn]$ after integration. For the third term, we assume the following integral is negligible as ∇_ Σs _ A is nearly zero on the tube wall:

$$[eqn]$$

where the divergence theorem is used and ** n ** is the outward normal in the cross section. Inside the tube and far from its boundary ∇Σ s _ _ A is replacable with ∇Σ s _ _ G(s) which is zero by definition. Similarly, we can use the divergence theorem to reduce the last term to

$$[eqn]$$

We assume this term is also negligible since u is already vanishingly small at the tube boundary and the normal component of its transverse gradient is negligible. These approximations reduce Relation to a one-dimensional diffusion equation in terms of probability density U(s, t) and the univariate PMF (or potential energy) G(s):

$$[eqn]$$

The thin transition tube, used to simplify Relation is a common assumption made for path-finding algorithms such as string method.? At this point, we have reduced our atomic model, first to a coarse variable space, and then we have focused upon a single dimension, namely the transition pathway of interest.

Relation is quite similar to the conventional Smoluchowski equation in a 1D Euclidean space. However, this is due to the fact that we assumed s to represent the geodesic distance along the MFEP path. A slightly more general relation can be derived based on (?) for an arbitrary parameter r that merely labels points on the MFEP path. The r-space is then endowed with a 1D metric h(r). We have ds ^2^ = h(r)dr ^2^ or $[eqn]$ . The Riemannian Smoluchowsky equation equation can therefore be written in terms of r as

$$[eqn]$$

which is slightly different from a conventional 1D position-dependent diffusion equation. Here D/h(r) is equivalent to the conventional 1D position-dependent diffusion constant and G(r) is the same as the conventional 1D PMF in terms of r with an extra term $[eqn]$ , which makes the Riemannian PMF independent of the choice of r.

Finally, we can be recast the Relation as

$$[eqn]$$

Now suppose that one has been able to identify the MFEP without necessarily quantifying the metric. Relation provides a framework to determine the dynamics of the system as long as it stays close to the MFEP transion tube. The PMF and metric along r (G(r) and h(r)) fully describe the diffusive dynamics and the arbitrary diffusion constant D simply determines the unit of h(r). The following section shows how Rel. ? can be discretized and statistically estimated, yielding a practical rate matrix ** R ** that can be used to determine thermodynamic and kinetic quantities such as free energy the G(r), the metric h(r), and the mean first passage time (MFPT). However, none of these analyses would be meaningful without first identifying the MFEP.

Determining the Rate of Transition

4

Starting from the Riemannian 1D Smoluchowski Rel. ?, we adopt Hummer’s Bayesian discretization scheme,? generalized to include the metric h(r). The Riemannian Smoluchowski Rel. of ? can be discretized such that r takes discrete values r _ n _ in the range [r _ A _, r _ B _] on a grid size Nδr (= r _ B _ – r _ A ) with mesh size δr = r _ n+1 – r _ n _. Applying a second-order finite difference scheme to the ? gives

$$[eqn]$$

where the rate element R _ n, m _ is given by

$$[eqn]$$

or for the nearest-neighbor by assuming, $[eqn]$ , we have

$$[eqn]$$ $$[eqn]$$

By construction, the rate matrix satisfies detailed balance, R _ nm _ e ^–βG _ m _ ^ = R _ mn _ e ^–βG _ n _ ^, and probability conservation, ∑_ m _ R _ nm _ = 0. The boundary conditions are reflective, with R 0,–1 = R _ N, N+1_ = 0. More generally, for any lag time Δt and any time t, we have

$$[eqn]$$

In this equation, ** U ** represents a discretized form of transition probabilities, which can be determined empirically. On the other hand, ** R ** is a matrix whose elements are defined by the expression given in Rel. ?. In compliance with the detailed balance condition, ** R ** must be a tridiagonal matrix such that R _ n,n _ = – R _ n,n+1_ – R _ n,n–1_. Now it can be seen that working in a Riemannian CV space is valid for the determination of quantities of interest such as rate, free energy, diffusion, and reaction pathways. However, it should be noted that this kinetic treatment assumes an effective diffusive (Markovian) description in the chosen CVs at an appropriate lag time Δt; if strong memory effects persist at all Δt, a CV diffusion model is not expected to yield quantitatively reliable rates.

The rate matrix R can be estimated self-consistently, ?−? ? by maximizing the likelihood of the observed transitions. Let a trajectory snapshot α start in bin i α at time t and be found in bin j α at time t + Δt. For a continuous-time Markov process the probability of that single observation equals the corresponding element of the time-propagator P(Δt) = exp­(RΔt). Assuming the observations are independent, the total likelihood of the data set is the product of those matrix elements: L = ∏α [e ^ ** R **Δt ^]_ i α, j α _. Therefore, G and h can be found both by maximizing the likelihood L.

From the optimized R, free energy between any two states n and m, ΔG _ nm _ = G _ n _ – G _ m _ can be calculated from the detailed balance relation:

$$[eqn]$$

and the metric at midpoint, $[eqn]$ , can be estimated as

$$[eqn]$$

which indicates the metric at the midpoint, $[eqn]$ , depends on the rate matrix elements and the metric at the adjacent grid points, h _ n , and h _ n+1. To obtain a numerical estimate, one may use this relationship within an iterative scheme or assume that h is slowly varying such that $[eqn]$ . Finally, the MFPT for a transition starting in bin A and reaching bin B along the MFEP, τ_ A→B _, which is the inverse of the rate constant k _ A→B _, can be written in closed form as the standard 1D Smoluchowski result: ?−? ?

$$[eqn]$$

Here r _ A _ and r _ B _ denote the positions of the bin A and B along the path coordinate r, while symbols ρ and σ are dummy integration variables along this path. The invariant PMF G(r) is obtained from the rate matrix (Rel. 38) and differs from the conventional PMF by the geometric term $[eqn]$ . Relation is obtained by first writing the Smoluchowski equation for diffusion in the geodesic coordinate s, where the invariant potential G(s) and the standard 1D MFPT expression apply, and then transforming back to the arbitrary path label r. This change of variables introduces the metric factors $[eqn]$ in the integrand. The inner integral in Rel. ? accumulates the equilibrium weight of all points between basin A and an intermediate position ρ, i.e., how much probability must flow through ρ, while the outer integral weights each intermediate position by its local resistance to motion. In practice, the double integral can be evaluated with standard quadrature techniques.

2-D Toy Model Example

5

To illustrate the practical application of the estimation methods proposed in Section and to compare the conventional and Riemannian frameworks, we examine a 2-D system evolving under an overdamped Langevin equation on the Muller-Brown potential surface.? The potential is defined as

$$[eqn]$$

with a scaling factor α = 0.1. The parameter sets (A _ i _/kcal mol^–1^, a _ i _, b _ i _, c _ i _, x _0i _, y _0i _) are (−200, −1.0, 0.0, −10.0, 1.0, 0.0), (−100, −1.0, 0.0, −10.0, 0.0, 0.5), (−170, – 6.5, 11.0, −6.5, −0.5, 1.5), and (15, 0.7, 0.6, 0.7, −1.0, 1.0) for i = 1, 2, 3, and 4, respectively. This potential is shown in FigureA.

Comparison of thermodynamic and kinetic quantities for a 2-D toy model, analyzed using conventional and Riemannian frameworks. (A) Muller-Brown potential energy surface. The gray line represents the true MFPT connecting the two primary energy minima. (B) Diffusion profiles along the arbitrarily parametrized path r in the log-polar (X, Y) spcae. The black line shows the constant diffusion coefficient D used in the Riemannian framework. The blue line shows the position-dependent diffusion coefficient D(r) = D/h(r). (C) The ground truth free energy profile, G(s), calculated along the true arc length s of the MFEP in the Cartesian (x, y) space. (D) Free energy profiles calculated along the same arbitrarily parametrized path r.

We consider two set of CVs (i) the Cartesian coordinate (x, y) and (ii) a non-Euclidean log-polar space (X, Y) where $[eqn]$ and Y = arctan­(y, x). To construct a 1D model of the system’s thermodynamics and kinetics, the first step is to identify the MFEP. The MFEP in the Cartesian coordinate space, (x, y) was previously identified for this toy model using the SMwST, yielding a coarse-grained pathway of 20 centers? shown in FigureA. The MFEP in the (X, Y) can be identify with Reimannian version of SMwST as discussed in our previous work,? but here our primary goal is to demonstrate that for a given path, even one parametrized by an arbitrary or geometrically incorrect coordinate, the Riemannian analysis framework proposed in Section can recover the thermodynamics and kinetics property the same as conventional one.

To construct a robust 1D model suitable for rate analysis, a higher-resolution path is necessary to satisfy the Markovian assumption. We prepared two distinct high-resolution paths, each discretized into 1000 bins. The first path serves as our “ground truth” and was generated by reparameterizing the initial MFEP by its true Euclidean arc length, s, in the (x, y) space. The second path, designed to test the analysis frameworks, was generated by first transforming the (x, y) coordinates into (X, Y), then reparameterized using an intentionally incorrect Euclidean metric in this space, defining an arbitrary path coordinate denoted as r, such that dr ^2^ = dX ^2^+dY ^2^. We not that the correct distance element in the (X, Y) is $[eqn]$ . Using these paths, we constructed ideal rate matrices for three distinct cases: (1) Using the ground truth path parametrized by the true arc length, s, in the (x, y) space, we built a rate matrix using the true potential, G(s) = V(x, y) and a constant diffusion D = k B T, δs = 1/1000 using

$$[eqn]$$

This case serves as our reference for the correct physical properties; (2) Using the arbitrarily parametrized path r in the (X, Y) space, we built a rate matrix based on the conventional framwork using A(r) = V(X, Y) – 2k B T X and a position-dependent diffusion coefficient, D(r) = D/h(r), where D = k B T and h(r) = exp­(2X), and δr = 1/1000. (3) Finally, using the same arbitrarily parametrized path r from case 2, we built an ideal rate matrix. As this matrix describes the same underlying physics as in case 2, it is constructed to be identical. We then applied the estimation procedure from Section to this rate matrix to determine the Riemannian PMF, G(r), and the 1D metric, h(r), namely via Rel. ? and ?. From the resulting rate matrices, we estimated the free energy and MFPT to compare the outcomes of the three frameworks.

FigureC,D show the calculated free energy profiles. Our analysis reveals that first, for a given, arbitrarily parametrized path in the non-Euclidean (X, Y) space, both the conventional and Riemannian analysis frameworks yield identical G(r) and MFPT of 1.9 × 10^6^ (time unit). These results confirm the mathematical equivalence of the two formalisms: the geometric information encoded in the Riemannian 1D metric h(r) is implicitly captured by the position-dependent diffusion D(r) in the conventional treatment. However, calculated MFPT is still different from the ground truth value of 1.3 × 10^6^ (time unit) obtained from the true MFEP in the Cartesian space. This discrepancy arises because the path analyzed in the (X, Y) space is not the true MFEP, but rather a nonoptimal pathway that results from incorrectly treating the curved space as Euclidean. This demonstrates the central advantage of the Riemannian framework: while both analysis methods are self-consistent, only a Riemannian approach to path-finding can guarantee the identification of the correct, physically meaningful transition path, thereby leading to accurate predictions of kinetic properties like the MFPT.

Biomolecular Example: Alanine

Dipeptide

6

While the toy model above highlights the advantages of the Riemannian framework in a controlled setting, we now demonstrate its use for a biomolecular system in which projected dynamics in CV space may not be strictly Markovian. As an illustrative example, we consider alanine dipeptide (N-acetyl-N′-methylalaninamide, NANMA) in vacuum, which exhibits well-characterized metastability in the (ϕ,ψ) dihedral space (FigureA). We focus on the rare transition between the C 7 ^ eq ^ and C 7 ^ ax ^ basins.

*Biomolecular example (alanine dipeptide): Comparison of thermodynamic and kinetic quantities for the C7

eq → C7

ax transition, analyzed using conventional and Riemannian frameworks. (A) Two-dimensional free-energy surface in the (ϕ, ψ) space; the gray curve shows a transition pathway estimated as MFEP, discretized into 20 equally spaced bins. (B) PMF profiles along the path coordinate (r = 0 and r = 1 correspond to C7

eq and C7

ax , respectively), computed using the conventional definition and the Riemannian definition. (C) Diffusion profiles inferred along the same discretized pathway: the conventional position-dependent diffusion D(r) and the Riemannian effective diffusion D/h(r). (D) The corresponding one-dimensional metric h(r) estimated along the pathway.*

A pathway connecting these basins was first identified using the string method with swarms of trajectories.? Treating this pathway as an approximation to the MFEP in the chosen CV space, we discretized it into a set of centers and estimated a lag-time dependent rate matrix, R _Δt _, from short unbiased MD trajectories initiated from these centers (see Supporting Information (SI) for simulation details). Because memory effects can be present after projection to CV space, the resulting kinetic model generally depends on the chosen lag time. Here, to illustrate the method and enable a controlled comparison between the conventional and Riemannian analyses, we select Δt = 10 fs, for which we hypothesize that the relevant slow modes have stabilized, consistent with an effective Markovian description on the discretized path. In general, the reliability of such a kinetic model can be assessed by monitoring the stability of implied time scales or metric over a range of lag times; for the present system, results for different lag time provided in Figure S1.

The free-energy profiles obtained from the conventional definition and from the invariant Riemannian definition are shown in FigureB and are nearly identical. This indicates that, for this particular transition pathway, the impact of the Riemannian framework on the PMF is small relative to the barrier height (or is effectively absorbed when the path is well optimized); therefore, the conventional and Riemannian PMFs need not differ substantially in all practical cases. The distinction emerges more clearly in the dynamical interpretation. In the conventional framework, variations in the estimated transition rates along the path are commonly interpreted as a position-dependent diffusion coefficient, D(r) (FigureC). In the Riemannian framework, the same variations are instead attributed to the geometric metric h(r), which quantifies local stretching (or resistance) along the path. Consistently, plotting the effective diffusion predicted by the Riemannian model, D(r)∝1/h(r), reproduces the shape of the conventional diffusion profile.

Finally, we computed the MFPT for the C 7 ^ eq ^ → C 7 ^ ax ^ transition. The conventional and Riemannian analyses yield comparable time scales (approximately 110 and 90 ns, respectively), consistent with their mathematical equivalence for a fixed discretized pathway. Importantly, the Riemannian decomposition provides a clearer physical interpretation by separating geometric effects (encoded in h) from energetics (encoded in G): what appears as variable diffusion in the conventional treatment can be understood as a consequence of curvature and stretching in CV space along the transition path.

Summary

and Conclusions

7

In this study, we presented a generalized Riemannian framework that rigorously addresses the noninvariance issues inherent to conventional definitions of thermodynamic and kinetic properties in CV spaces. By using Riemannian geometry, we introduced robust definitions of PMF, MFEP, and diffusion models that remain invariant under nonlinear coordinate transformations. Through detailed theoretical derivations, we demonstrated how conventional Euclidean assumptions can introduce coordinate-dependent artifacts, highlighting the necessity of a proper geometric treatment.

The central thermodynamic result is the relationship between the conventional PMF Ã(ζ) and the invariant (Riemannian) PMF A(ζ), which differ by an explicit metric correction term (eqs and ?). This expression provides a direct estimate of the potential magnitude of coordinate-induced distortions: when the metric varies weakly over the relevant region of CV space, conventional PMFs are often practically reliable up to an (approximately) constant offset, whereas stronger metric variations can lead to non-negligible changes in the apparent profile and barrier heights. On the kinetic side, we formulated the effective diffusive description in CV space in a coordinate-consistent manner (eqs and ?). Building on this, we presented a practical path-based strategy in which the dynamics are discretized along a pathway and represented by a nearest-neighbor rate model, from which free-energy differences and MFPTs can be computed (eqs, ?, ?, ?, ?, and ?).

Using 1D and 2-D representative examples, we explicitly illustrated the practical consequences of neglecting CV space geometry, revealing that a correct Riemannian approach accurately reproduces invariant free energy landscapes and consistent kinetic predictions, whereas conventional methods may yield biased estimates, depending on the choice (and parametrization) of the CVs. This framework thus not only resolves fundamental theoretical inconsistencies but also significantly enhances the interpretability and reliability of biomolecular simulations. Future work will explore the application of this Riemannian formalism to more complex biomolecular systems to further validate its broad applicability.

Supplementary Material