Replica-Permutation Method with the Suwa-Todo Algorithm beyond the Replica-Exchange Method

Enhanced sampling simulations to construct free-energy landscape of protein-partner substrate interaction

Molecular dynamics (MD) simulations using all-atom and explicit solvent models provide valuable information on the detailed behavior of protein-partner substrate binding at the atomic level. As the power of computational resources increase, MD simulations are being used more widely and easily. However, it is still difficult to investigate the thermodynamic properties of protein-partner substrate binding and protein folding with conventional MD simulations. Enhanced sampling methods have been developed to sample conformations that reflect equilibrium conditions in a more efficient manner than conventional MD simulations, thereby allowing the construction of accurate free-energy landscapes. In this review, we discuss these enhanced sampling methods using a series of case-by-case examples. In particular, we review enhanced sampling methods conforming to trivial trajectory parallelization, virtual-system coupled multicanonical MD, and adaptive lambda square dynamics. These methods have been recently developed based on the existing method of multicanonical MD simulation. Their applications are reviewed with an emphasis on describing their practical implementation. In our concluding remarks we explore extensions of the enhanced sampling methods that may allow for even more efficient sampling.

Mass-scaling replica-exchange molecular dynamics optimizes computational resources with simpler algorithm

We develop a novel method of replica-exchange molecular dynamics (REMD) simulation, mass-scaling REMD (MSREMD) method, which improves numerical stability of simulations. In addition, the MSREMD method can also simplify a replica-exchange routine by eliminating velocity scaling. As a pilot system, a Lennard-Jones fluid is simulated with the new method. The results suggest that the MSREMD method improves the numerical stability at high temperatures compared with the conventional REMD method. For the Nosé-Hoover thermostats, we analytically demonstrate that the MSREMD simulations can reproduce completely the same trajectories of the conventional REMD ones with shorter time steps at high temperatures. Accordingly, we can easily compare the computational costs of the REMD and MSREMD simulations. We conclude that the MSREMD method decreases the instability and optimizes the computational resources with simpler algorithm.

Simulated tempering based on global balance or detailed balance conditions: Suwa-Todo, heat bath, and Metropolis algorithms

Simulated tempering (ST) is a useful method to enhance sampling of molecular simulations. When ST is used, the Metropolis algorithm, which satisfies the detailed balance condition, is usually applied to calculate the transition probability. Recently, an alternative method that satisfies the global balance condition instead of the detailed balance condition has been proposed by Suwa and Todo. In this study, ST method with the Suwa-Todo algorithm is proposed. Molecular dynamics simulations with ST are performed with three algorithms (the Metropolis, heat bath, and Suwa-Todo algorithms) to calculate the transition probability. Among the three algorithms, the Suwa-Todo algorithm yields the highest acceptance ratio and the shortest autocorrelation time. These suggest that sampling by a ST simulation with the Suwa-Todo algorithm is most efficient. In addition, because the acceptance ratio of the Suwa-Todo algorithm is higher than that of the Metropolis algorithm, the number of temperature states can be reduced by 25% for the Suwa-Todo algorithm when compared with the Metropolis algorithm.

Dimerization process of amyloid-β(29-42) studied by the Hamiltonian replica-permutation molecular dynamics simulations

The amyloid-β peptides form amyloid fibrils which are associated with Alzheimer's disease. Amyloid-β(29-42) is its C-terminal fragment and a critical determinant of the amyloid formation rate. This fragment forms the amyloid fibril by itself. However, the fragment conformation in the fibril has yet to be determined. The oligomerization process including the dimerization process is also still unknown. The dimerization process corresponds to an early process of the amyloidogenesis. In order to investigate the dimerization process and conformations, we applied the Hamiltonian replica-permutation method, which is a better alternative to the Hamiltonian replica-exchange method, to two amyloid-β(29-42) molecules in explicit water solvent. At the first step of the dimerization process, two amyloid-β(29-42) molecules came close to each other and had intermolecular side chain contacts. When two molecules had the intermolecular side chain contacts, the amyloid-β(29-42) tended to have intramolecular secondary structures, especially β-hairpin structures. The two molecules had intermolecular β-bridge structures by coming much closer at the second step of the dimerization process. Formation of these intermolecular β-bridge structures was induced by the β-hairpin structures. The intermolecular β-sheet structures elongated at the final step. Structures of the amyloid-β(29-42) in the monomer and dimer states are also shown with the free-energy landscapes, which were obtained by performing efficient sampling in the conformational space in our simulations.

Decomposition-order effects of time integrator on ensemble averages for the Nosé-Hoover thermostat

Decomposition-order dependence of time development integrator on ensemble averages for the Nosé-Hoover dynamics is discussed. Six integrators were employed for comparison, which were extensions of the velocity-Verlet or position-Verlet algorithm. Molecular dynamics simulations by these integrators were performed for liquid-argon systems with several different time steps and system sizes. The obtained ensemble averages of temperature and potential energy were shifted from correct values depending on the integrators. These shifts increased in proportion to the square of the time step. Furthermore, the shifts could not be removed by increasing the number of argon atoms. We show the origin of these ensemble-average shifts analytically. Our discussion can be applied not only to the liquid-argon system but also to all MD simulations with the Nosé-Hoover thermostat. Our recommended integrators among the six integrators are presented to obtain correct ensemble averages.

Transformation between α-helix and β-sheet structures of one and two polyglutamine peptides in explicit water molecules by replica-exchange molecular dynamics simulations

Aggregation of polyglutamine peptides with β-sheet structures is related to some important neurodegenerative diseases such as Huntington's disease. However, it is not clear how polyglutamine peptides form the β-sheets and aggregate. To understand this problem, we performed all-atom replica-exchange molecular dynamics simulations of one and two polyglutamine peptides with 10 glutamine residues in explicit water molecules. Our results show that two polyglutamine peptides mainly formed helix or coil structures when they are separated, as in the system with one-polyglutamine peptide. As the interpeptide distance decreases, the intrapeptide β-sheet structure sometimes appear as an intermediate state, and finally the interpeptide β-sheets are formed. We also find that the polyglutamine dimer tends to form the antiparallel β-sheet conformations rather than the parallel β-sheet, which is consistent with previous experiments and a coarse-grained molecular dynamics simulation.

Enhanced exchange algorithm without detailed balance condition for replica exchange method

The replica exchange method (REM) is a powerful tool for the conformational sampling of biomolecules. In this study, we propose an enhanced exchange algorithm for REM not meeting the detailed balance condition (DBC), but satisfying the balance condition in all considered exchanges between two replicas. Breaking the DBC can minimize the rejection rate and make an exchange process rejection-free as the number of replicas increases. To enhance the efficiency of REM, all possible pairs--not only the nearest neighbor--were considered in the exchange process. The test simulations of the alanine dipeptide confirmed the correctness of our method. The average traveling distance of each replica in the temperature distribution was also increased in proportion to an increase in the exchange rate. Furthermore, we applied our algorithm to the conformational sampling of the 10-residue miniprotein, chignolin, with an implicit solvent model. The results showed a faster convergence in the calculation of its free energy landscape, compared to that achieved using the normal exchange method of adjacent pairs. This algorithm can also be applied to the conventional near neighbor method and is expected to reduce the required number of replicas.

Hamiltonian replica-permutation method and its applications to an alanine dipeptide and amyloid-β(29-42) peptides

We propose the Hamiltonian replica-permutation method (RPM) (or multidimensional RPM) for molecular dynamics and Monte Carlo simulations, in which parameters in the Hamiltonian are permuted among more than two replicas with the Suwa-Todo algorithm. We apply the Coulomb RPM, which is one of realization of the Hamiltonian RPM, to an alanine dipeptide and to two amyloid-β(29-42) molecules. The Hamiltonian RPM realizes more efficient sampling than the Hamiltonian replica-exchange method. We illustrate the protein misfolding funnel of amyloid-β(29-42) and reveal its dimerization pathways.

Transformation of a design peptide between the α-helix and β-hairpin structures using a helix-strand replica-exchange molecular dynamics simulation

We investigated the transformation between the α-helix and β-hairpin structures of an 18-residue design peptide, whose sequence is INYWLAHAKAGYIVHWTA. This peptide has both α-helix and β-hairpin structures in aqueous solution. For this purpose, we proposed the helix-strand replica-exchange method. This is one of the Hamiltonian replica-exchange methods in which we exchange parameters for umbrella potentials to enhance the α-helix or β-strand structure formation. We performed an all-atom helix-strand replica-exchange molecular dynamics (MD) simulation of this peptide in explicit water solvent with five replicas. Because the suitable umbrella potential was applied, the helix-strand replica-exchange MD simulation reproduced conformations closer to experimental conformations than a temperature replica-exchange MD simulation when the same numbers of the replicas were used, while the temperature replica-exchange MD simulation does not require bias along any specific order parameter. We calculated its free-energy landscape and revealed the transformation pathways between the α-helix and β-hairpin structures and the folding pathways from an extended structure. Although the fractions of the α-helix and β-hairpin structures are less than those obtained by the experiment, the free-energy difference between the two structures is calculated to be almost zero, which agrees with the experimental results.

Pressure-Induced Helical Structure of a Peptide Studied by Simulated Tempering Molecular Dynamics Simulations

It is known experimentally that an AK16 peptide forms more α-helix structures with increasing pressure while proteins unfold in general. In order to understand this abnormality, molecular dynamics (MD) simulations with the simulated tempering method for the isobaric-isothermal ensemble were performed in a wide pressure range from 1.0 × 10(-4) GPa to 1.4 GPa. From the results of the simulations, it is found that the fraction of the folded state decreases once and increases after that with increasing pressure. The partial molar volume change from the folded state to unfolded state increases monotonically from a negative value to a positive value with pressure. The behavior under high pressure conditions is consistent with the experimental results. The radius of gyration of highly helical structures decreases with increasing pressure, which indicates that the helix structure shrinks with pressure. This is the reason why the fraction of the folded state increases as pressure increases.