Free energy and kinetics of conformational transitions from Voronoi tessellated milestoning with restraining potentials
Milestoning is a method aimed at reconstructing the statistical properties of the long-time dynamics of a system by exploiting the crossing statistics of a set of hypersurfaces, called the \xe2\x80\x9cmilestones\xe2\x80\x9d, placed along the reaction coordinate [Faradjian and Elber, J. Chem. Phys. 2004, 120, 10880]. Recently, Vanden-Eijnden and Venturoli [J. Chem. Phys. 2009, 130, 194101] showed that when a complete Voronoi tessellation of the configurational space is available, milestoning can be reformulated to utilise the statistics from a series of independent simulations, each confined within a single cell via strict reflections at the boundaries. As a byproduct, this \xe2\x80\x9cVoronoi tessellated milestoning\xe2\x80\x9d method also permits to compute the free energy of the tessellation. Here, the method is extended to support the usage of differentiable restraining potentials to confine the trajectories within each cell.
1 Introduction<!>2.1 Soft walls restraining potentials<!>2.2 Free energy of the tessellation<!>2.3 Milestoning<!>3 Results and discussion<!>3.1 The Mueller potential<!>3.2 Solvated Alanine Dipeptide<!>4 Concluding remarks<!>
<p>Conformational changes in large biomolecules are complex and slow processes taking place on very long timescales, which often extend well beyond the reach of brute force molecular dynamics (MD) simulations. An emerging class of techniques attacks this problem by first trying to determine an optimal reaction pathway (or pathways) for the transition in a space of very high dimensionality, without making any a priori assumptions about the mechanism.1–5 Once such an optimal pathway has been identified, a strategy then consists in inferring the statistical properties of the long-time dynamics by exploiting the information harvested from relatively short independent trajectories.6–9 Of particular interest is the milestoning method introduced by Faradjian and Elber,8 in which the dynamical properties of a reactive event are reconstructed out of a series of short trajectories between a set of hypersurfaces (the "milestones"), disposed sequentially along the reaction pathway between the reactant and product.10–14 Recently, Vanden-Eijnden and Venturoli extended milestoning by identifying the edges of the cells of a complete Voronoi tessellation of configurational space as the milestones.9 It was shown that the "Voronoi tessellated milestoning" could be formulated to require only a series of independent simulations, each confined within a single Voronoi cell, leading to a simplification and increase in robustness of the original algorithm. The confinement was realized by a simple collision rule, via strict hard-wall reflections (velocity inversion) at the boundaries.4</p><p>Despite its formal simplicity, the hard-wall boundary condition involves modifications of the dynamical propagation algorithm at the heart of MD source codes. This makes its implementation in widely used biomolecular simulation MD packages (e.g., AMBER,15 CHARMM,16 GROMACS, 17 and NAMD18) somewhat cumbersome, and may even affect performance. With the aim of extending the range of applicability of the method, in this work we propose and test an alternative strategy, based on the introduction of continuous and differentiable restraining potentials to confine the system within a given cell. We illustrate how the formalism presented in Ref. 9 can still be used, with minor modifications, to compute the equilibrium probability and rates of transitions for the original dynamics. The possibility to use potentials will facilitate the application of the Voronoi tessellated milestoning because the introduction of user-defined external forces is straightfoward. We note that the confinement with potentials can also affect performance, because portions of trajectories that are transiently out of the cells are discarded in the analysis, but this effect can be minimized by proper tuning of the parameters in the potentials.</p><p>In the remaining of this letter, we first describe the details of the restraining potentials and the key quantities to be computed for free energy and milestoning calculations, and then we illustrate the method by applying it to two simple examples, a numerical toy model and solvated alanine dipeptide.</p><!><p>Let us consider a molecular system with coordinates x and potential energy U(x). In most biomolecular applications, the dimensionality of the system is very high, but most of the variables can be unimportant for the description of a reactive mechanism. For this reason, it is customary to introduce a smaller number of collective variables to characterize the process. Let us indicate the set of collective variables, functions of x, by {z(x)}. Suppose we are given a set of K points in z-space, (z1, z2, …, zK), that we call centroids. These support a complete and unique Voronoi tessellation of the original Cartesian space of the system, where each cell α is defined such that ||z(x) − zα|| < ||z(x) − zγ|| for all γ ≠ α; || · || indicates the norm in some metric, which shall be assumed Euclidean for the sake of simplicity.</p><p>It was shown previously,9 via a reformulation of the milestoning algorithm,8 that the statistical properties of the long-time dynamics of the system can be reconstructed from independent simulations, each confined inside one of the K cells. The key feature required for this result is that the confined dynamics must be equivalent to that from a long unbiased trajectory passing through the same set of cells.4,9 More specifically, the confinement must leave unperturbed the dynamical properties of the systems when it is in the interior of the cell, as well as the probability flux in and out of the cells. In Ref. 9, the confinement was realized with a strict reflection rule,4 i.e. by reversing the velocity of the atoms when a trajectory attempted to cross the hyperplane between two neighboring cells α and γ. An alternative strategy, explored here, is to introduce planar half-pseudo-harmonic restraining potentials to confine the system within each cell,</p><p>
(1)uα(x)=∑γ=1Kkw2(n^αγ·[z(x)−z¯αγ])2H(n^αγ·[z(x)−z¯αγ]), where the unit vector n̂αγ = (zγ − zα)/||zγ − zα|| is normal to the hyperplane separating the cells α and γ, z̄αγ = (zγ +zα)/2 is the mid-point between the centroids α and γ, and H is the Heaviside step-function. The real parameter kw determines the strength of the restraining potential. Hereafter, we will refer to the present approach as the soft walls restraints, to emphasize the distinction with the strict hard walls reflections used in Refs. 9 and 4.</p><p>An illustration of uα(x) is given in Figure 1. The restraining potential uα(x) is zero in the interior of the cell α, and it acts as a penalty function towards the cell boundary when the trajectory crosses any of the hyperplanes separating α from other cells.</p><p>It can be understood from the following argument that the fundamental requirements of the confined dynamics can be met when the restraints of Eq. (1) are used. First, since uα(x) = 0 inside the cell, using Eq. (1) guarantees that the equilibrium distribution inside the cell is the same as that of an unbiased simulation, apart for small errors at the boundaries related to time discretization. Moreover, this implies that we also have the correct fluxes in and out of the cell (otherwise, their effect would propagate in and spoil the distribution inside), and hence the pieces of trajectory of the restricted simulation are indistinguishable from pieces of an unbiased trajectory passing through the same cell.</p><p>It is possible to compute the equilibrium probability and rates of transitions for the original dynamics from the simulations confined within the Voronoi cells by the hard walls reflections.9 In the following sections, we show that, with minor modifications, the formalism can still be used with the soft walls restraints.</p><!><p>The conservation of probability flux through the boundaries of the cell gives a way to compute πα, the equilibrium probability of the system to be in cell α, and the associated free energy Gα = −β−1 logπα (where β = 1/kBT, with kB the Boltzmann constant and T the temperature), as summarized hereafter.4,9 The rate of escape from cell α to cell γ, conditional on the system being in the cell α, is defined as ναγ = Nαγ/Tα, where Nαγ is the number of collisions with the boundary separating the cells α and γ, and Tα is the total simulation time spent inside cell α. The normalized equilibrium probability πα solves the following equations involving the rates of escape ναγ: 
(2)∑γ=1γ≠αKπγνγα=∑γ=1γ≠αKπαναγ,∑α=1Kπα=1.</p><p>This equation expresses that, at statistical steady state, the total probability flux in cell α must be equal to the flux out of α by conservation of the total probability. In Ref. 9, Tα was simply the total simulation time with the system confined in cell α by means of hard walls. With soft walls, Eq. (2) is still valid, as long as one counts only the portion of trajectory spent inside the cell α as the time Tα.</p><!><p>Identifying the edges of the Voronoi tessellation as milestones, the dynamics of the system is reduced to that of a discrete state continuous-time Markov chain in the state space of the milestones indices.9 By indexing the milestones as i and j, this amounts to defining a rate matrix qi j, whose elements are given by</p><p>The factors Ni j and Ri in this equation can be expressed9 in terms of average properties extracted from simulations confined to the cell α, weighted by the the equilibrium probability πα of finding the system in cell α,</p><p>Here Tα is the duration of the simulation confined in cell α, 
NijÎ± is the number of transitions from edge i to edge j observed during this simulation, and 
RiÎ± is the total time that edge i was the last edge to be hit during this simulation (i.e. 
∑iRiα=Tα). The only requirement to use Eqs. (4) with soft walls is to prune the trajectory to its portion that is strictly inside the cell α, before computing 
NijÎ±,RiÎ± and TÎ± (i.e. the parts spent outside the cell in the soft walls region must be discarded).</p><p>The rate matrix qi j specifies completely the dynamics of the Markov process, and hence it can be used to compute many important quantities as for example the mean first passage times (MFPTs) from any milestone to any other.9,14 For instance, if 
τiN with i = 1, …, N − 1 denote the MFPTs from milestone i to milestone N (
τNN=0 by definition), then these MFPTs can be computed by solving the linear system of equations: 
(5)∑j=1N−1qijτjN=−1,i=1,…,N−1.</p><!><p>In the following, we illustrate the implementation of milestoning with soft walls restraints on two simple systems, the Mueller potential19 and the solvated alanine dipeptide. Prior to this, an important point concerning the position of the centroids supporting the Voronoi tessellation deserves a special attention. While Eqs. (3)–(5) can be applied to any set of hypersurfaces, it has been shown that the formalism gives exact MFPTs if the hypersurfaces used as milestones are chosen as the isocommittor surfaces of the reaction.14 From a practical view point, the isocommittor surfaces can be calculated approximately by the string method and its variants.3–5,20–22 In the string method, the transition path is represented as by an ordered sequence of K discrete "images", {z(1), z(2), …, z(K)}, in the space of collective variables z(x). This suggests9 to first use the string method to determine the transition path by optimizing the position of the K images, and then use these images as centroids to support the Voronoi tessellation. Naturally, the edges between the cells are hyperplanes perpendicular to the optimized path, and hence, approximations of the isocommittor surfaces.</p><!><p>We first consider a simple 2-dimensional system evolving on the Mueller potential energy surface according to Langevin dynamics (i.e. we take {z(x,y)} ≡ (x,y), the coordinates of the system). The same model was studied in Ref. 9, so for comparison we took the same centroids as generators of the Voronoi tessellation. These are 18 equidistant images (the red circles in Figure 2), computed with the finite-temperature string method.4 All simulations were performed by integrating Langevin equations of motion with the second order algorithm of Ref. 23, by taking β−1 = 20, friction 400, and a timestep Δt = 10−4. Hard walls and soft walls trajectories of 108 steps were generated in each cell.</p><p>Figure 2 shows the Mueller potential together with the Voronoi tessellation associated to the images from the string. The edges of the cells were taken as milestones. The successive positions of the system at every timestep along a piece of simulation with soft walls are also shown as white dots. Note that, as already pointed out before, the portion of trajectory inside the cell samples the equilibrium distribution associated to the underlying potential in the cell.</p><p>Figure 3 shows the free energy Gα = −β−1 logπα, where πα is the solution of Eq. (2), corresponding to the tessellation shown in Figure 2. Results are presented for soft walls (SW) and hard walls (HW) simulations, and numerical integration (Exact). The numbering of the cells goes from 1 to 18 from top-left to bottom-right in Figure 2. Table 1 summarizes results for MFPTs calculations considering the transition from the black dashed line to the black continuous line in Figure 2. The table also shows the mean and extremal values of the elements needed for the calculation of the qi j matrix defined in Eq. (3), computed from soft and hard walls simulations, and from a long free trajectory (2·109 steps). The excellent agreement shows that the confinements do not introduce biases in the calculations of qi j.</p><!><p>In order to test the implementation of soft (and hard) walls in a realistic biomolecular context, we examine the transition from helical to extended conformation of the alanine dipeptide solvated in explicit water. We describe the transition here by using two collective variables, the φ and ψ backbone dihedral angles, neglecting the role of the solvent degrees of freedom.24 We used the CHARMM16 code for all calculations, with CHARMM22 all atom force field25 and the TIP3P model26 for water molecules.</p><p>The free energy landscape in φ and ψ variables was computed with the single-sweep method,27 and it is shown for reference in Figure 4 (energies are in kcal/mol). A minimum free energy path in (φ, ψ) space was computed using the string method.3 Twenty images along this path are represented as circles in Figure 4, and they were used as centroids for the Voronoi tessellation. Note that in our calculation we focus on the reaction channel with φ < 0, where the main metastable states are. In principle, we could also investigate the kinetics of transitions between the states with φ > 0, but these are much less often visited, and hence their existence does not affect much the rates we compute. The edges of the Voronoi cells are represented as red lines in Figure 4 (taking into account periodicity). The dashed and continuous black lines are the start and ending milestones for MFPTs calculation. The hard walls confinement condition was implemented into the Nosé-Hoover (NH) dynamics subroutine of CHARMM.28 The successive values of dihedral angles at every timestep along pieces of simulations with hard walls are shown as white and gray dots. NH dynamics was also used for soft walls simulations. In this case, the forces on atoms coming from the potential (1) were implemented and added to the standard CHARMM forces.</p><p>Figure 5 shows the free energy corresponding to the tessellation shown in Figure 4. Results are for hard walls (HW) and soft walls (SW) simulations (with two different values of the penalty constant kw). In all cases, the trajectories in each cell were 1 ns long. Numbering of the cells goes from 1 to 20 from top to bottom in Figure 4. In order to assess convergence of the calculation, we monitored during the simulations the values of the escape rates from the cells, ναγ. Figure 6 shows the rates from cell 2 versus the length of the entire confined trajectory, for soft walls simulations with kw = 200 kcal/mol/rads2. Escape rates from other cells converge on the same time-scale.</p><p>MFPTs were computed considering the transition from the black dashed line to the black continuous line in Figure 4. Results are summarized in Table 2, together with the mean and extremal values of the elements needed for the calculation of the qi j matrix defined in Eq. (3). The excellent agreement shows again that the confinements do not introduce biases in the calculations of qi j.</p><p>The MFPT for the transition from dashed to continuous milestone in Figure 4 estimated from milestoning with hard and soft walls simulations is about 28 picoseconds. For comparison, the MFPT between the same two milestones was computed by direct counting of consecutive hittings, along a free, unbiased simulation of 80 ns, obtaining 29.97 ps. These values agree with results from other CHARMM calculations of solvated alanine dipeptide in similar simulation setup, where the MFPT for the same transition was estimated to be of about 30 ps.29,30</p><!><p>Milestoning with Voronoi tessellation9 is a method to reconstruct the dynamical properties of complex reactive systems by matching together informations obtained from multiple trajectories, each confined in a different cell of a Voronoi tessellation of configurational space. In this letter, we have illustrated how the formalism can still be applied, with minor modifications, when the confinement in the cells is realized via continuous and differentiable potentials. With respect to the original formulation, where the confinement is realized with strict reflections (velocity inversion at the boundary), the possibility to use potentials will facilitate the application of the method, because user-defined external forces are easy to introduce in most MD codes, and do not require modification of the dynamical propagators. However, the confinement with potentials can also affect performance, because portions of trajectories that are transiently out of the cells are discarded in the analysis. Note also that we did not optimize for efficiency here, which would require adjusting the penalty constant kw in Eq.(1).</p><p>As a final remark, let us point out that the method proposed in this letter, like the one in Ref. 9, can be trivially distributed on multiple processing nodes, since it is based on independent simulations with no required communication among them.</p><!><p>An illustration of the potential uα(x) used to constraint a trajectory in cell α, Eq. (1). Portions of three cells α, β, and γ with centroids (zα, zβ, zγ) are shown. The potential is zero in cell α (white region), and its units are arbitrary. Cell edges are represented as black lines. Note that the edge between cells β and γ (dashed black line) has no effect on the restraint, and it is shown only for the purpose of clarity.</p><p>Mueller potential with a set of milestones corresponding to the edges of Voronoi cells generated by images along the converged string from finite temperature string method. The successive positions of the system at every timestep along a piece of simulation with soft walls are also shown as white dots. MFPTs discussed in the text are calculated from the dashed to the continuous black line.</p><p>Free energy Gα = −β−1 log πα of the Voronoi tessellation from hard (HW) and soft (SW) walls simulations on the Mueller potential, compared with the exact one obtained by numerical integration. Numbering of the cells goes from top-left to bottom-right in Figure 2.</p><p>Backbone dihedral free energy surface of solvated Alanine dipeptide (units are kcal/mol), with a set of milestones (red lines) corresponding to the edges of Voronoi cells generated by images along the converged minimum free energy path (red circles), corresponding to the helical to extended transition. The successive values of dihedral angles at every timestep along pieces of simulations with hard walls are shown as white and gray dots. MFPTs discussed in the text are calculated from the dashed to the continuous black line.</p><p>Free energy of the Voronoi tessellation from hard and soft walls (HW and SW) sampling for the solvated alanine dipeptide. Numbering of the cells goes from 1 to 20 counting from top to bottom in Figure 4.</p><p>Rates of escape from cell 2 (i.e. associated to the second string image counting from top to bottom in Figure 4), ν21, top line, and ν23, bottom line, as a function of trajectory length for the alanine dipeptide soft walls simulation with kw = 200 kcal/mol/rads2.</p><p>Kinetics for the Mueller potential for the transition from the black dashed milestone to the black continuous milestone in Figure 2. The agreement of all quantities shows that the dynamical properties of the trajectory inside the cells are left unperturbed by the two different confinement methods, i.e. it is as if the confinement was absent. Errors on MFPTs were obtained from five different runs.</p><p>Kinetics for the solvated Alanine dipeptide for the transition from the black dashed milestone to the black continuous milestone in Figure 4.</p>