Simulation of peptide folding with explicit water--a mean solvation method

Reformulation of the self-guided molecular simulation method

Self-guided molecular/Langevin dynamics (SGMD/SGLD) simulation methods were developed to enhance conformational sampling through promoting low frequency motion of molecular systems and have been successfully applied in many simulation studies. Quantitative understanding of conformational distribution in SGLD has been achieved by separating microscopic properties according to frequency. However, a missing link between the guiding factors and conformational distributions makes it highly empirical and system dependent when choosing the values of the guiding parameters. Based on the understanding that molecular interactions are the source of energy barriers and diffusion friction, this work reformulates the equation of the low frequency motion to resemble Langevin dynamics. This reformulation leads to new forms of guiding forces and establishes a relation between the guiding factors and conformational distributions. We call simulations with these new guiding forces the generalized self-guided molecular/Langevin dynamics (SGMDg/SGLDg). In addition, we present a new way to calculate low frequency properties and an efficient algorithm to implement SGMDg/SGLDg that minimizes memory usage and inter-processor communication. Through example simulations with a skewed double well system, an argon fluid, and a cryo-EM map flexible fitting case, we demonstrate the guiding effects on conformational distributions and conformational searching.

Origin of pKa Shifts of Internal Lysine Residues in SNase Studied Via Equal-Molar VMMS Simulations in Explicit Water

Protein internal ionizable groups can exhibit large shifts in pK_a values. Although the environment and interaction changes have been extensively studied both experimentally and computationally, direct calculation of pK_a values of these internal ionizable groups in explicit water is challenging due to energy barriers in solvent interaction and in conformational transition. The virtual mixture of multiple states (VMMS) method is a new approach designed to study chemical state equilibrium. This method constructs a virtual mixture of multiple chemical states in order to sample the conformational space of all states simultaneously and to avoid crossing energy barriers related to state transition. By applying VMMS to 25 variants of staphylococcal nuclease with lysine residues at internal positions, we obtained the pK_a values of these lysine residues and investigated the physics underlining the pK_a shifts. Our calculation results agree reasonably well with experimental measurements, validating the VMMS method for pK_a calculation and providing molecular details of the protonation equilibrium for protein internal ionizable groups. Based on our analyses of protein conformation relaxation, lysine side chain flexibility, water penetration, and the microenvironment, we conclude that the hydrophobicity of the microenvironment around the lysine side chain (which affects water penetration differently for different protonation states) plays an important role in the pK_a shifts.

A Virtual Mixture Approach to the Study of Multistate Equilibrium: Application to Constant pH Simulation in Explicit Water

Chemical and thermodynamic equilibrium of multiple states is a fundamental phenomenon in biology systems and has been the focus of many experimental and computational studies. This work presents a simulation method to directly study the equilibrium of multiple states. This method constructs a virtual mixture of multiple states (VMMS) to sample the conformational space of all chemical states simultaneously. The VMMS system consists of multiple subsystems, one for each state. The subsystem contains a solute and a solvent environment. The solute molecules in all subsystems share the same conformation but have their own solvent environments. Transition between states is implicated by the change of their molar fractions. Simulation of a VMMS system allows efficient calculation of relative free energies of all states, which in turn determine their equilibrium molar fractions. For systems with a large number of state transition sites, an implicit site approximation is introduced to minimize the cost of simulation. A direct application of the VMMS method is for constant pH simulation to study protonation equilibrium. Applying the VMMS method to a heptapeptide of 3 ionizable residues, we calculated the pKas of those residues both with all explicit states and with implicit sites and obtained consistent results. For mouse epidermal growth factor of 9 ionizable groups, our VMMS simulations with implicit sites produced pKas of all 9 ionizable groups and the results agree qualitatively with NMR measurement. This example demonstrates the VMMS method can be applied to systems of a large number of ionizable groups and the computational cost scales linearly with the number of ionizable groups. For one of the most challenging systems in constant pH calculation, SNase Δ+PHS/V66K, our VMMS simulation shows that it is the state-dependent water penetration that causes the large deviation in lysine66’s pKa.

Computer simulation plays an important role to understand molecular systems and has been applied to problems of increasing complexity. Multistate equilibrium is a fundamental concept behind the structure and function of biological systems. Due to the limit in computing resources and lack of good alternative methods, computer simulation has been conducted for systems in a single state, sampling from one state to another to infer equilibrium properties. This sequential approach has been successful in many cases such as protonation equilibrium with implicit solvation model. However, state transition is difficult when explicit solvent is used for more accurate solvation description. Many efforts have been dedicated to overcome this difficulty. Analogous to real multistate systems, we proposed a virtual mixture of multiple states (VMMS) to directly simulate the equilibrium. State transitions are replaced by changes in state molar fractions. Mimicking a test tube environment, all states are simulated in parallel to equilibrate with each other. Application to constant pH simulation in explicit water demonstrates the capability of this method. It is expected that the VMMS method will find more applications in biological problems related to the equilibrium of competing states.

Efficient and Unbiased Sampling of Biomolecular Systems in the Canonical Ensemble: A Review of Self-Guided Langevin Dynamics

This review provides a comprehensive description of the self-guided Langevin dynamics (SGLD) and the self-guided molecular dynamics (SGMD) methods and their applications. Example systems are included to provide guidance on optimal application of these methods in simulation studies. SGMD/SGLD has enhanced ability to overcome energy barriers and accelerate rare events to affordable time scales. It has been demonstrated that with moderate parameters, SGLD can routinely cross energy barriers of 20 kT at a rate that molecular dynamics (MD) or Langevin dynamics (LD) crosses 10 kT barriers. The core of these methods is the use of local averages of forces and momenta in a direct manner that can preserve the canonical ensemble. The use of such local averages results in methods where low frequency motion "borrows" energy from high frequency degrees of freedom when a barrier is approached and then returns that excess energy after a barrier is crossed. This self-guiding effect also results in an accelerated diffusion to enhance conformational sampling efficiency. The resulting ensemble with SGLD deviates in a small way from the canonical ensemble, and that deviation can be corrected with either an on-the-fly or a post processing reweighting procedure that provides an excellent canonical ensemble for systems with a limited number of accelerated degrees of freedom. Since reweighting procedures are generally not size extensive, a newer method, SGLDfp, uses local averages of both momenta and forces to preserve the ensemble without reweighting. The SGLDfp approach is size extensive and can be used to accelerate low frequency motion in large systems, or in systems with explicit solvent where solvent diffusion is also to be enhanced. Since these methods are direct and straightforward, they can be used in conjunction with many other sampling methods or free energy methods by simply replacing the integration of degrees of freedom that are normally sampled by MD or LD.

Force-momentum-based self-guided Langevin dynamics: a rapid sampling method that approaches the canonical ensemble

The self-guided Langevin dynamics (SGLD) is a method to accelerate conformational searching. This method is unique in the way that it selectively enhances and suppresses molecular motions based on their frequency to accelerate conformational searching without modifying energy surfaces or raising temperatures. It has been applied to studies of many long time scale events, such as protein folding. Recent progress in the understanding of the conformational distribution in SGLD simulations makes SGLD also an accurate method for quantitative studies. The SGLD partition function provides a way to convert the SGLD conformational distribution to the canonical ensemble distribution and to calculate ensemble average properties through reweighting. Based on the SGLD partition function, this work presents a force-momentum-based self-guided Langevin dynamics (SGLDfp) simulation method to directly sample the canonical ensemble. This method includes interaction forces in its guiding force to compensate the perturbation caused by the momentum-based guiding force so that it can approximately sample the canonical ensemble. Using several example systems, we demonstrate that SGLDfp simulations can approximately maintain the canonical ensemble distribution and significantly accelerate conformational searching. With optimal parameters, SGLDfp and SGLD simulations can cross energy barriers of more than 15 kT and 20 kT, respectively, at similar rates for LD simulations to cross energy barriers of 10 kT. The SGLDfp method is size extensive and works well for large systems. For studies where preserving accessible conformational space is critical, such as free energy calculations and protein folding studies, SGLDfp is an efficient approach to search and sample the conformational space.

Toward canonical ensemble distribution from self-guided Langevin dynamics simulation

This work derives a quantitative description of the conformational distribution in self-guided Langevin dynamics (SGLD) simulations. SGLD simulations employ guiding forces calculated from local average momentums to enhance low-frequency motion. This enhancement in low-frequency motion dramatically accelerates conformational search efficiency, but also induces certain perturbations in conformational distribution. Through the local averaging, we separate properties of molecular systems into low-frequency and high-frequency portions. The guiding force effect on the conformational distribution is quantitatively described using these low-frequency and high-frequency properties. This quantitative relation provides a way to convert between a canonical ensemble and a self-guided ensemble. Using example systems, we demonstrated how to utilize the relation to obtain canonical ensemble properties and conformational distributions from SGLD simulations. This development makes SGLD not only an efficient approach for conformational searching, but also an accurate means for conformational sampling.

Quantitative Assessment of Force Fields on Both Low-Energy Conformational Basins and Transition-State Regions of the (ϕ-ψ) Space

The free energy surfaces (FESs) of alanine dipeptide are studied to illustrate a new strategy to assess the performance of classical molecular mechanics force field on the full range of the (ϕ-ψ) conformational space. The FES is obtained from metadynamics simulations with five commonly used force fields and from ab initio density functional theory calculations in both gas phase and aqueous solution. The FESs obtained at the B3LYP/6-311+G(2d,p)//B3LYP/6-31G(d,p) level of theory are validated by comparison with previously reported MP2 and LMP2 results as well as with experimentally obtained probability distribution between the C5-β (or β-PPII) and αR states. A quantitative assessment is made for each force field in three conformational basins, LeRI (C5-β-C7eq), LeRII (β2-αR), and LeRIII(αL-C7ax-αD) as well as three transition-state regions linking the above conformational basins. The performance of each force field is evaluated in terms of the average free energy of each region in comparison with that of the ab initio results. We quantify how well a force field FES matches the ab initio FES through the calculation of the standard deviation of a free energy difference map between the two FESs. The results indicate that the performance varies largely from region to region or from force field to force field. Although not one force field is able to outperform all others in all conformational areas, the OPLSAA/L force field gives the best performance overall, followed by OPLSAA and AMBER03. For the three top performers, the average free energies differ from the corresponding ab initio values from within the error range (<0.4 kcal/mol) to ∼1.5 kcal/mol for the low-energy regions and up to ∼2.0 kcal/mol for the transition-state regions. The strategy presented and the results obtained here should be useful for improving the parametrization of force fields targeting both accuracy in the energies of conformers and the transition-state barriers.

Monte carlo folding of trans-membrane helical peptides in an implicit generalized Born membrane

An efficient Monte Carlo (MC) algorithm using concerted backbone rotations is combined with a recently developed implicit membrane model to simulate the folding of the hydrophobic transmembrane domain M2TM of the M2 protein from influenza A virus and Sarcolipin at atomic resolution. The implicit membrane environment is based on generalized Born theory and has been calibrated against experimental data. The MC sampling has previously been used to fold several small polypeptides and been shown to be equivalent to molecular dynamics (MD). In combination with a replica exchange algorithm, M2TM is found to form continuous membrane spanning helical conformations for low temperature replicas. Sarcolipin is only partially helical, in agreement with the experimental NMR structures in lipid bilayers and detergent micelles. Higher temperature replicas exhibit a rapidly decreasing helicity, in agreement with expected thermodynamic behavior. To exclude the possibility of an erroneous helical bias in the simulations, the model is tested by sampling a synthetic Alanine-rich polypeptide of known helicity. The results demonstrate there is no overstabilization of helical conformations, indicating that the implicit model captures the essential components of the native membrane environment for M2TM and Sarcolipin.

Breaking non-native hydrophobic clusters is the rate-limiting step in the folding of an alanine-based peptide

The formation mechanism of an alanine-based peptide has been studied by all-atom molecular dynamics simulations with a recently developed all-atom point-charge force field and the Generalize Born continuum solvent model at an effective salt concentration of 0.2M. Thirty-two simulations were conducted. Each simulation was performed for 100 ns. A surprisingly complex folding process was observed. The development of the helical content can be divided into three phases with time constants of 0.06-0.08, 1.4-2.3, and 12-13 ns, respectively. Helices initiate extreme rapidly in the first phase similar to that estimated from explicit solvent simulations. Hydrophobic collapse also takes place in this phase. A folding intermediate state develops in the second phase and is unfolded to allow the peptide to reach the transition state in the third phase. The folding intermediate states are characterized by the two-turn short helices and the transition states are helix-turn-helix motifs-both of which are stabilized by hydrophobic clusters. The equilibrium helical content, calculated by both the main-chain Phi-Psi torsion angles and the main-chain hydrogen bonds, is 64-66%, which is in remarkable agreement with experiments. After corrected for the solvent viscosity effect, an extrapolated folding time of 16-20 ns is obtained that is in qualitative agreement with experiments. Contrary to the prevailing opinion, neither initiation nor growth of the helix is the rate-limiting step. Instead, the rate-limiting step for this peptide is breaking the non-native hydrophobic clusters in order to reach the transition state. The implication to the folding mechanisms of proteins is also discussed.

Molecular dynamics simulations of beta-hairpin folding

Molecular dynamics simulations of beta-hairpin folding have been carried out with a solvent-referenced potential at 274 K. The model peptide V4DPGV4 formed stable beta-hairpin conformations and the beta-hairpin ratio calculated by the DSSP algorithm was about 56% in the 50-ns simulation. Folding into beta-hairpin conformations is independent of the initial conformations. The simulations provided insights into the folding mechanism. The hydrogen bond often formed in a beta-turn first, and then propagated by forming more hydrogen bonds along the strands. Unfolding and refolding occurred repeatedly during the simulations. Both the hydrogen bonding and the hydrophobic interaction played important roles in forming the ordered structure. Without the hydrophobic effect, stable beta-hairpin conformations did not form in the simulations. With the same energy functions, the alanine-based peptide (AAQAA)3Y folded into helical conformations, in agreement with experiments. Folding into an alpha-helix or a beta-hairpin is amino acid sequence-dependent.