Ewald summation and multiple time step methods for molecular dynamics simulation of biological molecules

Methods by which to determine conditions for a molecular dynamics (MD) simulation of biological molecules were investigated. Derivation of the optimal parameters of the Ewald summation was described so as to give same precision to the real space, the reciprocal space summations and the van der Waals interaction. Later, the procedure by which to determine the condition of the multiple time step method by RESPA (REference System Propagator Algorithm; Tuckerman et al., 1992, J. Chem. Phys., 97, 1990) was described as exemplified by MD simulations of a solvated β-sheet peptide. The conservation of the total energy in a microcanonical ensemble was measured to investigate the stability of the simulation conditions. The most feasible respective combinations of the time steps were: 0.25 fs for bond, angle and torsion interactions; 2 fs for van der Waals interaction and Ewald real-space summation; and 4 fs for Ewald reciprocal-space summation. Though it retained an acceptable accuracy, this condition accelerated the simulation ten-fold compared to that in which a simple velocity-Verlet method with a time step of 0.25 fs was used. The update of the correction term due to excluded neighbors was then investigated. Better results were obtained when the correction was updated with the real-space than when it was updated with the reciprocal-space summation. Finally, an MD simulation as long as 50 ps performed under the optimal Ewald and RESPA parameters was thus determined. The trajectory showed a good stability, indicating the feasibility of the parameters.     

A multiple time step algorithm compatible with a large number of distance classes and an arbitrary distance dependence of the time step size for the fast evaluation of nonbonded interactions in molecular simulations

A new algorithm is introduced to perform the multiple time step integration of the equations of motion for a molecular system, based on the splitting of the nonbonded interactions into a series of distance classes. The interactions between particle pairs in successive classes are updated at a progressively decreasing frequency. Unlike previous multiple time-stepping schemes relying on distance classes, the present algorithm sorts interacting particle pairs by their next update times rather than by their update frequencies. For this reason, the proposed scheme is extremely flexible with respect to the number of classes that can be employed (up to hundred or more) and the distance dependence of the relative time step size (arbitrary integer function of the distance). It can also easily be adapted to classes defined based on a criterion other than the interparticle distance (e.g., interaction magnitude). Different variants of the algorithm are tested in terms of accuracy and efficiency for simulations of a pure water system (6167 molecules) under truncated-octahedral periodic boundary conditions, and compared to the twin-range method standardly used with GROMOS96 (short- and long-range cutoff distances of 0.8 and 1.4 nm, pair list and intermediate-range interactions updated every five steps). In particular, multiple time-stepping schemes with an accuracy comparable to that of the twin-range method can be designed, that permit to increase the effective (long-range) cutoff distance from 1.4 to 3.0 nm with a performance loss of only about a factor 2. This result is quite encouraging, considering the benefits of doubling the cutoff radius in the context of (bio-)molecular simulations.     

Computationally efficient canonical molecular dynamics simulations by using a multiple time‐step integrator algorithm combined with the particle mesh Ewald method and with the fast multipole method

An efficient implementation of the canonical molecular dynamics simulation using the reversible reference system propagator algorithm (r‐RESPA) combined with the particle mesh Ewald method (PMEM) and with the macroscopic expansion of the fast multipole method (MEFMM) was examined. The performance of the calculations was evaluated for systems with 3000, 9999, 30,000, 60,000, and 99,840 particles. For a given accuracy, the optimal conditions for minimizing the CPU time for the implementation of the Ewald method, the PMEM, and the MEFMM were first analyzed. Using the optimal conditions, we evaluated the performance and the reliability of the integrated methods. For all the systems examined, the r‐RESPA with the PMEM was about twice as fast as the r‐RESPA with the MEFMM. The difference arose from the difference in the numerical complexities of the fast Fourier transform in the PMEM and from the transformation of the multipole moments into the coefficients of the local field expansion in the MEFMM. Compared with conventional methods, such as the velocity‐verlet algorithm with the Ewald method, significant speedups were obtained by the integrated methods; the speedup of the calculation was a function of system size, and was a factor of 100 for a system with 3000 particles and increased to a factor of 700 for a system with 99,840 particles. These integrated calculations are, therefore, promising for realizing large‐scale molecular dynamics simulations for complex systems. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 201–217, 2000     

GENESIS 1.1: A hybrid‐parallel molecular dynamics simulator with enhanced sampling algorithms on multiple computational platforms

GENeralized‐Ensemble SImulation System (GENESIS) is a software package for molecular dynamics (MD) simulation of biological systems. It is designed to extend limitations in system size and accessible time scale by adopting highly parallelized schemes and enhanced conformational sampling algorithms. In this new version, GENESIS 1.1, new functions and advanced algorithms have been added. The all‐atom and coarse‐grained potential energy functions used in AMBER and GROMACS packages now become available in addition to CHARMM energy functions. The performance of MD simulations has been greatly improved by further optimization, multiple time‐step integration, and hybrid (CPU + GPU) computing. The string method and replica‐exchange umbrella sampling with flexible collective variable choice are used for finding the minimum free‐energy pathway and obtaining free‐energy profiles for conformational changes of a macromolecule. These new features increase the usefulness and power of GENESIS for modeling and simulation in biological research. © 2017 Wiley Periodicals, Inc.     

Molecule‐specific determination of atomic polarizabilities with the polarizable atomic multipole model

Recently, many polarizable force fields have been devised to describe induction effects between molecules. In popular polarizable models based on induced dipole moments, atomic polarizabilities are the essential parameters and should be derived carefully. Here, we present a parameterization scheme for atomic polarizabilities using a minimization target function containing both molecular and atomic information. The main idea is to adopt reference data only from quantum chemical calculations, to perform atomic polarizability parameterizations even when relevant experimental data are scarce as in the case of electronically excited molecules. Specifically, our scheme assigns the atomic polarizabilities of any given molecule in such a way that its molecular polarizability tensor is well reproduced. We show that our scheme successfully works for various molecules in mimicking dipole responses not only in ground states but also in valence excited states. The electrostatic potential around a molecule with an externally perturbing nearby charge also exhibits a near‐quantitative agreement with the reference data from quantum chemical calculations. The limitation of the model with isotropic atoms is also discussed to examine the scope of its applicability. © 2012 Wiley Periodicals, Inc.     

Increasing the time step with mass scaling in Born‐Oppenheimer ab initio QM/MM molecular dynamics simulations

Born‐Oppenheimer ab initio QM/MM molecular dynamics simulation with umbrella sampling is a state‐of‐the‐art approach to calculate free energy profiles of chemical reactions in complex systems. To further improve its computational efficiency, a mass‐scaling method with the increased time step in MD simulations has been explored and tested. It is found that by increasing the hydrogen mass to 10 amu, a time step of 3 fs can be employed in ab initio QM/MM MD simulations. In all our three test cases, including two solution reactions and one enzyme reaction, the resulted reaction free energy profiles with 3 fs time step and mass scaling are found to be in excellent agreement with the corresponding simulation results using 1 fs time step and the normal mass. These results indicate that for Born‐Oppenheimer ab initio QM/MM molecular dynamics simulations with umbrella sampling, the mass‐scaling method can significantly reduce its computational cost while has little effect on the calculated free energy profiles. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009     

Influence of variation of a side chain on the folding equilibrium of a β‐peptide: Limitations of one‐step perturbation

In a recent study (Lin et al., Helv. Chim. Acta 2011, 94 , 597), the one‐step perturbation method was applied to tackle a challenging computational problem, that is, the calculation of the folding free enthalpies ΔG_F,U of six hepta‐β‐peptides with different, Ala, Val, Leu, Ile, Ser, or Thr, side chains in the fifth residue. The ΔG_F,U values obtained using one‐step perturbation based on a single molecular dynamics simulation of a judiciously chosen reference state with soft‐core atoms in the side chain of the fifth residue showed an overall accuracy of about k_BT for the four peptides with nonpolar side chains, but twice as large deviations were observed for the peptides with polar side chains. Here, alternative reference‐state Hamiltonians that better cover the conformational space relevant to these peptides are investigated, and post simulation rotational sampling of the χ₁ and χ₂ torsional angles of the fifth residue is carried out to sample different orientations of the side chain. A reference state with rather soft atoms yields accurate ΔG_F,U values for the peptides with the Ser and Thr side chains, but it failed to correctly predict the folding free enthalpy for one peptide with a nonpolar side chain, that is, Leu. Based on the results and those of earlier studies, possible ways to improve the accuracy of the efficient one‐step perturbation technique to compute free enthalpies of folding are discussed. © 2013 Wiley Periodicals, Inc.     

On the choice of a reference state for one‐step perturbation calculations between polar and nonpolar molecules in a polar environment

One‐step perturbation is an efficient method to estimate free energy differences in molecular dynamics (MD) simulations, but its accuracy depends critically on the choice of an appropriate, possibly unphysical, reference state that optimizes the sampling of the physical end states. In particular, the perturbation from a polar moiety to a nonpolar one and vice versa in a polar environment such as water poses a challenge which is of importance when estimating free energy differences that involve entropy changes and the hydrophobic effect. In this work, we systematically study the performance of the one‐step perturbation method in the calculation of the free enthalpy difference between a polar water solute and a nonpolar “water” solute molecule solvated in a box of 999 polar water molecules. Both these polar and nonpolar physical reference states fail to predict the free enthalpy difference as obtained by thermodynamic integration, but the result is worse using the nonpolar physical reference state, because both a properly sized cavity and a favorable orientation of the polar solute in a polar environment are rarely, if ever, sampled in a simulation of the nonpolar solute in such an environment. Use of nonphysical soft‐core reference states helps to sample properly sized cavities, and post‐MD simulation rotational and translational sampling of the solute to be perturbed leads to much improved free enthalpy estimates from one‐step perturbation. © 2012 Wiley Periodicals, Inc.     

Difficulties with multiple time stepping and fast multipole algorithm in molecular dynamics

Numerical experiments are performed on a 36,000-atom protein–DNA–water simulation to ascertain the effectiveness of two devices for reducing the time spent computing long-range electrostatics interactions. It is shown for Verlet-I/r-RESPA multiple time stepping, which is based on approximating long-range forces as widely separated impulses, that a long time step of 5 fs results in a dramatic energy drift and that this is reduced by using an even larger long time step. It is also shown that the use of as many as six terms in a fast multipole algorithm approximation to long-range electrostatics still fails to prevent significant energy drift even though four digits of accuracy is obtained. © 1997 John Wiley & Sons, Inc. J Comput Chem 18 : 1785–1791, 1997     

Screening and identification of multiple constituents and their metabolites of Zhi‐zi‐chi decoction in rat urine and bile by ultra‐high‐performance liquid chromatography quadrupole time‐of‐flight mass spectrometry

Zhi‐zi‐chi decoction (ZZCD) is a classical formula widely used in Chinese clinical application. In the present study, a novel and efficient strategy has been developed for screening and identification of multiple constituents and their metabolites of ZZCD using ultra‐high‐performance liquid chromatography combined with triple time‐of‐flight mass spectrometry. The novel approach of an online data acquisition method dependent on multiple mass defect filter and dynamic background subtraction is combined with multiple data processing techniques. First, a total of 109 potential bioactive compounds were detected in ZZCD. Based on the same instrumental conditions, 100 compounds were found in rat biofluids after oral administration of ZZCD, including 61 original compounds of ZZCD as well as 39 metabolites. Conjugations with sulfate, glucuronate and amino acids were found as the predominant metabolic reaction of ZZCD. As more xenobiotics were detected in urine than those in bile were, it demonstrated that multiple components of ZZCD have undergone comprehensive renal excretion. This study reported the urinary and biliary excretion in rats after oral administration of ZZCD for the first time. The present study expands our knowledge about the constituents and metabolism of ZZCD, which could be very useful for further pharmacological and clinical studies of ZZCD.     

Symplectic algorithm for use in computing the time‐independent Schrödinger equation

The structure of benzoic acid as monomer was studied by semiempirical, ab initio, and density functional methods using several basis sets. The performance of these methods in calculating and describing the vibrational frequencies of benzoic acid and several derivatives was determined. The cyclic dimer form of benzoic acid was also reproduced. Two new procedures of scaling the frequencies were presented. For the ring modes, specific scale equations and scale factors were used from benzene molecule. For the carboxylic group, scaling equations and specific scale factors at different levels were also determined to be used in benzoic acid derivatives. A reassignment of several bands was done. A comparison of the cost/effective method and procedure of scaling was carried out. A significant reduction of the error in the predicted frequencies was obtained over the one‐factor standard scaling procedure. © 2002 Wiley Periodicals, Inc. Int J Quantum Chem, 2002     

Phenylene Ethynylene‐Tethered Perylene Bisimide Folda‐Dimer and Folda‐Trimer: Investigations on Folding Features in Ground and Excited States

In this work, we have elucidated in detail the folding properties of two perylene bisimide (PBI) foldamers composed of two and three PBI units, respectively, attached to a phenylene ethynylene backbone. The folding behaviors of these new PBI folda‐dimer and trimer have been studied by solvent‐dependent UV/Vis absorption and 1D and 2D NMR spectroscopy, revealing facile folding of both systems in tetrahydrofuran (THF). In CHCl₃ the dimer exists in extended (unfolded) conformation, whereas partially folded conformations are observed in the trimer. Temperature‐dependent ¹H NMR spectroscopic studies in [D₈]THF revealed intramolecular dynamic processes for both PBI foldamers due to, on the one hand, hindered rotation around C?N imide bonds and, on the other hand, backbone flapping; the latter process being energetically more demanding as it was observed only at elevated temperature. The structural features of folded conformations of the dimer and trimer have been elucidated by different 2D‐NMR spectroscopy (e.g., ROESY and DOSY) in [D₈]THF. The energetics of folding processes for the PBI dimer and trimer have been assessed by calculations applying various methods, particularly the semiempirical PM6‐DH2 and the more sophisticated B97D approach, in which relevant dispersion corrections are included. These calculations corroborate the results of NMR spectroscopic studies. Folding features in the excited states of these PBI foldamers have been characterized by using time‐resolved fluorescence and transient absorption spectroscopy in THF and CHCl₃, exhibiting similar solvent‐dependent behavior as observed for the ground state. Interestingly, photoinduced electron transfer (PET) process from electron‐donating backbone to electron‐deficient PBI core for extended, but not for folded, conformations was observed, which can be explained by a fast relaxation of excited PBI stacks in the folded conformation into fluorescent excimer states.     

Accurate integration over atomic regions bounded by zero‐flux surfaces

The approach for the integration over a region covered by zero‐flux surface is described. This approach based on the surface triangulation technique is efficiently realized in a newly developed program TWOE . The elaborated method is tested on several atomic properties including the source function. TWOE results are compared with those produced by using well‐known existing programs. Absolute errors in computed atomic properties are shown to range usually from 10^?6 to 10^?5 au. The demonstrative examples prove that present realization has perfect convergence of atomic properties with increasing size of angular grid and allows to obtain highly accurate data even in the most difficult cases. It is believed that the developed program can be bridgehead that allows to implement atomic partitioning of any desired molecular property with high accuracy. © 2012 Wiley Periodicals, Inc.     

Molecular dynamics of cis‐1‐(2‐hydroxy‐5‐ methylphenyl)ethanone oxime and N‐(2‐hydroxy‐4‐methylphenyl)acetamide in solution studied by NMR spectroscopy

The formation of hydrogen bonds and molecular dynamics for the molecules cis‐1‐(2‐hydroxy‐5‐methylphenyl)ethanone oxime ( I ) and N‐(2‐hydroxy‐4‐methylphenyl)acetamide ( II ) have been investigated in solution using NMR. The results confirm the formation of O? H···O, O? H···N and O···H? N type inter‐ and intramolecular hydrogen bonds. Spin‐lattice relaxation times (T₁), activation energy of molecular dynamics and energy of intramolecular hydrogen bonds have been determined. Copyright © 2010 John Wiley & Sons, Ltd.     

Integrating open‐source software applications to build molecular dynamics systems

Three open‐source applications, NanoEngineer‐1, packmol, and mis2lmp are integrated using an open‐source file format to quickly create molecular dynamics (MD) cells for simulation. The three software applications collectively make up the open‐source software (OSS) suite known as MD Studio (MDS). The software is validated through software engineering practices and is verified through simulation of the diglycidyl ether of bisphenol‐a and isophorone diamine (DGEBA/IPD) system. Multiple simulations are run using the MDS software to create MD cells, and the data generated are used to calculate density, bulk modulus, and glass transition temperature of the DGEBA/IPD system. Simulation results compare well with published experimental and numerical results. The MDS software prototype confirms that OSS applications can be analyzed against real‐world research requirements and integrated to create a new capability. © 2014 Wiley Periodicals, Inc.     

Metabolic profiling of quercetin in rats using ultra‐performance liquid chromatography/quadrupole‐time‐of‐flight mass spectrometry

Quercetin, a kind of major flavonoid found in many traditional chinese medicines, is an effective substance for treatments such as lowering blood lipids. However, the studies on quercetin have been mainly focused on its pharmacological effect; the treatment of diseases on a material basis, particularly the metabolites derived from quercetin in vivo , has not been evaluated. In this study, we determined the levels, distributions and types of quercetin's metabolites in plasma, urine, feces and bile of rats after a single oral administration of quercetin at a dose of 80 mg/kg, using ultra‐performance liquid chromatography/quadrupole‐time‐of‐flight mass spectrometry (UPLC‐Q‐TOF/MS). A total of 36 metabolites of quercetin were identified, including 11 metabolites in plasma, 34 metabolites in urine, 12 metabolites in feces and 21 metabolites in bile. The results showed that phase I metabolites were reduction metabolites and phase II metabolites mainly included glucuronidation, sulfation and methylation metabolites. These results provide important information on the metabolism of quercetin, which will be helpful for its further development and utilization.