Predicting reaction coordinates in energy landscapes with diffusion anisotropy

TL;DR

This paper uses the SGOOP method to determine optimal reaction coordinates in energy landscapes with diffusion anisotropy, incorporating dynamical and static information for improved accuracy.

Contribution

It extends the SGOOP method to include dynamical observables, enabling accurate reaction coordinate prediction under diffusion anisotropy.

Findings

Incorporates dynamical observables into SGOOP

Accurately predicts reaction coordinates with anisotropic diffusion

Validates results against transmission coefficient calculations

Abstract

We consider a range of model potentials with metastable states undergoing molecular dynamics coupled to a thermal bath in the high friction regime, and consider how the optimal reaction coordinate depends on the diffusion anisotropy. For this we use our recently proposed method 'Spectral gap optimization of order parameters (SGOOP)' (Tiwary and Berne, Proc. Natl. Acad. Sci. 113 2839 2016). We show how available information about dynamical observables in addition to static information can be incorporated into SGOOP, which can then be used to accurately determine the 'best' reaction coordinate for arbitrary anisotropies. We compare our results with transmission coefficient calculations and published benchmarks where applicable or available respectively.

Figures

Equations

$$e^{-\beta F(q)}=\frac{\int \int dxdy\delta (q-q(x,y))e^{-\beta U(x,y)}}{\int \int dxdy{e}^{-\beta U(x,y)}}$$

$$\frac{\partial p(q,t)}{\partial t}=\frac{\partial }{\partial q}D(q)e^{-\beta F(q)}\frac{\partial }{\partial q}{e}^{\beta F(q)}p(q,t)\equiv -Hp$$

$$p(q,t)=p^0(q)+{\Sigma }_{i=1}^{\infty }{c}_i{\varphi }_i(q)e^{-{\lambda }_{i}t}$$

$$\left(\begin{array}[]{cc}D_{x}&0\\ 0&D_{y}\end{array}\right)$$

$$\frac{\partial p(q,t)}{\partial t}=D_{\theta }\frac{\partial }{\partial q}{e}^{-\beta F(q)}\frac{\partial }{\partial q}{e}^{\beta F(q)}p(q,t)$$

$$D_{\theta }=D_{x}co{s}^2(\theta )+D_{y}si{n}^2(\theta )$$

$$\frac{\Delta p_n(t)}{\Delta t}={\Sigma }_m{k}_{mn}{p}_m(t)-{\Sigma }_m{k}_{nm}{p}_n(t)\equiv {\Sigma }_m{\mathbf{K}}_{nm}{p}_m(t)$$

$$k_{mn}=\frac{D_{\theta }}{2d^2}\sqrt{\frac{p_n^0}{p_m^0}}$$

$$⟨N⟩=\sum _{\begin{array}{c}(m,n)\\ \forall |m-n|=1\end{array}}{p}_m^0{k}_{mn}{N}_{mn}=\sum _{\begin{array}{c}(m,n)\\ \forall |m-n|=1\end{array}}{p}_m^0{k}_{mn}$$

$$k_{mn}=\lambda \sqrt{\frac{p_n^0}{p_m^0}}=\left\{\frac{⟨N⟩}{\sum \sqrt{p_n^0{p}_m^0}}\right\}\sqrt{\frac{p_n^0}{p_m^0}}$$

$$U(x,y)=$$

$$-18(e^{-2(x-.8{)}^2-2(y-1.2{)}^2})$$

$$-16(e^{-2(x-.5{)}^2-2(y+.3{)}^2})+0.5(x^6+y^6)$$

$$\left(\begin{array}[]{cc}D_{x}&0\\ 0&\delta D_{x}\end{array}\right)$$

Full text

I Introduction

The notion that the potential energy surface (PES), or often the potential of mean force (PMF), governs the choice and the form of the reaction pathway is a central tenet in chemical physics. Given the PES (PMF), numerous very effective and successful techniques existHenkelman, Uberuaga, and Jónsson (2000); Branduardi, Gervasio, and Parrinello (2007); Weinan, Ren, and Vanden-Eijnden (2007) that can identify the minimum energy (free energy) pathway between a collection of metastable states, which is then taken to be the reaction pathway. The reaction coordinate (RC) is proverbially defined as an abstract low-dimensional coordinate that best captures progress along this reaction pathway. We define the RC more rigorously in Sec. II.1.

However, such a reaction pathway and reaction coordinate derived from a purely static framework can be significantly perturbed by the stochastic nature of the environment, and specifically by how the various degrees of freedom couple to the bath and relax in response to perturbations in it van der Zwan and Hynes (1982); Berezhkovskii and Szabo (2005); Hynes (1985); Kl/osek-Dygas et al. (1989); Johnson and Hummer (2012); Peters et al. (2013); Echeverria, Makarov, and Papoian (2014). This coupling can be quantified through either a friction coefficient, or through its inverse, the diffusion constant. In the limit of strong coupling to the environment, one then expects that the diffusion will differ for different degrees of freedom, this anisotropy will manifest itself in the reaction pathway. A certain reactive mode derived from the PES/PMF might no longer be the path of least resistance between two metastable states, and one needs to consider the work done not just against the static forces from PES/PMF, but the dynamic effect of the environment. In fact, several workers have pointed out that the effect of diffusion anisotropy can be so strong as to lead to the phenomenon of “saddle point avoidance”Northrup and McCammon (1983); Berezhkovskii, Berezhkovskii, and Zitzerman (1989), where the saddle points between metastable states could altogether be avoided while moving between them.

Especially when one deals with free energies (as contrasted to potential energies) which are often defined as functions of qualitatively different collective variables, the diffusivities for the various coordinates will not be the same, and could possibly even be position and temperature dependent. For instance, for the very important problem of protein-ligand (un)bindingCopeland (2016); Tiwary et al. (2015a); Casasnovas et al. (2017), these variables could be ligand-host relative displacement, their individual conformations, their hydration states and other variables Tiwary and Berne (2016a); Agmon and Rabinovich (1992). These variables can naturally be expected to possess very different diffusivities. For example, at low temperatures the protein fluctuations can be completely frozen out, with the protein effectively trapped in one metastable conformation, while ligand unbinding still happensSteinbach et al. (1991). In such cases, the free energy barrier is thus not the only determinant of the dynamics, and the RC which encapsulates the dominant slow fluctuations in the system will thus be a function of the diffusion anisotropy. In fact, while the free energy can be calculated as a function of any arbitrary order parameter, it is the most informative when expressed as a function of the optimal RC. This work is restricted to simple model potentials, as we are mainly interested in understanding how the optimal RC depends on diffusional anisotropy—a question pertinent to a range of very practical problems involving for example protein conformational transitions, drug un-binding and the efficiency of enhanced sampling.Tiwary and Berne (2016b, a); Tiwary et al. (2015b)

In the current work, we tackle the problem of how diffusion anisotropy affects the RC for different model potential energy landscapes with high and arbitrary number of energy barriers. To identify the RC and its dependence on diffusion anisotropy, we propose a practical numerical solution using our recent method “Spectral gap optimization of order parameters (SGOOP)”.Tiwary and Berne (2016b) For several potential surfaces we demonstrate how dynamical information can be incorporated into SGOOP which can then be used to predict the RC for arbitrary diffusion anisotropies. We look at three 2-d model potentials undergoing Langevin dynamics with different values of the diffusion anisotropy, and predict how the RC changes with anisotropy for these different cases. Where possible, we validate the SGOOP-predicted RC with published results which use different methods for RC calculationBerezhkovskii and Szabo (2005), and also with extensive calculations of the reactive flux and the transmission coefficients along different putative RCs. The method SGOOP, which is described in Sec. II, considers the inverse problem of what can we infer about the RC given limited information about the static and dynamic behavior of the system. The dynamic information in SGOOP can be explicitly implemented if one knows the diffusion tensor, or through a Maximum Caliber frameworkPressé et al. (2013); Dixit et al. (2015), where one effectively fits the response function to the environment depending on known dynamical observables. These dynamical observables are clarified in Sec. II.5. We use both the methods and find excellent agreement between the results for all three potentials. Given the robustness and reliability of SGOOP shown in this work, we expect it will be very useful for complex systems with qualitatively different collective variables that respond differently to the environment.

II Theory

II.1 Reaction coordinate

In spite, and possibly due to, its ubiquity across the chemical kinetics literature, the reaction coordinate (RC) has multiple definitions, often depending on the intended use. Here we clarify our interpretation of the RC. For a given multidimensional complex system undergoing a certain dynamics, we define RC as a low-dimensional (say one or two dimensional) variable such that the multidimensional dynamics of the full system in terms of movement between different metastable states, can be “best” mapped into Markovian dynamics between various metastable states viewed as a function of the low-dimensional RC. In other words, the RC should have two features: (1) it should be able to demarcate between the various important metastable states that the full system possesses, and (2) there should be a separation of timescales between the relaxation times in the various metastable states, and the relatively slow transition times between them. In Section II.3 we define the timescale separation more carefully through the notion of the spectral gapTiwary and Berne (2016b); Rohrdanz et al. (2011); Coifman et al. (2008a, b).

II.2 Potential of mean force along putative reaction coordinate

Without loss of generality, we consider a 2-dimensional energy landscape, in the space of Cartesian or generalized variables $x$ and $y$, with corresponding PES or PMF respectively given by $U(x,y)$. On this landscape we consider a putative RC $\hat{u}=\hat{x}cos(\theta)+\hat{y}sin(\theta)$, oriented at an angle $\theta$ from $\hat{x}$. For example, in Fig. 1(a) the magenta colored dashed line denotes a putative RC rotated by an angle $\theta$ from the x-axis (solid black line). Note that the relevant parameter for $\hat{u}$ parametrized as such is its direction, and it is invariant to translation. This could easily be generalized to more than 2 variables with the use of more than 1 angle. The PMF $F(q)$ where $q$ measures distance along the unit vector $\hat{u}$ i.e. $\vec{q}=q\hat{u}$ can be defined as

[TABLE]

where $\beta={(k_{B}T)^{-1}}$ is the inverse temperature. The normalization in the denominator allows us to define an associated stationary probability density $p^{0}(q)=e^{-\beta F(q)}$.

The PMF can be calculated analytically/numerically for different choices of $\hat{u}$ for the model potentials in this work. For more complex higher dimensional potentials it can be calculated through a range of sampling techniques as was shown in Ref.Tiwary and Berne (2016b, a). In this work, the primary interest is to see how the dynamical effects of diffusion anisotropy can be incorporated into SGOOP, and as such as far as the PMF is concerned, we evaluate it analytically/numerically. We found no perceivable differences in the results by using a numerically estimated PMF. However, as explained later, we calculate dynamical observables analytically as well as through sampling.

II.3 Dynamics and spectral gap along putative reaction coordinate

Along a putative RC $\hat{u}=\hat{x}cos(\theta)+\hat{y}sin(\theta)$, we assume that the time-dependent probability density $p(q,t)$ is governed by a Smoluchowski Equation with position dependent diffusivity:

[TABLE]

Under fairly general conditions, the operator $H$ has a discrete spectrum of non-negative eigenvalues $\{\lambda_{i}\}$ with $\lambda_{0}=0<\lambda_{1}\leq\lambda_{2}\leq...$, and corresponding eigenvectors $\phi_{i}(q)$, and the general solution of Eq. 2 is given by

[TABLE]

For a system with $n\geq 1$ slow processes (for example corresponding to existence of $n$ metastable states), the eigenspectrum of the operator $H$ will have a gap between the eigenvalues $\lambda_{n}$ and $\lambda_{n+1}$, i.e. $\lambda_{n}\ll\lambda_{n+1}$. Please note that strictly speaking we consider the spectrum of eigenvalues of the associated master equation described in the following section which has similar properties. We refer to the quantity $e^{-\lambda_{n}}-e^{-\lambda_{n+1}}$ as the spectral gap. The larger the spectral gap, the larger will be the timescale separation between the slow and fast processes when the dynamics is viewed as a function of the putative RC $\hat{u}$.

We consider a general diffusivity tensor in the Cartesian basis given by

[TABLE]

i.e. one without any off-diagonal terms. Interesting complications would naturally arise if these cross-terms were not 0, but they can always be dealt with by diagonalizing the diffusion tensor and working with a different PES/PMF. We define the anisotropy parameter as $\delta={D_{y}/D_{x}}$. We consider different values of $\delta$ ranging from $\delta\ll 1$ to $\delta\gg 1$, where for instance $\delta\ll 1$ represents the case when diffusion along $y$-axis is much slower than diffusion along the $x$-axis.

For simplicity, we also assume that the terms $D_{x}$ and $D_{y}$ have no positional dependence. Under this assumption, the diffusivity $D_{\theta}$ along any putative RC $\hat{u}$ rotated at angle $\theta$ from the $x-$axis can be removed from inside the partial derivative in Eq. 2, which can now be written as:

[TABLE]

In fact for the two-dimensional case given the diffusivity tensor above we can easily calculate $D_{\theta}$:

[TABLE]

II.4 Master equation along putative reaction coordinate

We now spatially discretize the putative RC $\hat{u}$ by defining a grid $\{q_{n}\}$ along it, where $n$ takes integral values. Let $p_{n}(t)$ denote the instantaneous probability of being on any of these grid points $n$ at time $t$ (also discretized in intervals $\Delta t$), and $k_{mn}$ be the time-independent rate of transition from grid point $m$ to $n$ per unit time $\Delta t$. We can then write a master equation governing the flow of probability from one grid point to another:

[TABLE]

The matrix $\bf K$, where ${\bf K}_{nm}=k_{mn}$, is the entirety of all these rates. The diagonal terms in this matrix are obtained through $\Sigma_{n}k_{mn}=0$. Under very general conditions the master equation possesses a unique stationary state solution, which we denote $p^{0}_{n}$ in accordance with Sec. II.2. We further assume that detailed balance is satisfied i.e. $p^{0}_{n}k_{nm}=p^{0}_{m}k_{mn}$ which can be shown to be true for closed isolated physical systems, but is in general a decent assumption.

Under these conditions, the master equation Eq.6 will also have a solution through expansion in eigenfunctions similar to Eq.3. For convenience we use the same notation for the eigenfunction solution to the Smoluchowski equation as well as the master equation, and denote the eigenspectrum as $\{\lambda_{i}\}$ with $\lambda_{0}=0<\lambda_{1}\leq\lambda_{2}\leq...$. It is the spectral gap arising from this eigenspectrum that we use in SGOOP and in this work.

II.5 Relations between transition rates and stationary probability

The transition rates $k_{mn}$ are underconstrained by the information presented so far including the detailed balance condition. We now introduce further relations connecting the transition rates $k_{mn}$ to the stationary probabilities $p^{0}_{n}$. Such a relation is fundamental to SGOOP, which we will summarize in the next section Sec. II.6. It can be derived in 2 ways which we describe here. In this work we use both these relations and assess how SGOOP performs with each of them.

The first way to connect the transition rates $k_{mn}$ to the stationary probabilities $p^{0}_{n}$ is by discretizing the Smoluchowski Equation (Eq. 2) and then comparing it with the master equation (Eq. 6). This is described in Ref. Bicout and Szabo, 1998 and below we give the final result:

[TABLE]

where $d=|u_{k}-u_{k-1}|$ $\forall k$ is the grid spacing, $D_{\theta}$ is the diffusivity along the putative RC $\hat{u}$, which as explained in Sec. II.3 we assume to depend on $\theta$ but not vary along $\hat{u}$. This approach is practical if one knows beforehand or calculates the diffusivity along various degrees of freedom for instance by using the approach of Ref. Straub, Berne, and Roux (1990); Hummer (2005), and then uses Eq. 5.

A second approach connecting the transition rates $k_{mn}$ to the stationary probabilities $p^{0}_{n}$ is through the Maximum Caliber framework. While fairly general treatments are possible in this framework, here we start by assuming that apart from the stationary state information all one knows is the average value $\langle N\rangle$ of some dynamical observable $N$ over the full phase space. For instance $N$ could be the average number of first-nearest neighbor transitions observed in time $\Delta t$. Then $\langle N\rangle$ is given by

[TABLE]

where we have used $N_{mn}=1$ $\forall|m-n|=1$ and [math] otherwise. We then seek an equation akin to Eq. 7 with the only known constraints being (1) principle of detailed balance (2) measurement of $\langle N\rangle$. The principle of Maximum Caliber maximizes the path entropy subject to these constraintsPressé et al. (2013); Dixit et al. (2015), and yields the same square-root of probability dependence, but now instead of using the explicit diffusivity, one works in the framework of Lagrange multipliers:

[TABLE]

where we have calculated $\lambda$ by using the first equality of Eq. 9 in Eq. 8. Note that similar expressions could be derived using a very wide range of dynamical observables and not just the one we use here.

With either Eq. 7 or 9, accompanied with either explicit diffusivities or measurement of $\langle N\rangle$ using molecular dynamics (MD) respectively, we can complement our knowledge of the stationary probability $p_{0}$ to calculate the full transition matrix $\bf K$, where ${\bf K}_{nm}=k_{mn}$ as defined above. The eigenvalues of this matrix directly give us the spectral gap.

A technical but crucial point must be mentioned here. In order to define and calculate $\langle N\rangle$ we have limited ourselves to first neighbors. We could go further than that, but then we subject ourselves to noisy and poorly sampled observations (i.e. transitions between far off grid points on $\{u_{k}\}$). As such it is desirable to stick to first or first few neighbors in this context. However when we evaluate the eigenspectrum of the matrix $\bf K$,which could have come either through Eq. 7 or Eq. 9, we are faced with a similar problem - how sparse should the matrix $\bf K$ be? Here, in contrast to the calculation of $\langle N\rangle$, it is desirable to have not just first neighbor terms as non-0 but ideally terms in the whole matrix. From practical perspectives for both these operations we start with first neighbor, move on to second neighbors, and so on, and stop when the spectral gap does not change.

We finish this subsection by making the interesting observation that equations very similar to Eqs. 7 and 9 have recently been used in the context of machine learning as well Sohl-Dickstein, Battaglino, and DeWeese (2011, 2009).

II.6 Summarizing Spectral gap optimization of order parameters (SGOOP)

Given the full theoretical framework of the preceding sections, it is now easy to describe SGOOP Tiwary and Berne (2016b), which is a method to construct an optimal RC as a combination (linear or non-linear) of a set of many candidate order parameters or collective variables $\Psi=(\Psi_{1},\Psi_{2},...,\Psi_{d})$. This set is assumed to be known a priori. The working definition of the optimal RC that SGOOP follows was described in Secs. II.1 and II.3: (1) demarcate between various metastable states that the full system possesses, and (2) maximize the spectral gap of the eigenvalues of the associated master equation. Once again, let $\{\lambda\}$ denote this set of eigenvalues, with $\lambda_{0}=0<\lambda_{1}\leq\lambda_{2}\leq...$. We refer to the quantity $e^{-\lambda_{n}}-e^{-\lambda_{n+1}}$ as the spectral gap, where $n$ is the number of barriers apparent from the free energy estimate that are higher than a user-defined threshold (typically $\gtrsim k_{B}T$). This calculation is done without any pre-knowledge of the number of metastable states the system might have, and as such makes it easy to use SGOOP for systems with an arbitrary number of metastable states.

The first input to SGOOP is an estimate of the stationary probability density (or equivalently the free energy) of the system as a function of any putative RC $f(\Psi)$. For the model systems in this work $\Psi=\{x,y\}$ and the stationary density for any putative RC can be calculated analytically or numerically. For more complex systems Tiwary and Berne (2016b, a) SGOOP uses biased simulationTiwary and Parrinello (2014) performed along a sub-optimal trial RC given by some linear or non-linear function $f_{0}(\Psi)$, and then a posteriori reweighting it along any putative RC $f(\Psi)$. The ability to perform computationally easy a posteriori reweighting without repeating the simulations is crucial. The second input to SGOOP is some form of dynamical information. Depending on whether Eq. 7 or Eq. 9 is used, this would respectively be the diffusivity tensor from Sec. II.2, or the average number of transitions $\langle N\rangle$ from Sec. II.5. Both these pieces of information are then used through either Eq. 7 or Eq. 9 to obtain the eigenspectrum and hence the spectral gaps of the various trial CVs $f(\Psi)$. Through a post-processing optimization procedure we then find the optimal RC as the $f(\Psi)$ which gives the maximal spectral gap.

For the simple case with just one barrier in two dimensions, our method amounts to maximizing the difference between the first two nonzero eigenvalues of the master equation, which is in the spirt of other methods such as Ref. Berezhkovskii and Szabo (2013). In this work we implement SGOOP in two ways - by explicitly considering the diffusion tensor in the discretized Smoluchowski equation and through a Maximal Caliber framework, and demonstrate excellent agreement between both. The novelty in SGOOP partly comes from its ability to deal with more than one barrier which we show explicitly in the Results section, as well as through its use of the Maximal Caliber framework when the diffusion tensor is not known.

III Results

Previous applications of SGOOP have been limited to considering only the stationary density information and ignoring dynamical effects. In this work we consider three model 2-dimensional potentials undergoing Langevin dynamics coupled to a bath, with distinct dynamical effects due to diffusion anisotropy, and demonstrate how SGOOP can be used to accurately predict the effects of diffusion anisotropy on the RC. The first two are 2-state potentials, namely the Berezhkovskii-Szabo (BS) potential Berezhkovskii and Szabo (2005) and the Deleon-Berne (DB) potential De Leon and Berne (1981). The third potential is a 3-state potential which we describe later. The set of candidate order parameters is given by $\Psi=\{x,y\}$, and for each potential, we identify the optimal RC as a linear combination $f(x,y)=ax+by$ for different diffusion anisotropies $\delta={D_{y}/D_{x}}$. The This amounts to calculating the optimal angle $\theta$ from $\hat{x}$ for RC given by $\hat{u}=\hat{x}cos(\theta)+\hat{y}sin(\theta)$ for $\delta=$ 0.1, 1 and 10. For all the potentials, Newton’s laws of motion were integrated per Langevin dynamicsBussi and Parrinello (2007) with a timestep 0.01 unit and different values of anisotropic friction coefficients $\gamma_{x},\gamma_{y}$. Since $\delta={D_{y}/D_{x}}$, we have the relation $\delta={\gamma_{x}/\gamma_{y}}$. Corresponding to $\delta=$ 0.1, 1 and 10, we took ($\gamma_{x},\gamma_{y}$) = (2,20), (5,5) and (20,2) respectively in inverse time units. All of these simulations correspond to the overdamped friction regime Tiwary and Berne (2016c).

We consider two different 2-state potentials because of the structure of the diffusion tensor. In one of them (DB), the two directions $x$ and $y$ are parallel/perpendicular to the line through the saddle connecting the two minima, while in the other (BS), both the directions are at a tilt with respect to this line. This simple difference, as we show later, leads to an arguably profound difference in how the RC depends on diffusion anisotropy. For the 2-state potentials we also perform explicit calculations of the transmission coefficient given by the ratio of the rate constant $k_{MD}$ from long MD simulations, and the TST rate constant $k_{TST}^{\theta}$ along various $\hat{u}(\theta)$.

III.1 Berezhkovskii-Szabo (BS) potential

The first potential we consider is a 2-state potential (Fig. 1(a)) introduced by Berezhkovskii and Szabo Berezhkovskii and Szabo (2005). We start with this potential because the authors in Ref. Berezhkovskii and Szabo, 2005 have explicitly calculated the dependence of the RC on the anisotropy parameter $\delta$ through a different method. Thus it serves as an excellent benchmark for our work. In Fig. 1(a) the magenta colored dashed line denotes a putative RC $\hat{u}=\hat{x}cos(\theta)+\hat{y}sin(\theta)$, at an angle $\theta$ from the x-axis (solid black line). The potential is detailed in Ref. Berezhkovskii and Szabo, 2005. Here we perform Langevin dynamics at inverse temperature $\beta=1$.

In Fig. 1 we give the detailed results for this potential. Fig. 1(b) shows the results for the isotropic diffusion case. Here the use of SGOOP with Eq. 7 or Eq. 9 gives the same optimal RC ($\theta$=32∘) in very good agreement with the benchmark calculation in Ref. Berezhkovskii and Szabo (2005). We also provide here our calculations of the transmission coefficients along various putative RCs $\hat{u}(\theta)$, given by the ratio of the rate constant $k_{MD}$, and the TST rate constant $k_{TST}^{\theta}$ along $\hat{u}(\theta)$. It is very interesting that the scaled spectral gap through SGOOP follows roughly the dependence of the transmission coefficient on $\theta$. Fig. 1(c) gives the same information as (b) but for the anisotropic cases. Here as well the agreement with the benchmark calculation in Ref. Berezhkovskii and Szabo (2005) is excellent, with the optimal RC at $\theta$=5∘ and $\theta$=70∘ respectively for $\delta$=10 and 0.1 respectively. Naturally, the transmission coefficient ($\kappa=k_{MD}/k_{TST}^{\theta}$) would have the same qualitative behavior including location of maximal value at $\theta$=32∘, irrespective of the extent of anisotropy. This is because $k_{TST}^{\theta}$ does not depend on the anisotropy in $D$ but only on the static details through the PMF for the putative RC corresponding to $\theta$, and $k_{MD}$ is a measurement independent of the choice of RC i.e. $\theta$. This leads to the possibly obvious but important observation that while maximizing the value of the transmission coefficient still gives the best TST description (i.e. variational transition state theory), it might not give the best reaction coordinate as it does not take into account dynamical effects through the anisotropy of the diffusion coefficients. Finally, in Fig. 1 (d), we provide the respective 1-d PMFs along three putative RCs. It is interesting to see here how the effect of a smaller barrier - which would be ideal for a 2-state problem from the purpose of optimizing transmission coefficient - is compensated by the slower diffusion along a certain direction.

III.2 Deleon-Berne (DB) potential

The second potential we consider is a 2-state potential (Fig. 2(a)) introduced by Deleon and Berne De Leon and Berne (1981). The potential is detailed in Ref. Tiwary and Berne, 2016c. Here we perform Langevin dynamics at inverse temperature $\beta=5$. In Fig. 2 we give the detailed results for this potential. Fig. 2(b) shows the results for the isotropic diffusion case. Here the use of SGOOP with Eq. 7 or Eq. 9 gives the same optimal RC ($\theta$=$90^{\circ}$) in very good agreement with the transmission coefficient calculation. It is again very interesting that the scaled spectral gap through SGOOP follows roughly the dependence of the transmission coefficient on $\theta$. Fig. 2(c) gives the same information as (b) for the anisotropic cases.

For this potential, unlike the Berezhkovskii-Szabo potential, changing the diffusion anisotropy does not change the precise location of the optimal RC. However, it does change the tolerance in quality of RC one has in terms of deviating from the optimal location. Thus as can be seen in Fig. 2(d), the spectral gap versus $\theta$ profiles for the $\delta$=10 and 0.1 cases are broader and narrower respectively than the $\delta$=1 profile in Fig. 2(b). Thus when the diffusion along $y$ is slower ($\delta=$0.1), the penalty for deviating from the optimal RC $\theta$ = $90^{\circ}$ is much higher, and correspondingly it is much lower in the case when diffusion along $x$ is slower ($\delta=$10). We further discuss this point in Sec. IV. In Fig. 2 (d), we provide the respective 1-d PMFs along three putative RCs. As can be seen here, the RC with $\theta$ = 90*∘* gives the highest barrier and thus corresponds to the optimal RC for the isotropic diffusion case.

One interesting observation that can be made from the two 2-state problems we looked at pertains to the qualitative difference in their response to the diffusion anisotropy. The Berezhkovskii-Szabo (BS) potential showed a very clear dependence of the optimal RC on extent of diffusion anisotropy, while the Deleon-Berne (DB) potential had the same optimal RC location irrespective of the extent of anisotropy. Apart from the difference in the precise functional form, both of these are simple 2-state potentials, fairly similar to each other except that in the BS potential, the diffusion tensor non-zero components $x$ and $y$ are not aligned along the normal mode of the PES. In the DB potential, one of the diffusion tensor non-zero components ($y$) is aligned along the normal mode of the PES.

III.3 3-state potential

The final potential we consider is a 3-state potential given by the following functional form:

[TABLE]

Here we perform Langevin dynamics at inverse temperature $\beta=1.25$. In Fig. 3 we give the detailed results for this potential. Fig. 3(b) shows the results for the isotropic diffusion case. Here the use of SGOOP with Eq. 7 or Eq. 9 gives the same optimal RC ($\theta$=81.5∘).

Note that there are two other peaks in this profile, signifying the presence of multiple barriers and thus different families of possible RC. These peaks are however not as dominant as the optimal RC. Fig. 3(c) gives the same information as (b) for the anisotropic cases. For this potential, when we change the diffusion anisotropy to $\delta=$10 (i.e. diffusion along $x$ is slower), we see a flip in the relative strengths of the peaks at $\theta$=81.5∘ and $\theta$=21.8∘, and the optimal RC switches to $\theta$=21.8∘. For the opposite end when $\delta=$0.1 (i.e. diffusion along $y$ is slower), the optimal RC is same as in the isotropic case. The same results are found through the use of Eq. 7 or Eq. 9.

In Fig. 3 (d), we provide the respective 1-d PMFs along the two putative RCs obtained for $\theta$=81.5∘ (blue line) and $\theta$=21.8∘ (red line). In contrast to the other two potentials, now we have two barriers. Similar interpretations can be made here as for the 2-state potentials. The RC for $\delta$=1 has higher barriers (corresponding to both the eigenvalues $\lambda_{1}$ and $\lambda_{2}$). However for the $\delta$=10 case it is the effect of having slower diffusion along $x$ that leads to the putative RC with lower free energy barriers becoming the optimal RC.

III.4 Additional treatment of the Berezhkovskii-Szabo (BS) potential

In this subsection, we revisit the BS potential from yet another perspective which has been used by workers including Hynes, Szabo and others van der Zwan and Hynes (1982); Berezhkovskii and Szabo (2005), and has been described in those papers. Unlike SGOOP, this treatment is applicable only if (1) the reactants and products are separated by a single saddle point, and (2) the dynamics is in the high friction limit. The key idea here is to consider the matrix $DV$ where $D$ is the full diffusion tensor, and $V$ is matrix of second derivatives of the multi-dimensional potential at its saddle point, with det $V<$ 0. Both these matrices could in general have off-diagonal cross terms. By then solving for the only negative eigenvalue and associate eigenvector of the operator $DV$, we get the reaction coordinate for the given potential. We do this exercise for various diffusion anisotropies as quantified by the value of the parameter $\delta$ in the diffusion tensor

[TABLE]

In Fig. 4 we provide the rotation of the optimal RC from $x-$axis for the BS potential as obtained by diagonalizing the $DV$ matrix and calculating the orientation of the eigenvector with the only negative eigenvalue. In this figure we have also overlaid the estimate of the optimal RC from SGOOP and the two estimates are in excellent agreement.

IV Discussion

In this work we have demonstrated how using our recent method SGOOP Tiwary and Berne (2016b) we can reliably and accurately calculate the effect of diffusion anisotropy on the reaction coordinate for different model potential energy landscapes with high and arbitrarily many energy barriers, and arbitrary diffusion anisotropies. Where possible, we validated the SGOOP-predicted RC with published results using different methods for RC calculationBerezhkovskii and Szabo (2005), and also with extensive calculations of the transmission coefficients along different putative RCs. We showed how the dynamic information in SGOOP can be explicitly implemented either through explicit knowledge of the diffusion tensor, or through a Maximum Caliber frameworkPressé et al. (2013); Dixit et al. (2015). We found excellent agreement between the results for all three potentials using both these methods. Unlike other methods for optimizing RC that we are aware of, SGOOP was also shown to give an approximate quantification of how tolerant is the system with respect to deviations from the optimal RC.

Clearly, knowing the optimal RC is important from the perspective of kinetics calculations. However we believe that diffusion anisotropy and how it is manifested in the interplay between statics and dynamics might be an overlooked but vital problem in free energy sampling as well (barring few notable exceptions such as Ref. Poon and Peters (2015)), where one is often faced with a choice of many qualitatively different order parameters or collective variables Valsson, Tiwary, and Parrinello (2016); Tiwary and van de Walle (2016); McCarty and Parrinello (2017). Barring very specialized conditions, these variables will naturally have different diffusivities. In such situations, a method like SGOOP should be very useful for accurate and efficient free energy sampling. A free energy profile sampled along the optimal RC will be easier to compute through umbrella sampling or metadynamics, and features in this profile will be of more direct use as well for drawing physical conclusions. This will be the subject of future investigations probing the thermodynamics and kinetics of complex systems.

**ACKNOWLEDGMENTS

**

This work was supported by grants from the National Institutes of Health [NIH-GM4330] and the Extreme Science and Engineering Discovery Environment (XSEDE) [TG-MCA08X002]. The authors thanks Purushottam Dixit for helpful discussions.