Single molecule pulling with large time steps

Implementation of Telescoping Boxes in Adaptive Steered Molecular Dynamics

Long-time dynamical processes, such as those involving protein unfolding and ligand interactions, can be accelerated and realized through steered molecular dynamics (SMD). The challenge has been the extraction of information from such simulations that generalize for complex nonequilibrium processes. The use of Jarzynski's equality opened the possibility of determining the free energy along the steered coordinate, but sampling over the nonequilibrium trajectories is slow to converge. Adaptive steered molecular dynamics (ASMD) and other related techniques have been introduced to overcome this challenge through the use of stages. Here, we take advantage of these stages to address the numerical cost that arises from the required use of very large solvent boxes. We introduce telescoping box schemes within adaptive steered molecular dynamics (ASMD) in which we adjust the solvent box between stages and thereby vary (and optimize) the required number of solvent molecules. We have benchmarked the method on a relatively long α-helical peptide, Ala₃₀, with respect to the potential of mean force and hydrogen bonds. We show that the use of telescoping boxes introduces little numerical error while significantly reducing the computational cost.

Orientational and Folding Thermodynamics via Electric Dipole Moment Restraining

The projection of molecular processes onto a small set of relevant descriptors, the so-called reaction coordinates or collective variables (CVs), is a technique nowadays routinely employed by the biomolecular simulation community. In this work, we implemented two CVs to manipulate the orientation (i.e., angle) (μ⃗_a) and magnitude (|μ⃗|) of the electric dipole moment. In doing so, we studied the thermodynamics of water orientation under the application of external voltages and the folding of two polypeptides at zero-field conditions. The projection of the free-energy [potential of mean force (PMF)] along water orientation defined an upper limit of around 0.3 V for irrelevant thermodynamic effects. On the other hand, sufficiently strong μ⃗_a restraints applied on 12-alanine (Ala₁₂) triggered structural effects because of the alignment of local dipoles; for lower restraints, a full-body rotation is achieved. The manipulation of |μ⃗| produced strong perturbations on the secondary structure of Ala₁₂, promoting an enhanced sampling to its configurational space. Rigorous free-energy calculations in the form of 2-D PMFs for deca-alanine showed the utility of |μ⃗| as a reaction coordinate to study folding in small α helices. As a whole, we propose that the manipulation of both components of the dipole moment, μ⃗_a and |μ⃗|, provides thermodynamics insights into the structural conformation and stability of biomolecules. These new CVs are implemented in the Colvars module, available for NAMD and LAMMPS.

Constrained Unfolding of a Helical Peptide: Implicit versus Explicit Solvents

Steered Molecular Dynamics (SMD) has been seen to provide the potential of mean force (PMF) along a peptide unfolding pathway effectively but at significant computational cost, particularly in all-atom solvents. Adaptive steered molecular dynamics (ASMD) has been seen to provide a significant computational advantage by limiting the spread of the trajectories in a staged approach. The contraction of the trajectories at the end of each stage can be performed by taking a structure whose nonequilibrium work is closest to the Jarzynski average (in naive ASMD) or by relaxing the trajectories under a no-work condition (in full-relaxation ASMD--namely, FR-ASMD). Both approaches have been used to determine the energetics and hydrogen-bonding structure along the pathway for unfolding of a benchmark peptide initially constrained as an α-helix in a water environment. The energetics are quite different to those in vacuum, but are found to be similar between implicit and explicit solvents. Surprisingly, the hydrogen-bonding pathways are also similar in the implicit and explicit solvents despite the fact that the solvent contact plays an important role in opening the helix.

On the improvement of free-energy calculation from steered molecular dynamics simulations using adaptive stochastic perturbation protocols

The potential of mean force (PMF) calculation in single molecule manipulation experiments performed via the steered molecular dynamics (SMD) technique is a computationally very demanding task because the analyzed system has to be perturbed very slowly to be kept close to equilibrium. Faster perturbations, far from equilibrium, increase dissipation and move the average work away from the underlying free energy profile, and thus introduce a bias into the PMF estimate. The Jarzynski equality offers a way to overcome the bias problem by being able to produce an exact estimate of the free energy difference, regardless of the perturbation regime. However, with a limited number of samples and high dissipation the Jarzynski equality also introduces a bias. In our previous work, based on the Brownian motion formalism, we introduced three stochastic perturbation protocols aimed at improving the PMF calculation with the Jarzynski equality in single molecule manipulation experiments and analogous computer simulations. This paper describes the PMF reconstruction results based on full-atom molecular dynamics simulations, obtained with those three protocols. We also want to show that the protocols are applicable with the second-order cumulant expansion formula. Our protocols offer a very noticeable improvement over the simple constant velocity pulling. They are able to produce an acceptable estimate of PMF with a significantly reduced bias, even with very fast perturbation regimes. Therefore, the protocols can be adopted as practical and efficient tools for the analysis of mechanical properties of biological molecules.

Hysteresis and nonequilibrium work theorem for DNA unzipping

We study by using Monte Carlo simulations the hysteresis in unzipping and rezipping of a double stranded DNA (dsDNA) by pulling its strands in opposite directions in the fixed force ensemble. The force is increased at a constant rate from an initial value g(0) to some maximum value g(m) that lies above the phase boundary and then decreased back again to g(0). We observed hysteresis during a complete cycle of unzipping and rezipping. We obtained probability distributions of work performed over a cycle of unzipping and rezipping for various pulling rates. The mean of the distribution is found to be close (the difference being within 10%, except for very fast pulling) to the area of the hysteresis loop. We extract the equilibrium force versus separation isotherm by using the work theorem on repeated nonequilibrium force measurements. Our method is capable of reproducing the equilibrium and the nonequilibrium force-separation isotherms for the spontaneous rezipping of dsDNA.

Helix propensities calculations for amino acids in alanine based peptides using Jarzynski's equality

Jarzynski's equality (Jarzynski, Phys Rev E 1997; 56:5018 and Jarzynski, Phys Rev Lett 1997; 78:2690) relates equilibrium free energy differences between two states A and B to the work done when the system is driven repeatedly and irreversibly from an equilibrium state A to equilibrium state B. We present calculations of helix propensities using a novel procedure based on this equality. In particular, a work probability distribution is built based on combinations of multi-step trajectories that give representative work distributions without requiring computing an unreasonable large number of trajectories between states. A small number of trajectories (15) are used to construct a distribution that contains millions of work values. This distribution is used to calculate DeltaG(AB) using Jarsynski's equation. To apply and test this method, we used as a model system a dodeca-alanine helix, analyzing its extension using mechanical force. This helix was used as the basis of a host guest system in which two of the 12 residues are substituted by some other amino acid (as the guests). The differences between the unfolding free energies of the substituted peptides and the all-alanine peptide provided values for DeltaDeltaG that can be interpreted as the helix propensities of each amino acid. Results show good correlation with the experimental measurements of Baldwin and coworkers (Chakrabartty et al., Protein Sci 1994; 3:843-852).

Measuring the convergence of Monte Carlo free-energy calculations

The nonequilibrium work fluctuation theorem provides the way for calculations of (equilibrium) free-energy based on work measurements of nonequilibrium, finite-time processes, and their reversed counterparts by applying Bennett's acceptance ratio method. A nice property of this method is that each free-energy estimate readily yields an estimate of the asymptotic mean square error. Assuming convergence, it is easy to specify the uncertainty of the results. However, sample sizes have often to be balanced with respect to experimental or computational limitations and the question arises whether available samples of work values are sufficiently large in order to ensure convergence. Here, we propose a convergence measure for the two-sided free-energy estimator and characterize some of its properties, explain how it works, and test its statistical behavior. In total, we derive a convergence criterion for Bennett's acceptance ratio method.

Optimum protocol for fast-switching free-energy calculations

Free-energy differences computed from fast-switching simulations or measurements according to the Jarzynski equation are independent of the particular protocol specifying how the control parameter is changed in time. In contrast, the average work carried out on the system as well the accuracy of the resulting free energy strongly depend on the protocol. Recently, Schmiedl and Seifert [Phys. Rev. Lett. 98, 108301 (2007)] found that protocols that minimize the average work for a given duration of the switching process have discrete steps at the beginning and the end. Here we determine numerically the protocols that minimize the statistical error in the free energy estimate and find that such minimum error protocols have similar discrete jumps. Our analysis shows that the reduction in computational effort achieved by the use of steplike protocols can be considerable. Such large savings of computing time, however, typically occur for parameter ranges in which an application of the Jarzynski equation is impractical due to large statistical errors arising from the exponential work average.

Characteristic of Bennett's acceptance ratio method

A powerful and well-established tool for free-energy estimation is Bennett's acceptance ratio method. Central properties of this estimator, which employs samples of work values of a forward and its time-reversed process, are known: for given sets of measured work values, it results in the best estimate of the free-energy difference in the large sample limit. Here we state and prove a further characteristic of the acceptance ratio method: the convexity of its mean-square error. As a two-sided estimator, it depends on the ratio of the numbers of forward and reverse work values used. Convexity of its mean-square error immediately implies that there exists a unique optimal ratio for which the error becomes minimal. Further, it yields insight into the relation of the acceptance ratio method and estimators based on the Jarzynski equation. As an application, we study the performance of a dynamic strategy of sampling forward and reverse work values.

Using bijective maps to improve free-energy estimates

We derive a fluctuation theorem for generalized work distributions, related to bijective mappings of the phase spaces of two physical systems, and use it to derive a two-sided constraint maximum likelihood estimator of their free-energy difference which uses samples from the equilibrium configurations of both systems. As an application, we evaluate the chemical potential of a dense Lennard-Jones fluid and study the construction and performance of suitable maps.