A variational conformational dynamics approach to the selection of collective variables in metadynamics

TL;DR

This paper introduces a novel method combining metadynamics and time-lagged independent component analysis to improve the selection of collective variables, enhancing the study of rare events and complex free energy landscapes.

Contribution

The paper presents a new variational conformational dynamics approach that optimizes collective variables for metadynamics using conformational analysis techniques.

Findings

More efficient collective variables are obtained.

Improved convergence in free energy landscape exploration.

Enhanced understanding of system physics.

Abstract

In this paper we combine two powerful computational techniques, well-tempered metadynamics and time lagged independent component analysis. The aim is to develop a new tool for studying rare events and exploring complex free energy landscapes. Metadynamics is a well-established and widely used enhanced sampling method whose efficiency depends on an appropriate choice of collective variables. Often the initial choice is not optimal leading to slow convergence. However by analyzing the dynamics generated in one such a run with a time-lagged independent component analysis and the techniques recently developed in the area of conformational dynamics, we obtain much more efficient collective variables, that are also better capable of illuminating the physics of the system. We demonstrate the power of this approach in two paradigmatic examples.

Figures

Equations

$$\underset{t\to \infty }{\lim }V(s,t)=-\left(1-\frac{1}{\gamma }\right)F(s)$$

$$F(s)=-\frac{1}{\beta }\log \int dR\delta (s-s(R))e^{-\beta U(R)}$$

$$⟨O(R)⟩=\underset{T\to \infty }{\lim }\frac{{\int }_0^{T}O(R_t)e^{\beta (V(s(R_t))-c(t))}dt}{{\int }_0^T{e}^{\beta (V(s(R_t))-c(t))}dt}$$

$$c(t)=-\frac{1}{\beta }\log \frac{\int ds{e}^{-\beta (F(s)+V(s,t))}}{\int ds{e}^{-\beta F(s)}}$$

$$d\tilde{t}=e^{\beta (V(s(R_t))-c(t))}dt$$

$$⟨O(R)⟩=\underset{T_{\tilde{t}}\to \infty }{\lim }\frac{1}{T_{\tilde{t}}}{\int }_0^{T_{\tilde{t}}}O(R_{\tilde{t}})d\tilde{t}$$

$${\psi }_i(R)=\sum _{k=1}^N{b}_{ik}{O}_k(R)$$

$$\mathbf{C}(\tau ){\mathbf{b}}_i=\mathbf{C}(0){\lambda }_i(\tau ){\mathbf{b}}_i$$

$$\text{C}(0)$$

$$\text{C}(\tau )$$

$$w(t)=\frac{e^{\beta (V(s(t))-c(t))}}{{\sum }_t{e}^{\beta (V(s(t))-c(t))}}$$

$$\tau =\sum _{t=t_0}^{{\tau }^{\prime }}{e}^{\beta (V(s(t))-c(t))}\Delta t$$

$$t_i^{\ast }=-\frac{\tau }{\log |{\lambda }_i|}$$

$$s_i(R)=\sum _{k=1}^{M<N}{b}_{ik}{O}_k(R)$$

Full text

Molecular Dynamics (MD) simulations have become pervasive in contemporary science, and are extensively used in fields as diverse as chemistry, biology and material science. Yet in spite of many successes, the limited time scale that can be explored in MD simulations severely limits their scope and power. In many cases the time scale problem results from the presence of metastable states separated by kinetic bottlenecks that render the transitions between such states so rare that they take place on a time scale that far exceeds what can be afforded by ordinary simulation methods.

For this reason, a plethora of enhanced simulation methods have been suggested. Starting from the pioneering work of Torrie and ValleauTorrie and Valleau (1977), a wide class of such methods rely on the identification of appropriate collective variables (CVs)Darve and Pohorille (2001); Huber, Torda, and Gunsteren (1994); Laio and Parrinello (2002); Wang and Landau (2001); Barducci (2008); Valsson, Tiwary, and Parrinello (2016). The CVs are a set of functions $s(R)$ of the atomic coordinates $R$ that describe those degrees of freedom whose sampling needs to be accelerated. This latter goal is achieved by adding to the interaction potential $U(R)$ an appropriately constructed bias $V(s(R))$, designed so as to accelerate sampling. Here we shall focus on one such method, namely well tempered metadynamics (WTMetaD)Barducci (2008); Valsson, Tiwary, and Parrinello (2016) that is enjoying increasing popularity and that offers the possibility in one of its variants, usually referred to as infrequent metadynamicsTiwary and Parrinello (2013), to calculate transition rates from metastable state to metastable state. However in WTMetaD the rate of convergence depends on an appropriate choice of CVs especially when it comes to calculating rates. In an important paper this has been pointed out by Tiwary and BerneTiwary and Berne (2016) who have proposed the spectral gap optimization of order parameters (SGOOP) method, that is based on maximum caliber and that has been shown to improve the quality of an initial CV guess in a spectacular way.

This work offers an alternative to SGOOP based on the signal processing technique Time-lagged Independent Component Analysis (TICA)Molgedey and Schuster (1994); Pérez-Hernández et al. (2013); Schwantes and Pande (2013). As in the case of Ref. 9 we start from a CV that is able to push the system from one metastable state to another albeit sluggishly and ameliorate it so as to greatly improve its efficiency. It has been shown in the important work of Noé and coworkersPérez-Hernández et al. (2013) that TICA provides an optimal solution of the variational approach to conformational dynamics (VAC), hence identifying slow order parameters which may serve as optimal CVs in a metadynamics simulation. We shall refer to this combination as the variational conformational dynamics approach to metadynamics (VAC-MetaD).

We recall first that in WTMetaD the bias potential $V(s(R),)$ is built on the fly by periodically adding a small repulsive Gaussian whose amplitude decreases as the simulation progresses. A remarkable feature of WTMetaD is that this stochastic, apparently out of equilibrium process, can be described by an ordinary differential equation whose asymptotic solution for the bias tends rigorously to the following equilibrium resultDama, Parrinello, and Voth (2014):

[TABLE]

where $F(s)$ is the free energy associated to the s CV given within an irrelevant constant by:

[TABLE]

where $\beta=1/k_{B}T$ is the inverse temperature and $\gamma$ is the so-called bias factor.

One of the consequences of the existence behind WTMetaD of an ordinary differential equation is that the reweighting of the trajectories can be done on the fly and the equilibrium expectation value of an operator $\langle O(R)\rangle$ can be calculatedTiwary and Parrinello (2014) as an average over the metadynamics run as:

[TABLE]

where $R_{t}$ are the atomic positions at time $t$ and the time dependent energy offset $c(t)$ is:

[TABLE]

a quantity that tends asymptotically to the reversible work performed on the system by the bias. If we introduce a new time scale $\tilde{t}$ such that

[TABLE]

we can recognise that $\langle O(R)\rangle$ can be written as an average over the $\tilde{t}$ time

[TABLE]

where $T_{\tilde{t}}=\int_{0}^{T}e^{\beta(V(s(R_{t}))-c(t))}dt$ is the total elapsed $\tilde{t}$ time. The times $t$ and $\tilde{t}$ measure the metadynamics and Boltzmann sampling progress respectively. It follows from Eqn. 6 and from the convergency properties of WTMetaD that we can think of WTMetaD as an ergodic dynamics in $\tilde{t}$ time that samples the Boltzmann distribution. We note however that the $\tilde{t}$ dynamics cannot be directly related to the unbiased dynamics but depends on the choice of CVs. Poor CVs lead to long convergence times, while good CVs lead to much shorter ones. In fact the very purpose of biasing the system is to turn rare events into frequent ones, and the time $\tilde{t}$ should not be confused with the actual unbiased time, but a measure of the extent of metadynamics enhancement of the sampling.

We want now to take advantage of progress made in the field of conformational dynamicsSchütte et al. (1999); Prinz et al. (2011) and the realization that a propagator associated to an ergodic dynamics can be spectrally decomposed into eigenvalues $\lambda_{i}(t)$ and eigenfunctions $\psi_{i}(R)$. The highest eigenvalue $\lambda_{0}$ is one since the system evolves towards the invariant distribution, while the others decay with timeSchütte et al. (1999). If the first sets of M non trivial eigenvalues are separated by a gap from all others, the corresponding eigenfunctions can be identified as the CVs that describe the slowest modes of the system. In order to practically benefit from these theoretical results one needs an approximation method that can deal with the fact that the propagator is very high dimensional. Luckily an approximate evaluation, based on a variational principle similar to the Raleigh-Ritz principle of quantum mechanics, has been suggestedNoé and Nüske (2013); Nüske et al. (2014). In this approach the dynamics is projected into a low dimensional space spanned by the basis functions $O_{k}(R)\quad k=1,...,N$ and the eigenfunctions are approximated by a linear expansion:

[TABLE]

The eigenvalues and eigenfunctions that best approximate the exact eigenvalues and eigenvectors are given by the solution of the following generalized eigenvalue equation:

[TABLE]

where $\tau$ is the lag time and ${\bf C}(\tau)$ is the matrix of the dynamical correlation functions $C_{j,k}(\tau)=\langle O_{j}(t)O_{k}(t+\tau)\rangle$,while $\lambda(\tau)_{i}$ are the eigenvalues and ${\bf b}_{i}$ are the expansion coefficients of eigenfunctions.

As discussed above, it is possible to map WTMetaD into a dynamics that asymptotically samples the Boltzmann distribution. Thus we make the ansatz that also for the $\tilde{t}$ dynamics the properties at the basis of the spectral decomposition described above hold at least asymptotically since metadynamics in this limit explores ergodically the Boltzmann distribution.

If this is so we can approximate the slow modes of the systems with the solutions of the generalized eigenvalue Eqn. 8 in which the time correlation functions are expressed as a function of the scaled time $\tilde{t}$ . While the variational approach Eqn. 8 aims at identifying the slowest CVs by varying ${\bf b}_{i}$ so as to maximize the eigenvalues $\lambda(\tau)$, metadynamics aims at generating a biased dynamics in which the slowest processes are fast. Hence, a successful application of VAC-MetaD should result in a biased dynamics whose slowest processes are fast, corresponding to small leading eigenvalues, $\lambda_{1},\lambda_{2},...$. Thus our strategy will be of choosing first a set of CVs expressed in the space spanned by an appropriate set of functions $O_{k}(R)$ and then perform a WTMetaD run to calculate the correlation function in $\tilde{t}$ time, plug them into Eqn. 8, and use the eigenvectors of the highest eigenvalues as improved CVs. If we perform a new WTMetaD simulation driven by such CVs then we should see a much more efficient sampling, and due to the acceleration of the slowest modes, we should observe a decrease in the relaxation time associated with the leading eigenvalues. We must add however that since metadynamics needs some incubation time $\tau_{c}$ before reaching the asymptotic limit in which Eqn. 3 holds, only after time $\tau_{c}$ will the trajectory yield a stable estimate of the eigenvalues. Only in this limit the eigenvectors of the slowest modes in Eqn. 8 will be used as new CVs.

We now explicate how the method works in practice. As stated above the VAC-MetaD procedure aims to find from an initial set of candidate order parameters $\textbf{O}=\{O_{k}(R)\}$ the optimal linear combination of these order parameters to use as a CV in metadynamics. Initially we give equal weight to each order parameter and perform a short biased simulation with the initial CV as $s^{(0)}(R)=\frac{1}{\sqrt{N}}\sum^{N}_{k}O_{k}(R)$. In order to approximate the unbiased slow modes from the biased metadynamics simulation, we need to reweight trajectory samples in our computation of the correlation matricesWu et al. (2017). The two correlation matrices needed to solve Eqn. 8 are

[TABLE]

with $\textbf{O}(t)$ the set of candidate order parameters at time $t$ and $w(t)$ the WTMetaD weight

[TABLE]

In practice the correlation functions of Eqn. 9 are symmetrized to ensure that the $\lambda$s are real valuedWu et al. (2017). Diagonalization of the correlation matrix $\textbf{C}(0)$ would give the usual reweighted principle components. In the time-lagged matrix $\textbf{C}(\tau)$ the lag time $\tau$ is given by the sum of the rescaled time steps

[TABLE]

Insertion of the matrices given by Eqn. 9 into Eqn. 8 and solving the generalized eigenvalue equation for the eigenfunctions ${\bf b}_{i}$ gives a new set of N transformed basis functions $s_{i}={\bf b}^{T}_{i}\textbf{O}$, which are an approximation for the eigenfunctions of the dynamical propagator (Eqn 7). Each of the N eigenvalues has an associated relaxation time given at a fixed chosen lag time $\tau$ by

[TABLE]

One expects to find a gap in timescales between the the $M<N$ slow modes corresponding to large $t^{*}$ which provides a natural way to select a subset of the N basis function, those corresponding to the slowest relaxation times (largest eigenvalues), as new CVs to be biased with WTMetaD. Thus the procedure yields the desired set of M CVs as

[TABLE]

with $\textbf{b}_{i}$ the expansion coefficients in Eqn. 8.

We shall test this procedure by studying in vacuum the conformational landscape of two simple peptides that have often been used as testing ground for new methods, alanine dipeptide (Ace$-$Ala$-$Nme) and alanine tetrapeptide (Ace$-$Ala3$-$Nme) which are shown in Fig. 1. All simulations were performed with the GROMACS5.0.5 packageHess et al. (2008) using the Amber99-SB force field with an integration time step of 2 fs. Trajectories were generated in the NVT ensemble with the temperature maintained at 300 K using the stochastic velocity rescaling thermostatBussi, Donadio, and Parrinello (2007). All metadynamics calculations were performed within the PLUMED2 pluginTribello et al. (2014) with Gaussian hills deposited every 500 integration steps and an initial hill height of 1.2 kJ/mol, a width of 0.03 units and a bias factor $\gamma$ of 15.

We start with the much studied case of alanine dipeptide. Here we depart from the usual procedure and besides the Rahmachandran angles $\phi$ and $\psi$, we consider additionally the angle $\theta$ that is known to be part of the reaction coordinateBolhuis, Dellago, and Chandler (2000) (see Fig. 1 a). Each angle is transformed according to $\Theta_{k}=0.5+0.5\cos(\Theta_{k}-\Theta_{0})$ with the reference angle $\Theta_{0}=1.2$ radTiwary and Berne (2016) and $\boldsymbol{\Theta}=\{\phi,\psi,\theta\}$. The initial CV is then taken to be a linear combination with equal weights of the transformed angles $\phi$, $\psi$, and $\theta$, i.e. $s_{0}=c_{1}\Theta_{\phi}+c_{2}\Theta_{\psi}+c_{3}\Theta_{\theta}$ with the $\{c_{1},c_{2},c_{3}\}$ initially all equal and normalized to 1. An upper and lower restraint on the $\theta$ angle was introduced at $\pm 0.5$ radians. Following the usual TICA procedure we subtract the reweighted mean given by Eqn. 3 from the raw trajectory data and work in a zero-mean vector space. From an initial WTMetaD trajectory of 20 ns biasing $s_{0}$ we compute the two correlation matrices in Eqn. 9 with $\textbf{O}=\boldsymbol{\Theta}$. The trajectory of the initial CV $s_{0}$ is shown in Fig. 2 (top left). It can be seen that such a CV is able to induce transitions between the different conformers however such transitions are not very frequent. This is reflected by the fact that the highest $\lambda(\tau)$ decays more slowly than the other eigenvalues reflecting a difficulty of the chosen CV to promote transitions. The middle row of Fig. 2 shows the relaxation times $t^{*}_{i}$ associated with each eigenvalue given by Eqn. 12 at a lag time $\tau=800$ ps chosen within the regime for which the eigenvector coefficients are stable. Fig. 2 middle left clearly shows a slow process with a dominant large eigenvalue and a clear separation of time scales. The eigenvector associated with this highest eigenvalue, computed in the asymptotic regime where the eigenvectors are constant with respect to the lag time $\tau$, thus obtained is used as a new CV $s_{1}=\textbf{b}_{1}^{T}\boldsymbol{\Theta}$. If we use as CV this eigenvector associated with the highest $\lambda$ the exploration of the conformational space is greatly accelerated as depicted in Fig. 2 top right, and the relaxation rates of decay of the different $\lambda$s become comparable, reflecting the fact that the rare event has been made no longer rare (Fig. 2 middle right). It is interesting to note that in the optimized CV much of the weight is of the angle $\phi$ but both $\psi$ and $\theta$ are part of the optimal CV.

In Fig. 2 (bottom) we also show the free energies associated to the initial (left) and final (right) CV. The conformational landscape of alanine dipeptide is well known and characterized by a deeper basin in which conformers C5 and C7eq can easily interconvert and while conformer C7ax is higher in energy separated by a sizeable barrier from the first two. It can be seen that although the FES is fully converged, these physical pictures are hidden in the free energy representation provided by the initial CV but clearly evident when using the final CV.

We now turn to the second example for which we have applied our approach, that of alanine tetrapeptide (Ala3) shown in Fig. 1 b), a case already considered inValsson and Parrinello (2014); Tiwary and Berne (2016). The simulation is started as in Ref.  9 by taking as collective coordinate a linear combination with equal weights of the set of six dihedral angles $\boldsymbol{\Theta}=\{\phi_{1},\psi_{1},\phi_{2},\psi_{2},\phi_{3},\psi_{3}\}$ transformed as before so that $\Theta_{k}=0.5+0.5\cos(\Theta_{k}-\Theta_{0})$. As can be seen in Fig. 3 (left) the exploration of the conformational space is somewhat inefficient with rare conformational transitions. We project the $\tilde{t}$ dynamics in the space of the chosen angles and observe that the decay times of the two topmost eigenvalues are significantly slower than that of the others. This suggests to project the free energy surface in the space of the two topmost eigenvectors. As shown in Fig. 4, in this representation the seven different minima identify conformers that have a different arrangement of the three dihedral angles $\phi_{1}$, $\phi_{2}$, and $\phi_{3}$. The conformer in which all the three $\phi$s have all positive values is seldom visited and in this representation are projected into the top left tail of the central minimum.

An improvement in sampling efficiency is obtained if one uses the topmost eigenvector as a new CV (see the central panels in Fig. 3), but still it can be seen that a relatively slow process remains. However the behavior in time of the $\lambda$s cries out for the use of at least two CVs. If this is done, (see Fig. 3 far right panels) the improvement is amazing and it offsets the cost of using two CVs instead of one. In fact within a 20 ns run we get better converged results than the extensive parallel tempering run in Ref  24 that used an aggregated time equivalent to $8\times 500$ ns. These lead to a well converged and smooth free energy surface (see Fig. 4 right).

It is difficult at this stage to compare the relative practical merits of the SGOOP method to ours. Extensive test on a variety of applications will be needed. As far as we can tell, for the two cases examined here they seem to have comparable performances in the case of Ala3 when we use only the topmost eigenvector as collective variable. It is likely that the two methods will be complementary. However a difference in philosophy must be underlined. In SGOOP one reweighs a one-dimensional projection of the FES to find the optimal linear combinations of CVs. Here we are reweighting the simulation time to take advantage of the variational formalism of conformation dynamics whose solution provides an optimal estimate for the true dynamical propagator.

In summary, there is a growing interest in applying dimensionality reduction techniques on a larger candidate set of possible CVs to find a subset of generic good CVs that can be used for enhanced sampling. A widely used example is principle component analysis (PCA) which projects the data along the direction of largest variance. On the other hand, it is well known in the field of conformational dynamics that time-lagged independent component analysis on high dimensional coordinates from molecular dynamics can be useful for constructing Markov state models. In this letter we have taken this insight and combined it with well tempered metadynamics to find a set of optimal CVs for further enhanced sampling. We surmise finally that the method can be adapted to other sampling methods by appropriately changing the weights in Eqn. 9