Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2892875

Formats

Article sections

- Abstract
- 1. Introduction
- 2. Directional Milestoning – theory
- 3. Applications of Directional Milestoning
- 4. Discussions and conclusions
- REFERENCES

Authors

Related links

J Chem Theory Comput. Author manuscript; available in PMC 2011 January 1.

Published in final edited form as:

J Chem Theory Comput. 2010; 6(6): 1805–1817.

doi: 10.1021/ct100114jPMCID: PMC2892875

NIHMSID: NIHMS198913

See other articles in PMC that cite the published article.

Milestoning is a method for calculating kinetics and thermodynamics of long time processes typically not accessible for straightforward Molecular Dynamics (MD) simulation. In the Milestoning approach, the system of interest is partitioned into cells by dividing hypersurfaces (Milestones) and transitions are computed between nearby hypersurfaces. Kinetics and thermodynamics are derived from the statistics of these transitions. The original Milestoning work concentrated on systems in which a one-dimensional reaction coordinate or an order parameter could be identified. In many biomolecular processes the reaction proceeds via multiple channels or following more than a single order parameter. A description based on a one-dimensional reaction coordinate may be insufficient. In the present paper we introduce a variation that overcomes this limitation.

Following the ideas of Vanden-Eijnden and Venturoli on Voronoi cells that avoid the use of an order parameter (J. Chem. Phys. 2009, 130, 194101), we describe another way to “Milestone” systems without a reaction coordinate. We examine the assumptions of the Milestoning calculations of mean first passage times (MFPT) and describe strategies to weaken these assumptions. The method described in this paper, Directional Milestoning, arranges hypersurfaces in higher dimensions that “tag” trajectories such that efficient calculations can be done and at the same time the assumptions required for exact calculations of MFPTs are satisfied approximately.

In the original Milestoning papers trajectories are initiated from an equilibrium set of conformations. Here a more accurate distribution, that mimics the first hitting point distribution, is used. We demonstrate the usage of Directional Milestoning in conformational transitions of alanine dipeptide (in vacuum and in aqueous solution) and compare the correctness, efficiency, and statistical stability of the method with exact MD and with a related method.

Milestoning is a method to calculate kinetics and thermodynamics of molecular systems that evolve on long time scales typically not accessible for straightforward Molecular Dynamics (MD) simulation.^{1}^{–}^{8}

Straightforward Molecular dynamics can be used to compute rate of reactions. In these applications coordinates and velocities are initiated in the reactant state and the equations of motion are integrated until the product state is reached. While considerably promising, there are caveats: (i) the numerical integration of a typical biomolecular process is computationally demanding and may not be feasible; (ii) actual realizations of reactive trajectories are noisy, making their analysis difficult and may require significant filtering to recover useful signals.

In Milestoning, the conformational space between the reactant and the product is partitioned by a set of dividing hypersurfaces called Milestones (Fig. 1). An ensemble of initial conditions is prepared at each Milestone and trajectories are simulated from each initial point until another nearby Milestone is reached. These trajectories are significantly shorter and trivially parallelized compared to a reactive trajectory of the overall process. The efficiency of the algorithm is discussed in ^{1}.

A schematic arangement of Milestones (dashed lines) in a two-well potential. Also shown is a trajectory (dotted line) starting on a second Milestone and terminating on the first one.

In the original milestoning papers,^{1}^{,}^{3} a theory that relates the statistical properties of the short trajectories initiated on each Milestone and the overall rate was developed. In the present work we consider a variant of the Markovian limit of Milestoning,^{1}^{,}^{2} a method that uses only the first moments of local first passage time (LFPT) distributions. The advantage of the Markovian limit of Milestoning is that it is easier to implement and is statistically more stable. As we will show in Section 2.1 it calculates the overall mean first passage times (MFPT) accurately, given that certain assumptions are met. Milestoning in its complete settings (non-Markovian) provides a useful alternative if more detailed understanding of the reaction process is desired, for example if the reaction is non-exponential in time.

Vanden-Eijnden et. al.,^{4} considered reaction dynamics with overdamped Langevin. It was shown that if Milestones are chosen as isocommittor surfaces, i.e. surfaces for which the probability of reaching the product state before the reactant is constant, then Milestoning calculation of the MFPT using Brownian dynamics is exact. However, determination of exact isocommittor surfaces can be very difficult in practice.

Other limits in which Milestoning is expected to be accurate are available for systems near equilibrium. As outlined in the original Milestoning papers,^{1}^{,}^{3}, even when other surfaces are used (surfaces that are not isocommittors) Milestoning can still work well. If successive crossing events of Milestones are sufficiently separated in time to “lose” velocity memory Milestoning was illustrated to provide accurate results. This assumption is achieved in practice by placing Milestones sufficiently far from each other such that the average termination time of trajectories is at least a few hundred femtoseconds.^{1}

In Section 2 we propose a variant of Milestoning in the Markovian limit which we call Directional Milestoning (DiM) – the dividing hypersurfaces are redefined in more than one dimension to capture features of the reaction (e.g. multiple reaction channels or multiple collective variables) that at the same time maintain the concept of Milestone separation, e.g. trajectories initiated on any Milestone have time to “lose memory” before terminating on other Milestones.

The original Milestoning approach approximates the initial ensemble on each hypersurface by equilibrium distribution. To be exact the initial distribution at a Milestone must be the first hitting point distribution (FHPD). A first hitting point is a phase space point on the Milestone crossed for the first time by a trajectory arriving from a nearby hypersurface. The distribution of these phase space points is complex and a closed form of it is known only for overdamped Langevin dynamics in low dimensions.^{4}

In recent work,^{7} Vanden-Eijnden and Venturoli proposed a modification of Milestoning that avoids a generation of initial ensembles on each of the dividing surfaces. As we discuss later their approach is more accurate compared to the original Milestoning for the generation of the FHPD. Memory loss, however, is harder to control in the new approach. To improve the accuracy of the original Milestoning approach while retaining some of its advantages we propose in Section 2.4 another way to approximate FHPD, which is better than the original Milestoning.

In Section 3 we illustrate the Directional Milestoning (DiM) for the calculation of MFPT of a conformational transition of alanine dipeptide, both in vacuum and in water. We compare Directional Milestoning with exact Molecular Dynamics and with the related method Markovian Milestoning with Voronoi Tessellation (MMVT).^{7} We illustrate that as the complexity of the underlying energy surface increases, DiM becomes more effective. Discussions and conclusions are in presented Section 4.

We discuss below an extension of Milestoning that avoids the use of a reaction coordinate. Instead of placing hypersurfaces orthogonal to a one-dimensional curve as introduced in the original papers,^{1}^{,}^{3} we define the interfaces (Milestones) based on a set of coordinates (images) that sample the conformational space of the biophysical process under consideration. (Two of the images define the reactant and the product state.) These images may be obtained from long time simulations, high temperature trajectories, replica exchange simulations, etc., as discussed later in examples in the article. Having *N* images *X*_{1},…, *X _{N}* placed in the conformational space, we intuitively want to arrange Milestones as interfaces between the images, which is the approach taken in the Voronoi Tessellation of Markovian Milestoning.

$${M}_{i\to j}\equiv \left\{X|d{(X,{X}_{i})}^{2}=d{(X,{X}_{j})}^{2}+{\mathrm{\Delta}}_{i}^{2}\phantom{\rule{thinmathspace}{0ex}}\text{and}\phantom{\rule{thinmathspace}{0ex}}\forall k\phantom{\rule{thinmathspace}{0ex}}d(X,{X}_{j})\le d(X,{X}_{k})\right\},$$

(1)

where *d*(*X*, *Y*) is a distance function of images *X* and *Y* and Δ_{i} = min_{j ≠ i}*d*(*X _{i}*,

Example of Milestones according to definition (1). Conformational images are represented as black dots, Milestones related to the central image are displayed as dashed lines. A trajectory coming to the central region (gray, dotted) terminates on one of **...**

The proper selection of the conformational images *X*_{1},…, *X _{N}* will be explained in more detail in Section 2.3; for now we assume their arbitrary placement. If Δ

It is important to emphasize that the proposed placements of Milestones is not a tessellation. In accord with the definition of the original Milestoning, a trajectory is identified by the last Milestone that it passes and not by its actual current position. A memory is carried out in time until the trajectory crosses another interface (Milestone).

Trajectories from *X _{i}* to

In the rest of the manuscript we will use Roman subscripts do denote image index (as was done in the previous section) and Greek letters to denote Milestones. Consider the mean first passage time (MFPT) from any Milestone α to a given target Milestone β. We define it as follows: a trajectory is *assigned* to a Milestone α if the last Milestone it has passed through is α. *One-step transition* from a Milestone α to a Milestone β(β ≠ α) is a change of assignment of a trajectory from α to β. This step is clearly on a coarse Milestoning level and does not mean a single Molecular Dynamics step, which we will call a time-step. If such an event is possible we say that α connects to β. Note that by the definition given in equations (1), if α connects to β, the second index of α (e.g. *M*_{i → j}) must be equal to the first index of β (*M*_{j → k}). The first hitting point distribution on β, ρ_{β}(*p*), is the distribution of phase space points (denoted by *p*) at which an equilibrium trajectory passes through β numerous times while the previous Milestone it passes through was not β. In further discussion only the relative weight of trajectories that pass through β is important so we can choose to normalize ρ_{β}(*p*) such that ∫ ρ_{β}(*p*) *dp* = 1. We denote by τ_{αβ}(*p*) the mean time of all trajectories that start from the phase space point *p* in α and terminate on Milestone β (possibly crossing other Milestones on the way). Integrating the last entity over *p*, weighting it by the probability that *p* is a phase space point at which an equilibrium trajectory hits β for the first time, ∫ τ_{αβ}(*p*) ρ_{α}(*p*) *dp* τ_{αβ}, we obtain the MFPT from α to β.

Let the distribution of one-step transitions from α to β be *T*_{αβ}(*p*, *q*, *t*), where *p* is the phase space point at which a trajectory starts in α and *q* is the phase space point at which the trajectory changes its assignment to β after time *t*. *T*_{αβ}(*p*, *q*, *t*) is normalized in such a way that if we integrate over *t* and *q* we get conditional probability of a trajectory reaching β in one step given that it originates from *p* in α : ∫∫*T*_{αβ}(*p*, *q*, *t*) *dqdt* = *P*(β | α, *p*), or alternatively ∑_{β}∫∫*T*_{αβ}(*p*, *q*, *t*) *dqdt* = 1. Note that by the definition of trajectory assignment, *T*_{αα}(*p*, *q*, *t*) = 0 for all *p* and *q* (since a trajectory cannot change its assignment from α to α).

Assuming that the phase space point *p*(*t* + *dt*) can be determined from *p*(*t*) only, as is true for most microscopic dynamics (e.g. Newtonian, or Langevin dynamics, but not Generalized Langevin dynamics) we make the following argument: The MFPT from α to β, τ_{αβ}, is defined as the weighted average of termination times of trajectories from α to β. Each trajectory, starting at *p* in α jumps in one step to some other Milestone γ(γ ≠ α) at phase point *q* and then in multiple steps (possibly 0, if γ = β) continues to β. Consider all the trajectories that jump in one-step from *p* in α to *q* in γ exactly in time *t* and then eventually reach β (in potentially different total time). Since the microscopic dynamics is Markovian we can replace the contribution of these trajectories to τ_{αβ} by (*t* + τ_{γβ}(*q*)) weighted by sum of the weights of all of them (which is ρ_{α}(*p*)*T*_{αγ}(*p*, *q*, *t*)). By doing this for all possible combinations γ of and *q* we get the following equation:

$$\begin{array}{cc}\langle {\tau}_{\alpha \beta}\rangle \hfill & ={\displaystyle \sum _{\gamma}{\displaystyle \int {\displaystyle \int {\displaystyle \int {\rho}_{\alpha}(p){T}_{\alpha \gamma}}(p,q,t)\left(t+\langle {\tau}_{\gamma \beta}(q)\rangle \right)\mathit{\text{dp dq dt}}}}}\hfill \\ \hfill & ={\displaystyle \sum _{\gamma}{\displaystyle \int {\rho}_{\alpha}(p)\left({\displaystyle \int {\displaystyle \int {T}_{\alpha \gamma}}}(p,q,t)\phantom{\rule{thinmathspace}{0ex}}t\phantom{\rule{thinmathspace}{0ex}}\mathit{\text{dq dt}}\right)}\mathit{\text{dp}}}\hfill \\ \hfill & +{\displaystyle \sum _{\gamma}{\displaystyle \int \langle {\tau}_{\gamma \beta}(q)\rangle \phantom{\rule{thinmathspace}{0ex}}\left({\displaystyle \int {\displaystyle \int {\rho}_{\alpha}(p){T}_{\alpha \gamma}(p,q,t)\phantom{\rule{thinmathspace}{0ex}}\mathit{\text{dp dt}}}}\right)\phantom{\rule{thinmathspace}{0ex}}\mathit{\text{dq}}}}\hfill \end{array}$$

(2)

The first term of the above equation can be reduced as

$$\begin{array}{c}{\displaystyle \sum _{\gamma}{\displaystyle \int {\rho}_{\alpha}(p)\left({\displaystyle \int {\displaystyle \int {T}_{\alpha \gamma}(p,q,t)\phantom{\rule{thinmathspace}{0ex}}\mathit{\text{t dq dt}}}}\right)\mathit{\text{dp}}}}\hfill \\ ={\displaystyle \sum _{\gamma}{\displaystyle \int {\rho}_{\alpha}(p)\left(\frac{{\displaystyle \int {\displaystyle \int {T}_{\alpha \gamma}(p,q,t)\phantom{\rule{thinmathspace}{0ex}}\mathit{\text{t dq dt}}}}}{{\displaystyle \int {\displaystyle \int {T}_{\alpha \gamma}(p,q,t)\mathit{\text{dq dt}}}}}{\displaystyle \int {\displaystyle \int {T}_{\alpha \gamma}(p,q,t)\mathit{\text{dq dt}}}}\right)\phantom{\rule{thinmathspace}{0ex}}\mathit{\text{dp}}}}\hfill \\ ={\displaystyle \sum _{\gamma}{\displaystyle \int {\rho}_{\alpha}(p)\left(\langle {t}_{\alpha \gamma}(p)\rangle \phantom{\rule{thinmathspace}{0ex}}P(\gamma |\alpha ,p)\right)\mathit{\text{dp}}}}\hfill \\ ={\displaystyle \int {\rho}_{\alpha}(p)\left({\displaystyle {\sum}_{\gamma}\langle {t}_{\alpha \gamma}(p)\rangle \phantom{\rule{thinmathspace}{0ex}}P(\gamma |\alpha ,p)}\right)\mathit{\text{dp}}}\hfill \\ ={\displaystyle \int {\rho}_{\alpha}(p)\langle {t}_{\alpha}(p)\rangle}\phantom{\rule{thinmathspace}{0ex}}\mathit{\text{dp}}=\langle {t}_{\alpha}\rangle ,\hfill \end{array}$$

(3)

where *t*_{αγ}(*p*),*t*_{α}(*p*), and *t*_{α} are average times of one-step transitions from *p* α to γ, from *p* α to any other Milestone, and from α to any other Milestone (averaged over *p*), respectively. In the second term of equation (2) the average time from *q* γ to β is weighed by a factor depending on the phase space point *p* α ! To overcome this problem we use the following assumption: *The distribution at which any Milestone* γ *is hit does not depend on the Milestone to which the trajectory was assigned before the hit*:

$$\forall \alpha ,\gamma :\text{}{\rho}_{\gamma}(q)\propto {\displaystyle \int {\rho}_{\alpha}(p){T}_{\alpha \gamma}(p,q,t)}\mathit{\text{dp dt}}.$$

(4)

It is easier to illustrate the properties of equation (4) if we consider a one-dimensional arrangement of Milestones in which the forward and the backward Milestones occupy the same spatial coordinates. Consider a Milestone α that is pointing forward and is therefore denoted for the clarity of this discussion by α +. There are two Milestones that initiate trajectories that may terminate at α +. They are (α − 1)+ and (α − 1)−. Hence they occupy the same place in space but have their velocities pointing in the opposite directions. The assumption of equation (4) states that it does not matter if we start at (α − 1)+ or at (α − 1)−, both Milestones will generate the same hitting point distribution on α+. If the initial direction of the velocity de-correlates quickly there should be no difference in the results from Milestone (α − 1)+ and (α − 1)−. In this case the assumption formulated in equation (4) will be satisfied. Indeed, we observed empirically in reference ^{9} that even the usual Milestoning works well when the velocity de-correlates. This empirical formulation is now formulated mathematically. In higher dimension we will also require spatial de-correlation.

The multiplicative factor in the above equation is determined by the fact that if both sides of equation (4) are integrated over *q*, the left side equals to 1 and the right side to *P*(γ|α); the conditional probability that if a trajectory changes its assignment from α it changes to γ. Therefore using the above assumption the second term of equation (2) reduces to ∑_{γ}*P*(γ|α) τ_{αβ} and we obtain the final form for the MFPT:

$$\langle {\tau}_{\alpha \beta}\rangle =\langle {t}_{\alpha}\rangle +{\displaystyle \sum _{\gamma}P(\gamma |\alpha )\langle {\tau}_{\gamma \beta}\rangle .}$$

(5)

The set of equation (5) is extended by boundary conditions τ_{ββ} = 0, τ_{β} = 0, and ∀γ *P*(γ|β) = 0. It is a set of linear equations for all the τ_{αβ} that can be solved by any standard linear solver. The size of the problem (the number of Milestones) never exceeded a few hundreds in our hands. Equation (5) can be directly generalized for considering more than a single target Milestone (e.g. all incoming interfaces to the folded state of a peptide). Alternative equations equivalent to equation (5) were derived in ^{1}^{,}^{4}. These equations are independent of the type of microscopic dynamics that we use (e.g. overdamped Langevin or Newtonian). The system of linear equations (5) relates the overall rate (τ ‘s) with the local kinetics information (*t*_{α} and *P*(γ|α)). Milestoning collects this local information in a more effective way than running an ensemble of trajectories from α to β. On each Milestone α, *N*_{α} phase space points are sampled from the FHPD ρ_{α} (see Section 2.4 for details). As a second step, each of the sampled phase space points is propagated in time until a connected Milestone is reached. The termination times of these trajectories are typically several orders of magnitude shorter than the overall MFPT of the system. Furthermore the trajectories between Milestones are independent of each other and thus can be run in parallel. For each Milestone γ connected to α we record *N*_{αγ}- the number of trajectories that are initiated on α and terminated on γ. We also record α, the mean termination time of all *N*_{α} trajectories regardless of their terminal Milestone. The collected information {*N*_{α},α} is used to estimate the required entities for equation (5) as

$$\begin{array}{ccc}P(\gamma |\alpha )\cong {N}_{\alpha \gamma}/{N}_{\alpha}\hfill & \text{and}\hfill & \langle {t}_{\alpha}\rangle \cong {\overline{t}}_{\alpha}\hfill \end{array}.$$

(6)

In practice instead of using equation (6) we employ Bayesian inference on the collected data to calculate the MFPT supported by the data as well as an estimate of the statistical error due to the finite size of collected data. This procedure is described in detail in Appendix B.

The use of equation (5) for calculating MFPT depends on validity of the assumption expressed in equation (4). It has been shown in ^{4} that the assumption formulated in equation (4) holds if overdamped Langevin dynamics is used and the Milestones are chosen as isocommittor surfaces. To our knowledge there is no efficient algorithm that identifies exact isocommittor surfaces and scales moderately with system size. However, there are other ways of satisfying equation (4). Instead we base our strategy on selecting Milestones according to equation (1), making sure that Milestones are sufficiently separated to allow for a memory loss of trajectories as outlined in the arguments of reference.^{1} Consider a pair of connected Milestones *M*_{i → j}, *M*_{j → k} (defined by coordinate images *X _{i}*,

The first step of Milestoning is to sample the initial conditions on each Milestone α from the first hitting point distribution ρ_{α}(*p*). An analytical expression for ρ_{α}(*p*) is in general unknown. In ^{4} the authors provided the formula ρ_{α}(*x*) *e*^{−βV(x)}|*q*(*x*)| for the case of overdamped Langevin dynamics with Milestones being placed as isosurfaces of the committor function *q*(*x*). The last formula includes the gradient of committor function *q*(*x*) which is difficult to get in high dimensions.

Instead of computing ρ_{α}(*p*) exactly (no exact expression is available for Newtonian dynamic), we approximate it. First, phase space points are sampled from the equilibrium distribution at Milestone *M*_{i → j}. It can be done either by running an MD simulation constrained to the Milestone,^{1}^{,}^{3} or by employing the Umbrella Sampling technique (see appendix C and ^{10}). The second step involves filtering each of the sampled phase points to determine those that are indeed first hitting events of *M*_{i → j}. Exact verification tracks each of the sampled phase space points *p* back in time and tests termination on one of the incoming Milestones to the cell *X _{i}*(

To retain high efficiency we track the trajectory back in time only until it reaches an empirical test boundary that is placed at a distance *d* on the *X _{i}*’s side of the target Milestone

Illustration of sampling of the first hitting point distribution of trajectories initiated on the lower gray Milestone and terminating on the top (target) Milestone. The FHPD on the target Milestone (blue) is centered in the left basin, which is different **...**

In principle we can follow the trajectory back in time until one of the incoming Milestones to *X _{i}* (

To demonstrate an application of Directional Milestoning we compute MFPT of the transition between α helix and β sheet conformations in solvated alanine dipeptide (Fig. 4). The thermodynamics and kinetics of alanine dipeptide has been investigated in several studies ^{1}^{,}^{8}^{,}^{14}^{–}^{16}. In aqueous solution two dihedral angles, ϕ and ψ, shown in Fig. 4, are adequate coarse variables for the dynamics of the peptide. We therefore use a 2-norm distance in the reduced space of ϕ and ψ as the distance metric in the definition of Milestones (periodicity of the angles was taken into account in the calculation of a distance between two torsion angles).

The new module for Directional Milestoning was created in the program MOIL,^{17} and is available at https://wiki.ices.utexas.edu/clsb/wiki. The peptide molecule is solvated in a periodic box (20 Å)^{3} of 248 TIP3P water molecules. The OPLS force field,^{18} is used with electrostatics real space cutoff of 9 Å augmented with Particle Mesh Ewald summation. Van der Waals interactions are cut at a distance of 8 Å. All calculations were run in NVT ensemble at a temperature of 303 K by employing a weak Andersen thermostat that acts only on the center-of-mass motion of the water molecules.^{19} The probability of velocity re-sampling was set to 5 · 10^{−4} per fs. For a water box of this size an average of 13 water molecules had their velocities re-sampled in a 100 fs interval. This weak coupling does not change the transition rate obtained from NVE (Newtonian) simulations (with initial conditions sampled from the NVT ensemble). The free energy surface as a function of the two dihedral angles (ϕ, ψ) is shown in Fig. 5. It was calculated from statistics of a 340 ns long MD simulation. The white region of the map was not visited by the trajectory. There are two local free energy minima corresponding to an α helix conformation (ϕ, ψ = −100, −40) and to a β sheet conformation (ϕ, ψ = −100, −140).

Free energy profile of alanine dipeptide as a function of the two dihedral angles ϕ and ψ. It was calculated from statistics of a 340 ns long MD simulation. Images for DiM calculations are placed at the positions of the red numbers and **...**

The height of the free energy barrier between the two metastable regions at 303K is less than 2*k _{B}T* and the transitions between the metastable states are rapid on the trajectory time scale so the MFPTs can be estimated from straightforward MD simulations directly. We have performed five independent MD simulations of 68 ns. In each of the simulations more than 1000 transitions between the metastable regions occurred. The MFPT of α → β transition is 66.4 ps (± 2.7 ps) and that of the opposite transition is 53.8 ps (± 4.6 ps). We set up the Milestoning calculation by placing six images in the conformational space in the positions ϕ

Note that the employed dynamics is almost deterministic and thus a trajectory reflected from an interface (procedure required in MMVT) would approximately track itself back in time. Therefore we have slightly modified the MMVT protocol in a way suggested by Vanden-Eijnden in a private communication: instead of reversing the velocities of all the degrees of freedom at a cell interface, only the velocities of peptide atoms are reversed. This modification should not influence the statistics of observed fluxes through the interfaces since only the peptide degrees of freedom are used in the definition of cell boundaries.

Both methods, DiM and MMVT, perform well in this scenario, though MMVT is more efficient for this simple system. If enough sampling is done, both techniques provide reasonable estimates of MFPTs between the metastable regions, the systematic error is lower for MMVT (6 % and 10%) as compared to our method (10 % and 18 %). Analysis of MMVT on the same system was performed recently.^{8} A different force field was employed in [8] and the MFPT reported differs by a factor of two from our calculations; however the relative error of MMVT for the reported α → β transition is about 6%, which is comparable to our result. Results of β → α transition were not reported in ^{8}. Table 1 shows that MMVT needs about 2–3 times less CPU time compared to DiM to converge. DiM requires more computations in these setting since each interface of MMVT is effectively doubled for the two different directions. Furthermore, additional computation is needed in DiM to sample initial phase space points on each interface. In this one-dimensional set-up of Milestones with relatively large separation between Milestones and low free energy barrier MMVT is more efficient and as accurate as DiM. However, we will show below that with smaller separation between the interfaces, multi-dimensional arrangement of milestones, and rougher energy landscapes, DiM is better.

Even though previous Milestoning studies calculated accurately MFPTs on alanine dipeptide, memory effects in the system are not negligible. First hitting point distributions (in terms of ϕ angles) for the Milestones *M*_{4→5} and *M*_{6→5} are shown on Fig. 6. There is a noticeable difference between distributions of first hitting points on the Milestone *M*_{4→5} and on the Milestone *M*_{6→5}. As shown on the figure, the approximate sampling described in Section 2.4 distinguishes the first hitting point distributions arriving from different directions to the region of image *X*_{5} reasonably well.

Distributions of ϕ angle of the first hitting point conformations of the region of image *X*_{5} (located at ψ = 80°): distributions observed in a long MD simulation for conformations arriving to the hypersurface at *X*_{5} from the hypersurface **...**

In Table 2, we examine the use of directional Milestones on this system. The table shows that transitions between the six Milestones (if direction is not part of the description) are not Markovian. If no memory effects were present in the system then the probability of transiting to Milestone *i*+1 from Milestone *i* would not depend on the Milestone visited before *i*, i.e. the second and the forth columns of Table 2 would be the same within the error bars. We however see differences of up to 21% (for *i*=5) or by a factor of up to 2.2 (for *i*=1). One can see that the values of these relative probabilities estimated by Directional Milestoning (columns 3 and 5 in Table 2) are in good agreement with the true values.

In the second experiment we examine both methods (DiM and MMVT) on the same system with Milestones in more than one dimension. This experiment is performed to empirically illustrate that placing Milestones in a non-linear arrangement does not compromise accuracy of DiM calculations. Images are placed in a two-dimensional grid covering the accessible space at the target temperature (conformations with torsional angle ϕ < 0). For DiM, 18 images are placed in the positions marked on 1, …, 18 on Fig. 7a). Each image has 8 incoming Milestones and 8 outgoing Milestones (displayed in solid and dashed on Fig. 7a) respectively). We calculated the MFPT from *M*_{12→11} (or *M*_{10→11}) to the union of *M*_{10→9} and *M*_{8→9} for the β → α transition. The MFPTs from these two Milestones differ from each other by about 0.3 ps and we report their average in Table 3. The opposite transition (α → β) was defined in the equivalent way.

Placement of images on a two dimensional grid. (a) DiM settings: total of 18 images, located at position of numbers in the plot, are placed in a two dimensional grid. For two of the images, *X*_{11} and *X*_{14}, the outgoing (dashed) and incoming (solid) Milestones **...**

For MMVT the images were placed in slightly different positions than for DiM (see Fig. 7b) such that the Milestones inferred by the Voronoi Tessellation are in equivalent positions to those used in Directional Milestoning. For the α → β transitions, we calculated the MFPT of trajectories starting from the two white Milestones in Fig. 7b (*M*_{2↔3} and *M*_{8↔9}) and terminating at the union of the red Milestones. MFPT of the transitions from these two starting points differ by less than 0.2 ps so only their average is reported in Table 3. The β → α calculation was performed in the equivalent way (from the two central Milestones in the β sheet conformation (ψ = 140°) to the union of all the Milestones with ψ = −40°).

The results of both methods are listed in Table 3. The accuracy of Directional Milestoning is not compromised by multidimensionality; hence DiM works well for higher dimensions or higher connectivity of Milestones. The relative error of the MMVT method increased to 33 % (31 %). We think that this is mainly due to the corners between Milestones in the MMVT arrangement that cause rapid termination times between nearby Milestones and unwanted correlations between touching Milestones. Evidence of this can be seen in Fig. 8.

In vacuum there are two stable conformers *C*_{7eq} and *C*_{ax} of alanine dipeptide (Fig. 9). The state *C*_{7eq} is further split into two sub-states denoted by *C*_{7eq} and ${C}_{7\phantom{\rule{thinmathspace}{0ex}}\text{eq}}^{\prime}$ (located at *X*_{26} in Fig. 9b) separated by a small barrier. We calculate the MFPT of transition from *C*_{7eq} to *C*_{ax} at two different temperatures, 400 K and 350 K, using Langevin dynamics. This is performed by calculating MFPT starting from each of the incoming Milestones to *C*_{7eq} region (green on Fig. 9) and considering union of the incoming Milestones to the region *C*_{ax} (red on Fig. 9) as the final state. The MFPT is not sensitive to exact identity of the starting Milestone (variation of less than 2%) therefore an average MFPT from all green Milestones is considered. The friction constant of Langevin dynamics was set to 30 ps^{−1}.

The shown landscape is an adiabatic ϕ, ψ energy map. The energy is minimized while constraining the ϕ and ψ dihedrals to specified values. Placement of (a) 24 images, (b) 63 images in the conformational space based on the **...**

The images were generated by the following expansion. We start with the set of images *S* = {*X*_{1}, *X*_{2}}, where *X*_{1} is a conformation located at *C*_{ax} and *X*_{2} at *C*_{7eq}. Then we iteratively pick an image *X* from the set *S* and “expand” it: We launch trajectories starting from *X* with randomly initiated velocities and run each of these trajectories until it departs at least a pre-specified distance δ from *X*. Then we cluster the set of end points of these trajectories to existing images in *S* and potentially add new images to the set *S* if there are end points that are farther than δ from all images of *S*. We repeat this process until no new images are generated, i.e. we have tried launching trajectories from all images in *S* and all end coordinates are in *S*. There are three parameters in this algorithm: (i) the distance cutoff δ, (ii) the number of expanding trajectories *N _{e}*, and (iii) the clustering algorithm employed. For alanine dipeptide we have used expectation-maximization as a clustering algorithm,

By using different values for δ we obtained sets of images of size 24 (for δ_{1}) and 63 (for δ_{2}); both are shown on Fig. 9. The tessellations shown in black in this figure are only approximate since they are based on the Euclidean distance in (ϕ, ψ) space, where the real interfaces (Milestones) are defined using the RMSD distance. The MFPT of the transitions between the metastable conformations are significantly longer than those in the solvated peptide due to higher free energy barriers. Table 4 and Table 5 summarize the results of the Milestoning calculations in this system. At the high temperature (400 K) both methods, DiM and MMVT, predict accurate MFPT from *C*_{7eq} to *C*_{ax} (with systematic error of about 10%). MMVT needs to run about 1.5 µs MD simulations to obtain converged results, while DiM requires about 2.5 µs. Both of them provide significant speed up against straightforward MD simulation, even though a rough estimate of MFPT of the *C*_{7eq} to *C*_{ax} transition can be obtained by running about 11 independent MD simulations (equivalent to 4 µs of the total simulation time); however, both MMVT and DiM can be trivially parallelized to thousands of CPUs, shortening the actual time to perform the calculation.

When the temperature is lowered to 350 K (see Table 5) the *C*_{7eq} to *C*_{ax} transition is slower with MFPT of about 2.0 µs. As listed in Table 5, Directional Milestoning calculates the MFPT with systematic error of about 15% with as few as 7.5 µs of total simulation time. That is significant speedup compared to straightforward MD since DiM can be easily parallelized on thousands of processors. MMVT fails to calculate the MFPT accurately. The main reason for this failure is poor statistics. An important difference between DiM and MMVT is that DiM allocates computational resources to each Milestone, where MMVT allocates the computational resources to a cell. If a transition between two specific interfaces in a cell is needed to describe the reaction and the transition is significantly less likely than transitions between other interfaces of the cell, then sampling this transition using MMVT is inefficient. A simple realization of this effect is the existence of a barrier in the middle of the cell. In that case MMVT trajectory is likely to be confined to a one minimum, to collide with the same interface many times (hits that do not count for the statistics) and to record only a few successful transitions to the other minimum. In contrast DiM launches a large number of short trajectories. These trajectories terminate quickly, and contribute to the statistics.

In DiM, sampling is done (extensively) at the interfaces, so the probability of observing a transition between interfaces of interest is greatly enhanced, since at least one end of the transitional event is sampled extensively. A potential problem in DiM is a large number of interfaces that may make sampling expensive. To avoid sampling irrelevant interfaces (at a given temperature) trajectories are initiated at few initial interfaces and only interfaces that are hit at least once during the DiM calculation are sampled and launched. We stop the DiM calculation when the process converges (i.e. no new interfaces besides those already sampled are reached). For the MMVT calculation with 24 images, many cells cover a relatively large part of the conformational space with a rough energy landscape (see for example cell *X*_{6} on Fig. 9a). This arrangement may cause poor statistics for those regions since the trajectories spend most of their time in low free energy regions, rarely visiting interfaces higher in free energy. To increase the probability of having a double hit at the two desired surfaces, we run the same calculation with 63 images as well. But even when 63 images are used, the allocation of computational resources is highly unbalanced. For example, we consider the frequency of hitting the interfaces 49 → 47 and 33 → 47 (displayed in white on Fig. 9b) that are both important for the overall MFPT. In both, 49 and 33 cells, confined simulations of total time of 2.25 µs hit a cell boundary more than 2 · 10^{7} times. However, the interface 33 → 47 is hit only 17 times and the interface 49 → 47 only 7 times. In contrast DiM allocates an equal number of starting trajectories to each of the Milestones and transitions from Milestones located near the transition states are sampled as well as other Milestones. We have not experimented with any selection criterions for allocation of computational resources to different cells (in MMVT) or to different Milestones (in DiM) but both methods may benefit from selective allocation of resources to “important regions” of conformational space.

In this paper, we proposed a method to compute dynamics in high dimensions called Directional Milestoning. We have shown that the mean first passage times between Milestones can be calculated accurately given that the distribution at which a Milestone is hit does not depend on the previously assigned Milestone (the assumption formulated in equation (4)). Directional Milestoning arranges dividing hypersufaces in a special way, aiming to satisfy the above assumption: i) Milestones in DiM are made directional, so the local progress of the reaction (going from the region of *X _{i}* to

The algorithm, while based on the trajectory fragments of Milestoning, is a step in the direction of Transition Interface Sampling ^{11}^{–}^{13} (TIS) and the Forward Flux Sampling (FFS) methods^{22}^{,}^{23} compared to the original Milestoning. Here we use some trajectory tracking. The main difference between these methods and Directional Milestoning is that TIS and FFS are tracking trajectories all the way back to the reactant state. This tracking has the advantage of not relying on any assumption about the initial ensemble on an interface like is done in Milestoning. On the other hand, sampling of trajectories in TIS and FFS is computationally more expensive than in Milestoning because every attempted trajectory in these methods is tracked back to the reactant state where in (Directional) Milestoning a trajectory is tracked only until it reaches a different Milestone. Computations of trajectory fragments can be done in Milestoning in a massively parallel way. The PPTIS method uses a conceptually similar approach of trajectory fragments ^{24}.

An important distinction of Directional Milestoning compared to TIS, FFS, and the original Milestoning is that it allows for arbitrary arrangement of Milestones in conformational space, not necessarily following a linear arrangement along an order parameter or a reaction coordinate. A similar (arbitrary) arrangement of interfaces is used in the MMVT method,^{7} nonequilibrium umbrella sampling method,^{25}^{,}^{26} and Trajectory Parallelization and Tilting method.^{27} The last two techniques are using short trajectories in cells and balance the fluxes between cells. Recently the non-equilibrium umbrella sampling^{26} was illustrated to be more efficient than FFS ^{28} The Weighted Ensemble approach was also shown recently to work without a reaction coordinate ^{29}.

We have compared DiM with MMVT and showed that the performance of MMVT (in terms of effectiveness and correctness) is comparable to that of DiM in some of the examples, but that the correctness and/or effectiveness of MMVT can be compromised in systems with high free energy barriers, or in cells with two interfaces that are hard to reach. Another problem for straightforward implementation of MMVT is the existence of corners between Milestones along more than one dimension that contribute to termination times that are too short. So while DiM is in general somewhat slower than MMVT it provides reliable results more consistently, including cases in which MMVT fails.

We also would like to comment on the similarities (and the differences) of our approach to the Markov State Model (MSM -- for a recent study see ^{30}). In the applications of MSM that we are aware of, long to very long Molecular Dynamics trajectories at normal conditions are used to estimate transition times and population of different cells. MMVT and DiM are designed to avoid such long trajectories (at the cost of approximate matching of probability densities at the interfaces). Once a sample of conformational space is available (which can be done in numerous ways, reaction path calculations, replica exchange simulations, or high temperature trajectories) only very short Molecular Dynamics trajectories are required to estimate the local kinetics. These short trajectories can be trivially parallelized providing profound computational saving compared to straightforward Molecular Dynamics simulations. While significant progress has been made in parallelizing a single trajectory ^{31}, overhead still remains and special hardware that is frequently used is more expensive to buy and to maintain

This research was supported in part by NIH grant GM059796 and NSF grant 0833162 to RE.

Let *X _{i}* and

**Proof of Lemma A.1**: From def. (1) of *M*_{i → j}:

$$d{({X}_{i},A)}^{2}-d{({X}_{j},A)}^{2}={\mathrm{\Delta}}_{i}^{2}$$

By using the Pythagoras theorem for triangles and *X _{i}AB* and

$$\begin{array}{cc}d{({X}_{i},B)}^{2}-d{({X}_{j},B)}^{2}\hfill & =\left(d{({X}_{i},A)}^{2}+d{(A,B)}^{2}\right)-\left(d{({X}_{j},A)}^{2}+d{(A,B)}^{2}\right)\hfill \\ \hfill & =d{({X}_{i},A)}^{2}-d{({X}_{j},A)}^{2}={\mathrm{\Delta}}_{i}^{2}\hfill \end{array}$$

Consequence of Lemma 1: *M*_{i → j} is a hyperplane segment perpendicular to *X _{i}X_{j}*.

Let *S _{ij}* be the hyperplane perpendicular to the line segment

**Proof of Lemma 2**: Since both *S _{ij}* and

$$\begin{array}{c}{(d({X}_{i},{P}_{\mathit{\text{ij}}})+d({S}_{\mathit{\text{ij}}},{M}_{i\to j}))}^{2}-{(d({X}_{j},{P}_{\mathit{\text{ij}}})-d({S}_{\mathit{\text{ij}}},{M}_{i\to j}))}^{2}={\mathrm{\Delta}}_{i}^{2}\hfill \\ d({S}_{\mathit{\text{ij}}},{M}_{i\to j})=\frac{{\mathrm{\Delta}}_{i}^{2}}{4d({X}_{i},{P}_{\mathit{\text{ij}}})}=\frac{{\mathrm{\Delta}}_{i}^{2}}{2d({X}_{i},{X}_{j})}\hfill \end{array}$$

We describe an estimate of the statistical error of a milestoning calculation from a single set of collected data using Bayesian reasoning. As shown in Section 2, equation (5), repeated here as (B.1),

$$\langle {\tau}_{\alpha \beta}\rangle =\langle {t}_{\alpha}\rangle +{\displaystyle \sum _{\gamma}P(\gamma |\alpha )\langle {\tau}_{\gamma \beta}\rangle}$$

(B.1)

relates MFPTs (τ_{αβ}) and local kinetics entities (*t*_{α} and *P*(γ|α)). Milestoning aims to estimate *t*_{α} and *P*(γ|α) by launching *N*_{α} trajectories from each Milestone α. *N*_{αγ} of them terminate on the Milestone γ and the mean incubation time (time to termination) of all trajectories is _{α}. In Bayesian inference a statistical model of the transitions among Milestones is needed. We closely follow and extend notation used in the analysis of Markovian Milestoning with Voronoi Tesselations; for more details consult ^{7}. The same kinetic formulas (with different notation) are also available from ^{9}. We assume a continuous Markov jump process between the Milestones controlled by a transition matrix *Q* defined in the following way: let the probability distribution of the system over all the Milestones be, **ρ** = (ρ_{1},…,ρ_{N}) where ρ_{α} is the probability that the system is assigned to a Milestone α. Under continuous Markov jump process, **ρ** behaves as:

$$\dot{\mathbf{\rho}}=\mathbf{\rho}Q.$$

(B.2)

For transition matrix *Q*, by definition *q*_{αα} = −∑_{α ≠ β}*q*_{αβ} and it can be shown by simple algebra that *P*(β|α) = *q*_{βα}/∑_{γ ≠ α}*q*_{γα} and *t*_{α} =1/∑_{γ ≠ α} *q*_{γα} (for derivation see for example ^{1}^{,}^{2}^{,}^{7}). By plugging the last three identities to the linear system (B.1) it reduces to

$$Q\text{'}\langle \tau \rangle =-1,$$

(B.3)

where **τ** is the row vector (**τ _{1β}**, …,

$$L(\{{N}_{\alpha \gamma},{\overline{t}}_{\alpha}\}|Q)={\displaystyle \prod _{\alpha}{\displaystyle \prod _{\gamma \ne \alpha}{q}_{\alpha \gamma}^{{N}_{\alpha \gamma}}{e}^{-{q}_{\alpha \gamma}{N}_{\alpha}{\overline{t}}_{\alpha}}}}.$$

(B.4)

By using the Bayes’ rule the likelihood that the true transition matrix is *Q* given the collected data, *L*(*Q*|{*N*_{αγ}, _{α}}), is:

$$L(Q|\{{N}_{\alpha \gamma},{\overline{t}}_{\alpha}\})\propto {\displaystyle \prod _{\alpha}{\displaystyle \prod _{\gamma \ne \alpha}{q}_{\alpha \gamma}^{{N}_{\alpha \gamma}}{e}^{-{q}_{\alpha \gamma}{N}_{\alpha}{\overline{t}}_{\alpha}}P(Q)}},$$

(B.5)

where *P(Q)* is the prior probability distribution of *Q* without seeing any data (typically this is set to uniform if we do not have any prior knowledge about the system). Eq. (B.5) is typically used in maximum likelihood estimators, e.g. one estimates unknown entity *Q* with *Q**, the matrix that maximizes likelihood *L*(*Q*|{*N*_{αγ}, _{α}}). In this particular case, *Q** has form ${q}_{\alpha \gamma}^{*}={N}_{\alpha \gamma}/[{N}_{\alpha}{\overline{t}}_{\alpha}]$ what is in agreement with estimators given in Eq. (6) in the main text. Instead of using purely *Q** for calculations of MFPTs we can examine whole distribution of transition matrices according to Eq. (B.5) and understand what is the distribution of MFPTs consistent with the data collected. Therefore we typically sample a number of (typically 300) transition matrices from distribution (B.5) and look at the variance of MFPTs predicted by them. If standard deviation of MFPTs is large it suggests that more data about the system shall be collected. We report standard deviation obtained by this algorithm in the results of Section 3.

As described in Section 2.4 the equilibrium ensemble from a Milestone *M*_{i → j} is used to sample the first hitting point distribution on the Milestone *M*_{i → j}. The Milestone *M*_{i → j} is defined in equation (1) as ${M}_{i\to j}\equiv \{X|d{(X,{X}_{i})}^{2}=d{(X,{X}_{j})}^{2}+{\mathrm{\Delta}}_{i}^{2}\phantom{\rule{thinmathspace}{0ex}}\text{and}\phantom{\rule{thinmathspace}{0ex}}\forall k\phantom{\rule{thinmathspace}{0ex}}d(X,{X}_{j})\le d(X,{X}_{k})\}$, where {*X*_{1},…,*X _{K}*} is a set of images in the conformational space. In practice we work with the following approximation of

$$\begin{array}{c}{d}_{\mathit{\text{ij}}}(X)\equiv d{(X,{X}_{j})}^{2}-d{(X,{X}_{i})}^{2}+{\mathrm{\Delta}}_{i}^{2}\hfill \\ {M}_{i\to j}^{\text{'}}\equiv \{X|\forall k,d(X,{X}_{k})\ge d(X,{X}_{j})\phantom{\rule{thinmathspace}{0ex}}\Lambda -\lambda \le {d}_{\mathit{\text{ij}}}(X)\le 0\}\hfill \end{array}$$

(1.1)

Clearly as λ→0, ${M}_{i\to j}^{\text{'}}$ converges to *M*_{i → j}. We have used λ = 0.5° or λ = 0.01 Å for the calculations on alanine dipeptide.

To sample conformations in ${M}_{i\to j}^{\text{'}}$ from equilibrium distribution the following Umbrella Sampling protocol is employed. We run NVT trajectory of the system (using Andersen thermostat) with a modified potential function *U* and examine a conformation every few steps (every 100 – 400 fs for examples described in this paper). If an examined conformation belongs to ${M}_{i\to j}^{\text{'}}$ it is saved; otherwise it is discarded. If conformation is saved, corresponding velocities are sampled from Boltzmann distribution. The potential function *U* is modified to bias the system towards the region ${M}_{i\to j}^{\text{'}}$ in the following way:

$$\begin{array}{c}U\text{'}(X)=U(X)+{U}_{\mathit{\text{ij}}}^{1}(X)+{U}_{\mathit{\text{ij}}}^{2}(X)\hfill \\ {U}_{\mathit{\text{ij}}}^{1}(X)=\{\begin{array}{cc}{K}_{1}{d}_{\mathit{\text{ij}}}{(X)}^{2}\hfill & \text{if}\phantom{\rule{thinmathspace}{0ex}}{d}_{\mathit{\text{ij}}}(X)>0\hfill \\ {K}_{1}{\left({d}_{\mathit{\text{ij}}}(X)-\lambda \right)}^{2}\hfill & \text{if}\phantom{\rule{thinmathspace}{0ex}}{d}_{\mathit{\text{ij}}}(X)<-\lambda \hfill \\ \hfill 0\hfill & \text{otherwise}\hfill \end{array}\hfill \\ {U}_{\mathit{\text{ij}}}^{1}(X)=\{\begin{array}{cc}{K}_{2}{\left(d(X,{X}_{k})-d(X,{X}_{j})\right)}^{2}\hfill & \text{if}\phantom{\rule{thinmathspace}{0ex}}d(X,{X}_{k})<d(X,{X}_{j})\hfill \\ \hfill 0\hfill & \text{otherwise}\hfill \end{array}\hfill \end{array}$$

By definition for $X\in {M}_{i\to j}^{\text{'}}$, *U'*(*X*) = *U*(*X*) and therefore saved points from ${M}_{i\to j}^{\text{'}}$ are sampled with the true equilibrium probabilities. If on the other hand NVT trajectory of the system is outside of the region ${M}_{i\to j}^{\text{'}}$, the terms ${U}_{\mathit{\text{ij}}}^{1}$ and/or ${U}_{\mathit{\text{ij}}}^{2}$ force the system to return back to ${M}_{i\to j}^{\text{'}}$, the strength of this bias is controlled by force constants *K*_{1} and *K*_{2} (both are set to 10^{3} Kcal mol^{−1} rad^{−2} or 10^{4} Kcal mol^{−1} Å^{−2} for alanine dipeptide system).

1. West AMA, Elber R, Shalloway D. The Journal of Chemical Physics. 2007;126:145104. [PubMed]

2. Shalloway D, Faradjian AK. The Journal of Chemical Physics. 2006;124:054112. [PubMed]

3. Faradjian AK, Elber R. The Journal of Chemical Physics. 2004;120:10880. [PubMed]

4. Vanden-Eijnden E, Venturoli M, Ciccotti G, Elber R. The Journal of Chemical Physics. 2008;129:174102. [PMC free article] [PubMed]

5. Elber R. 2007;92:L85. [PubMed]

6. Kuczera K, Jas GS, Elber R. The Journal of Physical Chemistry A. 2009;113:7461. [PMC free article] [PubMed]

7. Vanden-Eijnden E, Venturoli M. The Journal of Chemical Physics. 2009;130:194101. [PubMed]

8. Maragliano L, Vanden-Eijnden E, Roux B. Journal of Chemical Theory and Computation. 2009;5:2589. [PMC free article] [PubMed]

9. West AMA, Elber R, Shalloway D. Journal of Chemical Physics. 2007;126

10. Torrie GM, Valleau JP. Journal of Computational Physics. 1977;23:187.

11. van Erp TS, Bolhuis PG. Journal of Computational Physics. 2005;205:157.

12. Moroni D, van Erp TS, Bolhuis PG. Physica A: Statistical Mechanics and its Applications. 2004;340:395.

13. Moroni D, Bolhuis PG, van Erp TS. The Journal of Chemical Physics. 2004;120:4055. [PubMed]

14. Ren W, Vanden-Eijnden E, Maragakis P, E W. The Journal of Chemical Physics. 2005;123:134109. [PubMed]

15. Maragliano L, Vanden-Eijnden E. The Journal of Chemical Physics. 2008;128:184110. [PubMed]

16. Ensing B, De Vivo M, Liu Z, Moore P, Klein ML. Accounts of Chemical Research. 2005;39:73. [PubMed]

17. Elber R, Roitberg A, Simmerling C, Goldstein R, Li H, Verkhivker G, Keasar C, Zhang J, Ulitsky A. Computer Physics Communications. 1995;91:159.

18. Jorgensen WL, Tirado-Rives J. Journal of the American Chemical Society. 2002;110:1657.

19. Juraszek J, Bolhuis PG. 2008;95:4246. [PubMed]

20. Hartley H. Biometrics. 1958;14:174.

21. Kabsch W. Acta Crystallographica Section A. 1976;32:922.

22. Valeriani C, Allen RJ, Morelli MJ, Frenkel D, ten Wolde PR. The Journal of Chemical Physics. 2007;127:114109. [PubMed]

23. Allen RJ, Frenkel D, ten Wolde PR. The Journal of Chemical Physics. 2006;124:024102. [PubMed]

24. Moroni D, Bolhuis PG, van Erp TS. Journal of Chemical Physics. 2004;120:4055. [PubMed]

25. Warmflash A, Bhimalapuram P, Dinner AR. The Journal of Chemical Physics. 2007;127:154112. [PubMed]

26. Dickson A, Warmflash A, Dinner AR. The Journal of Chemical Physics. 2009;130:074104. [PubMed]

27. Vanden-Eijnden E, Venturoli M. The Journal of Chemical Physics. 2009;131:044120. [PubMed]

28. Allen RJ, Frenkel D, ten Wolde PR. Journal of Chemical Physics. 2006;124:17. [PubMed]

29. Zhang BW, Jasnow D, Zuckerman DM. Journal of Chemical Physics. 2010;132 [PubMed]

30. Noe F, Schutte C, Vanden-Eijnden E, Reich L, Weikl TR. Proceedings of the National Academy of Sciences of the United States of America. 2009;106:19011. [PubMed]

31. Shaw DE, Deneroff MM, Dror RO, Kuskin JS, Larson RH, Salmon JK, Young C, Batson B, Bowers KJ, Chao JC, Eastwood MP, Gagliardo J, Grossman JP, Ho CR, Ierardi DJ, Kolossvary I, Klepeis JL, Layman T, McLeavey C, Moraes MA, Mueller R, Priest EC, Shan YB, Spengler J, Theobald M, Towles B, Wang SC. Communications of the Acm. 2008;51:91.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |