Takuya Takahashi

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.

On the use of mass scaling for stable and efficient simulated tempering with molecular dynamics

Simulated tempering (ST) is a generalized‐ensemble algorithm that employs trajectories exploring a range of temperatures to effectively sample rugged energy landscapes. When implemented using the molecular dynamics method, ST can require the use of short time steps for ensuring the stability of trajectories at high temperatures. To address this shortcoming, a mass‐scaling ST (MSST) method is presented in which the particle mass is scaled in proportion to the temperature. Mass scaling in the MSST method leads to velocity distributions that are independent of temperature and eliminates the need for velocity scaling after the accepted temperature updates that are required in conventional ST simulations. The homogeneity in time scales with changing temperature improves the stability of simulations and allows for the use of longer time steps at high temperatures. As a result, the MSST is found to be more efficient than the standard ST method, particularly for cases in which a large temperature range is employed.

Exhaustive Characterization of TCRâÂÂpMHC Binding Energy Estimated by the String Model and Miyazawa-Jernigan Matrix

Accurate calculations of protein–protein binding free energies based on rigorous models, which consider the binding complex structure in atomic detail, are computationally expensive and impracticable to apply to T cell repertoire formation that occurs in the thymus because this process involves the interactions among numerous combinations of T cell receptors (TCRs) and presented peptides. By comparison, an evaluation of binding free-energy using a combination of the string model and Miyazawa-Jernigan matrix is very efficient and was therefore applied to estimate interaction energies between T cell receptor–peptide–MHC (TCR–pMHC) complexes, which appeared to successfully explain the effects of binding capacity of MHC on repertoire–formation and the reason for the presence of elite-controllers of some viral infections. However, this evaluation method is overly simplified and requires more detailed considerations when applied to evaluating TCR-pMHC interactions. In this study, we examined this method exhaustively and revealed the limitations of the method. Following features necessitate cautious attitude when interpreting the calculation results: first, the apparent increase in the number of hot spots in accordance with an increase of educational epitope pool size does not mean an increased TCR specificity of surviving clones; second, strong binders to any TCR converge to some limited sequences that are determined by the physical nature of the Miyazawa-Jernigan matrix.

Momentum and Velocity Scaling Rules in Replica-Exchange Molecular Dynamics Simulations with Mass Manipulation

The authors have recently presented the mass-scaling replica-exchange molecular dynamics (MSREMD) method [J. Chem. Phys. 141, 114111 (2014)], in which all masses are scaled in proportion to temperature for better numerical stability with a large time step. In addition, the scaling of masses, performed in advance of the simulation, can be substituted for the velocity scaling necessary after every accepted replica-exchange attempt, and the algorithm can thereby be simplified. In this work, the authors present the mass-manipulating REMD (MMREMD) method, where arbitrary mass scaling is employed. Rules for momentum and velocity scaling are formalized for the MMREMD method, followed by a demonstration with two replicas. Adherence to these rules is crucial for sampling of the correct canonical distribution. The authors also review and discuss a previous study as a sample application of the MMREMD method.

Unfolding of α-helical 20-residue poly-glutamic acid analyzed by multiple runs of canonical molecular dynamics simulations

Elucidating the molecular mechanism of helix–coil transitions of short peptides is a long-standing conundrum in physical chemistry. Although the helix–coil transitions of poly-glutamic acid (PGA) have been extensively studied, the molecular details of its unfolding process still remain unclear. We performed all-atom canonical molecular dynamics simulations for a 20-residue PGA, over a total of 19 μs, in order to investigate its helix-unfolding processes in atomic resolution. Among the 28 simulations, starting with the α-helical conformation, all showed an unfolding process triggered by the unwinding of terminal residues, rather than by kinking and unwinding of the middle region of the chain. The helix–coil–helix conformation which is speculated by the previous experiments was not observed. Upon comparison between the N- and C-termini, the latter tended to be unstable and easily unfolded. While the probabilities of helix elongation were almost the same among the N-terminal, middle, and C-terminal regions of the chain, unwinding of the helix was enriched at the C-terminal region. The turn and 310-helix conformations were kinetic intermediates in the formation and deformation of α-helix, consistent with the previous computational studies for Ala-based peptides.

Influence of various parameters in the replica-exchange molecular dynamics method: Number of replicas, replica-exchange frequency, and thermostat coupling time constant

The replica-exchange molecular dynamics (REMD) method has been used for conformational sampling of various biomolecular systems. To maximize sampling efficiency, some adjustable parameters must be optimized. Although it is agreed that shorter intervals between the replica-exchange attempts enhance traversals in the temperature space, details regarding the artifacts caused by these short intervals are controversial. In this study, we revisit this problem by performing REMD simulations on an alanine octapeptide in an implicit solvent. Fifty different sets of conditions, which are a combination of five replica-exchange periods, five different numbers of replicas, and two thermostat coupling time constants, were investigated. As a result, although short replica-exchange intervals enhanced the traversals in the temperature space, they led to artifacts in the ensemble average of the temperature, potential energy, and helix content. With extremely short replica-exchange intervals, i.e., attempted at every time step, the ensemble average of the temperature deviated from the thermostat temperature by ca. 7 K. Differences in the ensembles were observed even for larger replica-exchange intervals (between 100 and 1,000 steps). In addition, the shorter thermostat coupling time constant reduced the artifacts found when short replica-exchange intervals were used, implying that these artifacts are caused by insufficient thermal relaxation between the replica-exchange events. Our results will be useful to reduce the artifacts found in REMD simulations by adjusting some key parameters.

Practical Estimation of TCR-pMHC Binding Free-Energy Based on the Dielectric Model and the Coarse-Grained Model

To evaluate free energy changes of bio-molecules in a water solution, ab initio molecular dynamics (MD) simulations such as Quantum Mechanical Molecular Mechanics (QM/MM) and MD are the most theoretically rigorous methods (Car and Parrinello 1985; Kuhne, Krack et al. 2007), although the calculation cost is far too large for large molecular systems that contain many electrons. Therefore, all-atom MD simulations based on classical mechanics (i.e., Newton’s equations) are used for the usual bio-molecular systems. As the conventional free energy perturbation (FEP) method based on all atom MD simulation is a strict method, to elucidate the molecular principles upon which the selectivity of a TCR is based, FEP simulations are used to analyse the binding free energy difference of a particular TCR (A6) for a wild-type peptide (Tax) and a mutant peptide (Tax P6A), both presented in HLA A2. The computed free energy difference is 2.9 kcal mol-1 and the agreement with the experimental value is good, although the calculation is very time-consuming and the simulation time is still insufficient for fully sampling the phase space. From this FEP calculation, better solvation of the mutant peptide when bound to the MHC molecule is important to the greater affinity of the TCR for the latter. This suggests that the exact and efficient evaluation of solvation is important for the affinity calculation (Michielin and Karplus 2002). Other FEP calculations of the wild-type and the variant human T cell lymphotropic virus type 1 Tax peptide presented by the MHC to the TCR have been performed using large scale massively parallel molecular dynamics simulations and the computed free energy difference using alchemical mutationbased thermodynamic integration, which agrees well with experimental data semiquantitatively (Wan, Coveney et al. 2005). However, the conventional FEP is still very timeconsuming when searching for so many unknown docking structures because all-atom MD for a large molecular system is a computationally hard task and MD simulations must be done not only in initial and final states but also in many intermediate states.

A Numerical Approach to Evaluate Electrostatic Stabilization of Protein Crystals

We developed a novel algorithm to numerically solve the Poisson-Boltzmann equations under a periodic boundary condition. By employing this algorithm to calculate the electrostatic potentials in two different types of protein crystals, a bovine-pancreatic-trypsin-inhibiter (BPTI) orthorhombic crystal and a pig-insulin cubic crystal, the energy contributions of the electrostatic interactions to the crystals stability were evaluated. At a high ionic strength, the condensed state of proteins in the crystal was electrostatically stabilized compared with that isolated in dilute solution. At a low ionic strength, on the other hand, the electrostatic interactions destabilized the crystalline state of both proteins, although the ionic strength dependence was different. In all of the solvent ionic strengths investigated, the attractive electrostatic interactions between charge pairs separated by less than 5 A on the respective protein molecules prominently stabilize the protein crystals. We also found a specific role for bound phosphate ions in the stabilization of the BPTI crystal based on comparison of the electrostatic energies of the two crystals with and without the ions.