Fluctuations in molecular dynamics simulations

Statistical fluctuations of a system about its equilibrium state, monitored in a molecular dynamics simulation, are an effective means of computing the thermodynamic and kinetic properties of interfaces in metals and alloys. In this work, three applications of fluctuation analyses are reviewed. First, the capillary fluctuation method can be used to extract the stiffness of grain boundaries, as well as the solid–liquid interfacial free energy and its small anisotropy. Second, both a random walk analysis and a computation of a time correlation function involving the Fourier amplitudes of the interface height can be utilized to derive the mobility of grain boundaries and crystal–melt interfaces. Finally, the probability distribution of premelted grain boundary widths as a function of undercooling can be used to obtain the so-called disjoining potential, which is the thermodynamic driving force responsible for premelting.     

Conformational analysis of compstatin analogues with molecular dynamics simulations in explicit water

The cyclic 13-residue peptide compstatin is a potential therapeutic agent against the unregulated activation of the complement system. A thorough knowledge of its structural and dynamical properties in solution may assist the design of improved complement inhibitors. NMR studies have suggested that the 5-8 segment of free compstatin folds into a critical for activity 5-8 beta turn and the rest of the peptide is mainly disordered. Earlier computational studies of compstatin analogues with a polar-hydrogen/generalized-Born approximation reproduced the 5-8 turn, but also indicated the formation of beta-hairpin or alpha-helical elements and the existence of interactions between certain charged or aromatic sidechains. However, these features are absent or partly present in the NMR spectra, due to extensive conformational averaging. In order to check the compstatin properties with a more rigorous model of the intra- and intermolecular interactions, we conduct here 98-ns all-atom/explicit-water simulations of three compstatin analogues with variable activity; a native analogue, the more active mutant V4W/H9A and the inactive mutant Q5G. The 5-8 beta-turn population is in good accord with NMR. For the systems studied here, the simulations suggest that the 5-8 turn population does not correlate strictly with activity, in agreement with earlier mutational studies. Furthermore, they show structural differences among the analogues outside the 5-8 region. The possible role of these differences in activity is discussed. The probability of beta-hairpin or alpha-helix elements is much smaller with respect to the polar-hydrogen/GB simulations, and the persistent Trp4-Trp7 or Asp6-Arg11 sidechain interactions of the earlier GB studies are not reproduced. The present simulations extend the NMR data and improve our understanding of the properties of compstatin and related analogues.     

Insights into the induced fit mechanism in antithrombin-heparin interaction using molecular dynamics simulations

Heparin was isolated in the beginning of the 20th century and until today remains as one of the most important drugs able to interfere with the haemostatic process. Due to the side effects produced by heparin therapy, new promising drugs have been developed, as the synthetic pentasaccharide (synthetically derived from the sequence GlcN-GlcA-GlcN-IdoA-GlcN). The anticoagulant activity of this compound is based on potentiation of antithrombin (AT) inhibitory activity upon serine proteinases of clotting cascade, a mechanism based on the conformational modification of AT. In this context, we present here a molecular dynamics (MD) study of the interaction between the synthetic pentasaccharide and AT. The obtained data correctly predicted an induced fit mechanism in AT-pentasaccharide interaction, showing a solvent-exposed P1 residue instead of a hided conformation. Also, the specific contribution of important amino acid residues to the overall process was also characterized, both in (2)S(0) and (1)C(4) conformations of IdoA residue, suggesting that there is no conformational requirement to the interaction of this residue with AT. Altogether, the results show that MD simulations could be used to characterize and quantify the interaction of synthetic compounds with AT, predicting its specific capacity to induce conformational changes in AT structure. Thus, MD simulations of heparin (and heparin-derived)-AT interactions are proposed here as a powerful tool to assist and support drug design of new antithrombotic agents.     

硅的平衡态分子动力学模拟

平衡态分子动力学模拟是研究既定系统向所期望的平衡态演化的一种方法,不仅可预测平衡态的各种性质,还为动力学加载过程提供合理的初始条件.本文主要研究Free、NVT、NVE平衡态分子动力学模拟中系统宏观量的演化过程;并讨论如何根据不同的初始条件,选择恰当的平衡态模拟方法.     

Factors responsible for the aggregation behavior of hydrophobic polyelectrolyte PEA in aqueous solution studied by molecular dynamics simulations

Self-association (i.e. interchain aggregation) behavior of atactic poly(ethacrylic acid) PEA in dilute aqueous solution as function of degree-of-neutralization by Na⁺ counter-ions (i.e. charge fraction f) was investigated by molecular dynamics simulations. Aggregation is found to occur in the range 0 ≤ f ≤0.7 in agreement with experimental results compared at specified polymer concentration C_p = 0.36 mol/l in dilute solution. The macromolecular solution was characterized and analysed for radius-of-gyration, torsion angle distribution, inter and intra-molecular hydrogen bonds, radial distribution functions of intermolecular and inter-atomic pairs, inter-chain contacts and solvation enthalpy. The PEA chains form aggregate through attractive inter-chain interaction via hydrogen bonding, in the range f < 0.7, in agreement with experimental observation. The numbers of inter-chain contacts decreases with f. A critical structural transition occurs at f = 0.7, observed via simulations for the first time, in R_g as well as inter-chain H-bonds. The inter-chain distance increases with f due to repulsive interactions between COO− groups on the chains. PEA-PEA electrostatic interactions dominant solvation enthalpy. The PEA solvation enthalpy becomes increasingly favorable with increase in f. The transition enthalpy change, in going from uncharged (acid) state to fully charged state (f = 1) is unfavorable towards aggregate formation.     

Solution conformation of phakellistatin 8 investigated by molecular dynamics simulations

Phakellistatin 8 is a cyclic decapeptide that inhibits cancer cell growth and has sequence and structure similar to antamanide. In molecular dynamics simulations of phakellistatin 8 in water, the decapeptide ring undergoes a conformational change from the saddle-like crystal structure to a more elongated conformation by a transition of the Tyr9 main chain from the alpha L to an extended structure. This is coupled to the loss of the NH9-O6 beta-turn hydrogen bond and the transient dissociation of the Pro7-Tyr9 side-chain packing. Furthermore, the water molecule acting as a transannular bridge forms an additional hydrogen bond with phakellistatin 8, namely with the NH group of Val5 besides those already present in the crystal structure, i.e., with the NH of Ile10 and the CO of Leu6. The alpha-turn hydrogen bond between the Phe4 amide hydrogen and the Ile10 carbonyl oxygen is always present. The solution conformations of the two cyclic decapeptides are similar, in particular in the region involving the NH4-O10 alpha turn of phakellistatin 8 and the NH5-O1 alpha turn of antamanide. The simulation results suggest that in aqueous solution the conformation of phakellistatin 8 is more extended than in the crystalline state, and on a nanosecond time scale phakellistatin 8 is more flexible than antamanide.     

Insights into the folding pathway of the Engrailed Homeodomain protein using replica exchange molecular dynamics simulations

Protein folding studies were carried out by performing microsecond time scale simulations on the ultrafast/fast folding protein Engrailed Homeodomain (EnHD) from Drosophila melanogaster. It is a three-helix bundle protein consisting of 54 residues (PDB ID: 1ENH). The positions of the helices are 8-20 (Helix I), 26-36 (Helix II) and 40-53 (Helix III). The second and third helices together form a Helix-Turn-Helix (HTH) motif which belongs to the family of DNA binding proteins. The molecular dynamics (MD) simulations were performed using replica exchange molecular dynamics (REMD). REMD is a method that involves simulating a protein at different temperatures and performing exchanges at regular time intervals. These exchanges were accepted or rejected based on the Metropolis criterion. REMD was performed using the AMBER FF03 force field with the generalised Born solvation model for the temperature range 286-373 K involving 30 replicas. The extended conformation of the protein was used as the starting structure. A simulation of 600 ns per replica was performed resulting in an overall simulation time of 18 μs. The protein was seen to fold close to the native state with backbone root mean square deviation (RMSD) of 3.16 ?. In this low RMSD structure, the Helix I was partially formed with a backbone RMSD of 3.37 ? while HTH motif had an RMSD of 1.81 ?. Analysis suggests that EnHD folds to its native structure via an intermediate in which the HTH motif is formed. The secondary structure development occurs first followed by tertiary packing. The results were in good agreement with the experimental findings.     

Structure-activity relationships of modified C-terminal endomorphin-2 analogues by molecular dynamics simulations

Motivated by recent experimental works on the modifications of endomorphin-2 (EM2, H-Tyr-Pro-Phe-Phe-NH2) to develop better painkiller, we performed structure-activity-relationship (SAR) studies to investigate modified C-terminal ligands by using molecular dynamics (MD) simulations. Specifically, instead of the CONH2 for the unmodified EM2's C-terminus, the analogue 2 with its C-terminus being CONHNH2 and analogue 3 with its C-terminus being COOMe are studied. First, a systematic conformer search was performed via the quantum chemical AM1 calculations. The cis/trans isomers of the lowest energy were hence selected as MD initial structures. We further showed that EM2s in water exhibited similar dihedral angles to those in DMSO, obtained from the NMR experiment. This similarity indicates the reliability of our MD simulations, and enables us to discuss related bioactivity. Our results showed that the interactions of the Tyr(1)-Phe(3) pair for cis-/trans-EM2s played a considerable role for structural stability. Furthermore, we utilized the chi(1) rotamers of individual aromatic side chains to examine the structural bioactivity. It is shown that this criterion to determine the conformational bioactivity toward mu-opioid receptor (MOR) is insufficient. Thus, we have further employed rotamer-combination approaches to examine the characteristics of SAR for cis-/trans-EM2s. Our results suggested that the bioactive chi(1) rotamers for Tyr(1)-Phe(3) pair remained to favor the [trans-trans] status for MOR selectivity. Therefore, based on the analysis of the chi(1) rotamers, it is suggested that the analogue 2 exhibit greater structural bioactivity for MOR than the analogue 3, and both of them be greater than unmodified EM2 for trans isomers.     

Phonon dispersion measured directly from molecular dynamics simulations

A method to measure the phonon dispersion of a crystal based on molecular dynamics simulation is proposed and implemented as an extension to an open source classical molecular dynamics simulation code LAMMPS. In the proposed method, the dynamical matrix is constructed by observing the displacements of atoms during molecular dynamics simulation, making use of the fluctuation–dissipation theory. The dynamical matrix can then be employed to compute the phonon spectra by evaluating its eigenvalues. It is found that the proposed method is capable of yielding the phonon dispersion accurately, while taking into account the anharmonic effect on phonons simultaneously. The implementation is done in the style of fix of LAMMPS, which is designed to run in parallel and to exploit the functions provided by LAMMPS; the measured dynamical matrices could be passed to an auxiliary postprocessing code to evaluate the phonons.

Program summary

Program title: FixPhonon, version 1.0Catalogue identifier: AEJB_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEJB_v1_0.htmlProgram obtainable from: CPC Program Library, Queen?s University, Belfast, N. IrelandLicensing provisions: GNU General Public licenseNo. of lines in distributed program, including test data, etc.: 105 393No. of bytes in distributed program, including test data, etc.: 3 231 800Distribution format: tar.gzProgramming language: C++Computer: AllOperating system: LinuxHas the code been vectorized or parallelized?: Yes. 1 to N processors may be usedRAM: Depends on problem, ≈1 kB to several MBClassification: 7.8External routines: MPI, FFT, LAMMPS version 15, January 2010 (http://lammps.sandia.gov/)Nature of problem: Atoms in solids make ceaseless vibrations about their equilibrium positions, and a collective vibration forms a wave of allowed wavelength and amplitude. The quantum of such lattice vibration is called the phonon, and the so-called “lattice dynamics” is the field of study to find the normal modes of these vibrations. In other words, lattice dynamics examines the relationship between the frequencies of phonons and the wave vectors, i.e., the phonon dispersion. The evaluation of the phonon dispersion requires the construction of the dynamical matrix. In atomic scale modeling, the dynamical matrices are usually constructed by deriving the derivatives of the force field employed, which cannot account for the effect of temperature on phonons, with an exception of the tedious “quasi-harmonic” procedure.Solution method: We propose here a method to construct the dynamical matrix directly from molecular dynamics simulations, simply by observing the displacements of atoms in the system thus making the constructing of the dynamical matrix a straightforward task. Moreover, the anharmonic effect was taken into account in molecular dynamics simulations naturally, the resultant phonons therefore reflect the finite temperature effect simultaneously.Restrictions: A well defined lattice is necessary to employ the proposed method as well as the implemented code to evaluate the phonon dispersion. In other words, the system under study should be in solid state where atoms vibrate about their equilibrium positions. Besides, no drifting of the lattice is expected. The method is best suited for periodic systems, although non-periodic system with a supercell approach is also possible, it will however become inefficient when the unit cell contains too many atoms.Additional comments: The readers are encouraged to visit http://code.google.com/p/fix-phonon for subsequent update of the code as well as the associated postprocessing code, so as to keep up with the latest version of LAMMPS.Running time: Running time depends on the system size, the numbers of processors used, and the complexity of the force field, like a typical molecular dynamics simulation. For the third example shown in this paper, it took about 2.5 hours on an Intel Xeon X3220 architecture (2.4G, quadcore).References:

• [1] 

  C. Campañá, M.H. Müser, Phys. Rev. B 74 (2006) 075420.

• [2] 

  L.T. Kong, G. Bartels, C. Campañá, C. Denniston, M.H. Müser, Comp. Phys. Commun. 180 (6) (2009) 1004–1010.

     

An expert protein loop refinement protocol by molecular dynamics simulations with restraints

Protein structure prediction (PSP) is a long standing problem in structural biology and bioinformatics. Within the PSP problem loop refinement is a major bottleneck. In this article we report the latest version of the CReF expert predictor system for the PSP problem with emphasis on loop refinement of the approximate 3-D structure 1ZDD_P of the Z34C mini protein predicted by CReF. We designed a loop refinement protocol based on seven molecular dynamics (MD) simulations runs at different temperatures. We found that, by letting the loop residues move freely during dynamics at 325 K and restraining the internal coordinates of the correctly predicted helical structures, while allowing them to move relative to each other, the refinement protocol was very effective in predicting an accurate loop conformation in the first 100 ps of a 1000 ps MD simulation. The quality of the predictions was confirmed by the RMSD between refined and experimental structures which varied from 0.6 to 1.3 Å. In addition, stereochemical analyses showed that 100% of all residues of the refined 1ZDD_P, including those in the loop, populates the most favorable core regions of the Ramachandran plot. Our study suggests that the proposed protocol may be suitable to refine more complex mini proteins with different classes and architectures.     

Reaction analysis and visualization of ReaxFF molecular dynamics simulations

ReaxFF MD (Reactive Force Field Molecular Dynamics) is a promising method for investigating complex chemical reactions in relatively larger scale molecular systems. The existing analysis tools for ReaxFF MD lack the capability of capturing chemical reactions directly by analyzing the simulation trajectory, which is critical in exploring reaction mechanisms. This paper presents the algorithms, implementation strategies, features, and applications of VARxMD, a tool for Visualization and Analysis of Reactive Molecular Dynamics. VARxMD is dedicated to detailed chemical reaction analysis and visualization from the trajectories obtained in ReaxFF MD simulations. The interrelationships among the atoms, bonds, fragments, species and reactions are analyzed directly from the three-dimensional (3D) coordinates and bond orders of the atoms in a trajectory, which are accomplished by determination of atomic connectivity for recognizing connected molecular fragments, perception of bond types in the connected fragments for molecules or radicals, indexing of all these molecules or radicals (chemical species) based on their 3D coordinates and recognition of bond breaking or forming in the chemical species for reactions. Consequently, detailed chemical reactions taking place between two sampled frames can be generated automatically. VARxMD is the first tool specialized for reaction analysis and visualization in ReaxFF MD simulations. Applications of VARxMD in ReaxFF MD simulations of coal and HDPE (high-density polyethylene) pyrolysis show that VARxMD provides the capabilities in exploring the reaction mechanism in large systems with complex chemical reactions involved that are difficult to access manually.     

Stability of C60 chains: molecular dynamics simulations

A linearly aligned structure of three C60 fullerene, interconnected by two benzorods of same size, have been investigated under heat treatment. The overall structure resembles a section of a beaded string. Nine different lengths of benzorods have been considered, and the effect on the thermal stability have been investigated by means of molecular dynamics method. It has been found that the structure is thermally stable up to elevated temperatures, and the linear alignment of the structure is persistent, up to the temperature of decomposition.     

Structural analysis of the inhibition of APRIL by TACI and BCMA through molecular dynamics simulations

APRIL (a proliferation-inducing ligand) is a member of the tumour necrosis factor (TNF) superfamily that binds the receptors (TNFRs) TACI and BCMA. Since it was discovered, a great amount of evidence has been reported about the involvement of APRIL in autoimmune diseases including systemic lupus erythematosus (SLE), rheumatoid arthritis (RA), Sjögren's syndrome (SS) and multiple sclerosis (MS). In addition, an important role of APRIL has been described in different types of tumour cell lines and in a variety of primary tumour tissues where, in contrast with the normal ones, high mRNA levels have been detected. Accordingly, the design of compounds mimicking the inhibition of APRIL by its receptors appears to be a promising way to treat autoimmune and cancer diseases. As a first step to achieve these goals and in order to better understand the key interactions involved in these systems, we report a structural analysis of the inhibition of human and murine APRIL by its human receptors TACI and BCMA obtained by molecular dynamics simulations. Although most of the key interactions can be obtained from the existing experimental information, new described interactions between human APRIL and its receptors can contribute to a better design of APRIL inhibitors.     

Molecular basis of agonicity and antagonicity in the androgen receptor studied by molecular dynamics simulations

Treatment of prostate cancer patients with antiandrogens is initially successful, though the therapy often becomes refractory over the time. This mechanism is not fully understood, but the presence of androgen receptor (AR) mutant forms which are activated by antiandrogens and other endogenous ligands, and overexpression of the receptor have been suggested. In an attempt to explain the molecular basis for agonicity and antagonicity in the androgen receptor, and the changes on biological activity of subtle modifications at the ligand and receptor (mutations) level, molecular dynamics simulations were performed on the androgen receptor wild type (WT), and T877A and W741 mutant forms, complexed with several non-steroidal androgens. The stabilizing role of residues from helices 3, 5, 11 and 12 was observed in non-steroidal androgens R-3, S-1, and R-bicalutamide and hydroxyflutamide in resistant mutations. In the AR WT antiandrogen R-bicalutamide complex, destabilization of M895 by both W741 and the sulfonyl linkage of the ligand may be responsible for reported antagonism. Changes in the ligand or mutations alleviating this effect were observed to stabilize the receptor in the active conformation, thus developing resistance to R-bicalutamide. The results presented provide a plausible explanation for the molecular basis of agonicity and antagonicity in the androgen receptor, and complement previous studies using static crystal structures, incorporating for the first time protein dynamics into the analysis. Thus, our results provide a valuable framework for the structure-based design of improved antiandrogens.     

Visualizations of molecular dynamics simulations of high-performance polycrystalline structural ceramics

Initiated by the Department of Defense (DOD) High Performance Computing Modernization Program (HPCMP), the Data Analysis and Assessment Center (DAAC), serves the needs of DOD HPCMP scientists by facilitating the analysis of an ever-increasing volume and complexity of data [1]. A research scientist and HPCMP user ran nanoscale molecular dynamics simulations using Large-scale Atomic/Molecular Massively Parallel Simulator code (LAMMPS) from Sandia National Labs. The largest simulation contained over 70 million atoms (Fig. 7). Data sets this large are required to study crack propagation and failure mechanisms that span multiple length scales with atomic resolution. The DAAC developed new methods to visualize the time evolution of data sets this large. The size and complexity of the molecular dynamics simulations and the analytics required the use of DOD HPCMP High Performance Computing (HPC) resources.     

Hybrid parallelization of molecular dynamics simulations to reduce load imbalance

The most widely used technique to allow for parallel simulations in molecular dynamics is spatial domain decomposition, where the physical geometry is divided into boxes, one per processor. This technique can inherently produce computational load imbalance when either the spatial distribution of particles or the computational cost per particle is not uniform. This paper shows the benefits of using a hybrid MPI+OpenMP model to deal with this load imbalance. We consider LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator), a prototypical molecular dynamics simulator that provides its own balancing mechanism and an OpenMP implementation for many of its modules, allowing for a hybrid setup. In this work, we extend the current OpenMP implementation of LAMMPS and optimize it and evaluate three different setups: MPI-only, MPI with the LAMMPS balance mechanism, and hybrid setup using our improved OpenMP version. This comparison is made using the five standard benchmarks included in the LAMMPS distribution plus two additional test cases. Results show that the hybrid approach can deal with load balancing problems better and more effectively (50% improvement versus MPI-only for a highly imbalanced test case) than the LAMMPS balance mechanism (only 43% improvement) and improve simulations with issues other than load imbalance.

     

MDWiZ: A platform for the automated translation of molecular dynamics simulations

A variety of popular molecular dynamics (MD) simulation packages were independently developed in the last decades to reach diverse scientific goals. However, such non-coordinated development of software, force fields, and analysis tools for molecular simulations gave rise to an array of software formats and arbitrary conventions for routine preparation and analysis of simulation input and output data. Different formats and/or parameter definitions are used at each stage of the modeling process despite largely contain redundant information between alternative software tools. Such Babel of languages that cannot be easily and univocally translated one into another poses one of the major technical obstacles to the preparation, translation, and comparison of molecular simulation data that users face on a daily basis. Here, we present the MDWiZ platform, a freely accessed online portal designed to aid the fast and reliable preparation and conversion of file formats that allows researchers to reproduce or generate data from MD simulations using different setups, including force fields and models with different underlying potential forms. The general structure of MDWiZ is presented, the features of version 1.0 are detailed, and an extensive validation based on GROMACS to LAMMPS conversion is presented. We believe that MDWiZ will be largely useful to the molecular dynamics community. Such fast format and force field exchange for a given system allows tailoring the chosen system to a given computer platform and/or taking advantage of a specific capabilities offered by different software engines.     

Denoising methods for thermomechanical decomposition for quasi-equilibrium molecular dynamics simulations

In this work, we applied robust denoising methods well established in the signal processing field for the thermomechanical decomposition of velocity data obtained from molecular dynamics (MD) simulations. In the decomposition, the atomic velocity was assumed to be the sum of the mechanical velocity and the thermal velocity, which can be linked to the stress and temperature field at the continuum scale, respectively. For the quasi-equilibrium process, with the thermal velocity treated as the Gaussian distributed stationary noise with zero mean and a positive variance that is linearly proportional to temperature and mechanical velocity as the clean signal, the velocity decomposition can be recasted into a denoising problem, for which powerful denoising methods have been developed to estimate a clean signal from noisy data. We investigated the widely-used linear parametric real-domain, linear nonparametric Fourier-based, and nonlinear nonparametric wavelet-based denoising methods, first on their theoretical properties and then made comparsion among them for denosing some synthetic noisy 1-dimensional (1D) data generated from MD simulations. The nonlinear wavelet-based thresholding estimator possessed better optimality properties than the other estimators, and also outperformed the other estimators in the synthetic data test. A further test comparing the various denoising methods for an adiabatic shear crack nucleation and propagation process simulated using MD simulations showed better performance by the wavelet-based denoisng method. Results from this work reveal good potential of applying wavelet-based denoising method to the study of thermomechanical processes simulated using MD simulations.     

Accelerating molecular dynamics simulations using Graphics Processing Units with CUDA

Molecular dynamics is an important computational tool to simulate and understand biochemical processes at the atomic level. However, accurate simulation of processes such as protein folding requires a large number of both atoms and time steps. This in turn leads to huge runtime requirements. Hence, finding fast solutions is of highest importance to research. In this paper we present a new approach to accelerate molecular dynamics simulations with inexpensive commodity graphics hardware. To derive an efficient mapping onto this type of computer architecture, we have used the new Compute Unified Device Architecture programming interface to implement a new parallel algorithm. Our experimental results show that the graphics card based approach allows speedups of up to factor nineteen compared to the corresponding sequential implementation.