Harvesting graphics power for MD simulations

We discuss an implementation of molecular dynamics (MD) simulations on a graphic processing unit (GPU) in the NVIDIA CUDA language. We tested our code on a modern GPU, the NVIDIA GeForce 8800 GTX. Results for two MD algorithms suitable for short-ranged and long-ranged interactions, and a congruential shift random number generator are presented. The performance of the GPU's is compared to their main processor counterpart. We achieve speedups of up to 40, 80 and 150 fold, respectively. With the latest generation of GPU's one can run standard MD simulations at 10⁷ flops/$.     

Scrutiny of the mechanism of small molecule inhibitor preventing conformational transition of amyloid-β42 monomer: insights from molecular dynamics simulations

Alzheimer’s disease (AD) is a progressive neurodegenerative disorder that is characterized by loss of intellectual functioning of brain and memory loss. According to amyloid cascade hypothesis, aggregation of amyloid-β₄₂ (Aβ₄₂) peptide can generate toxic oligomers and their accumulation in the brain is responsible for the onset of AD. In spite of carrying out a large number of experimental studies on inhibition of Aβ₄₂ aggregation by small molecules, the detailed inhibitory mechanism remains elusive. In the present study, comparable molecular dynamics (MD) simulations were performed to elucidate the inhibitory mechanism of a sulfonamide inhibitor C1 (2,5-dichloro-N-(4-piperidinophenyl)-3-thiophenesulfonamide), reported for its in vitro and in vivo anti-aggregation activity against Aβ₄₂. MD simulations reveal that C1 stabilizes native α-helix conformation of Aβ₄₂ by interacting with key residues in the central helix region (13–26) with hydrogen bonds and π–π interactions. C1 lowers the solvent-accessible surface area of the central hydrophobic core (CHC), KLVFF (16–20), that confirms burial of hydrophobic residues leading to the dominance of helical conformation in the CHC region. The binding free energy analysis with MM–PBSA demonstrates that Ala2, Phe4, Tyr10, Gln15, Lys16, Leu17, Val18, Phe19, Phe20, Glu22, and Met35 contribute maximum to binding free energy (?43.1 kcal/mol) between C1 and Aβ₄₂ monomer. Overall, MD simulations reveal that C1 inhibits Aβ₄₂ aggregation by stabilizing native helical conformation and inhibiting the formation of aggregation-prone β-sheet conformation. The present results will shed light on the underlying inhibitory mechanism of small molecules that show potential in vitro anti-aggregation activity against Aβ₄₂.     

Systolic Loop Methods for Molecular Dynamics Simulation,Generalised for Macromolecules

Abstract

Systolic loop programs have been shown to be very efficient for molecular dynamics simulations of liquid systems on distributed memory parallel computers. The original methods address the case where the number of molecules simulated exceeds the number of processors used. Simulations of large flexible molecules often do not meet this condition, requiring the three- and four-body terms used to model chemical bonds within a molecule to be distributed over several processors. This paper discusses how the systolic loop methods may be generalised to accommodate such systems, and describes the implementation of a computer program for simulation of protein dynamics. Performance figures are given for this program running typical simulations on a Meiko Computing Surface using different number of processors.     

Influence of dielectric constant on protonation equilibria of maleic acid and L-asparagine in acetonitrile water-mixtures

Abstract

The protonation constants of maleic acid and L-asparagine have been studied pH-metrically in various concentrations (0–50% v/v) of acetonitrile–water mixtures maintaining an ionic strength of 0.16 mol L^-1 at 30⁰C. The protonation constants have been calculated using the computer program MINIQUAD75 and are selected based on statistical parameters. Linear variation of step-wise protonation constants (log K) with the reciprocal of the dielectric constant of the solvent mixture has been attributed to the dominance of the electrostatic forces.     

Comparative study of stability and transport of molecules through cyclic peptide nanotube and aquaporin: a molecular dynamics simulation approach

Abstract

The structural stability and transport properties of the cyclic peptide nanotube (CPN) 8?×?[Cys–Gly–Met–Gly]₂ in different phospholipid bilayers such as POPA (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphatidic acid), POPE (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoethanolamine), POPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine), POPG (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoglycerol) and POPS (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoserine) with water have been investigated using molecular dynamics (MD) simulation. The hydrogen bonds and non-bonded interaction energies were calculated to study the stability in different bilayers. One µs MD simulation in POPA lipid membrane reveals the stability of the cyclic peptide nanotube, and the simulations at various temperatures manifest the higher stability of 8?×?[Cys–Gly–Met–Gly]₂. We demonstrated that the presence of sulphur-containing amino acids in CPN enhances the stability through disulphide bonds between the adjacent rings. Further, the water permeation coefficient of the CPN is calculated and compared with human aquaporin-2 (AQP2) channel protein. It is found that the coefficients are highly comparable to the AQP2 channel though the mechanism of water transport is not similar to AQP 2; the flow of water in the CPN is taking place as a two-line 1–2–1–2 file fashion. In addition to that, the transport behavior of Na⁺ and K⁺ ions, single water molecule, urea and anti-cancer drug fluorouracil were investigated using pulling simulation and potential of mean force calculation. The above transport behavior shows that Na⁺ is trapped in CPN for a longer time than other molecules. Also, the interactions of the ions and molecules in Cα and mid-Cα plane were studied to understand the transport behavior of the CPN. Abbreviations AQP2 Aquaporin-2

CPN Cyclic peptide nanotube

MD Molecular dynamics

POPA 1-Palmitoyl-2-oleoyl-sn-glycero-3-phosphatidic acid

POPE 1-Palmitoyl-2-oleoyl-sn-glycero-3-phosphoethanolamine

POPG 1-Palmitoyl-2-oleoyl-sn-glycero-3-phosphoglycerol

POPS 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoserine

Communicated by Ramaswamy H. Sarma     

FF12MC: A revised AMBER forcefield and new protein simulation protocol

Specialized to simulate proteins in molecular dynamics (MD) simulations with explicit solvation, FF12MC is a combination of a new protein simulation protocol employing uniformly reduced atomic masses by tenfold and a revised AMBER forcefield FF99 with (i) shortened C? H bonds, (ii) removal of torsions involving a nonperipheral sp³ atom, and (iii) reduced 1–4 interaction scaling factors of torsions ? and ψ. This article reports that in multiple, distinct, independent, unrestricted, unbiased, isobaric–isothermal, and classical MD simulations FF12MC can (i) simulate the experimentally observed flipping between left‐ and right‐handed configurations for C14–C38 of BPTI in solution, (ii) autonomously fold chignolin, CLN025, and Trp‐cage with folding times that agree with the experimental values, (iii) simulate subsequent unfolding and refolding of these miniproteins, and (iv) achieve a robust Z score of 1.33 for refining protein models TMR01, TMR04, and TMR07. By comparison, the latest general‐purpose AMBER forcefield FF14SB locks the C14–C38 bond to the right‐handed configuration in solution under the same protein simulation conditions. Statistical survival analysis shows that FF12MC folds chignolin and CLN025 in isobaric–isothermal MD simulations 2–4 times faster than FF14SB under the same protein simulation conditions. These results suggest that FF12MC may be used for protein simulations to study kinetics and thermodynamics of miniprotein folding as well as protein structure and dynamics. Proteins 2016; 84:1490–1516. © 2016 The Authors Proteins: Structure, Function, and Bioinformatics Published by Wiley Periodicals, Inc.     

88 Comparison of monovalent and divalent ion distributions around a DNA duplex with molecular dynamic simulation and Poisson–Boltzmann approach

Ion interactions with nucleic acids (both DNA and RNA) are an important and evolving field of investigation. Positively charged cations may interact with highly negatively charged nucleic acids via simple electrostatic interactions to help screen the electrostatic repulsion along the nucleic acids and assist their folding and/or compaction. Cations may also bind at specific sites and become integral parts of the structures, possibly playing important enzymatic roles. Two popular methods for computationally exploring a nucleic acid’s ion atmosphere are atomistic molecular dynamics (MD) simulations and the Poisson–Boltzmann (PB) equation. In general, monovalent ion results obtained from MD simulations and the PB equation agree well with experiment. However, Bai et al. (2007) observed discrepancies between experiment and the PB equation while examining the competitive binding of monovalent and divalent ions, with more significant discrepancies for divalent ions. The goal of this project was to thoroughly investigate monovalent (Na⁺) and divalent (Mg²⁺) ion distributions formed around a DNA duplex with MD simulations and the PB equation. We simulated three different cation concentrations, and matched the equilibrated bulk ion concentration for our theoretical calculations with the PB equation. Based on previous work, our Mg²⁺ ions were fully solvated, the expected state of Mg²⁺ ions when interacting with a duplex, when the production simulations began and remained throughout the simulations (Kirmizialtin, 2010; Robbins, 2012). Na⁺ ion distributions and number of Na⁺ ions within 10?Å of the DNA obtained from our two methods agreed well. However, results differed for Mg²⁺ ions, with a lower number of ions within the cut-off distance obtained from the PB equation when compared to MD simulations. The Mg²⁺ ion distributions around the DNA obtained via the two methods also differed. Based on our results, we conclude that the PB equation will systematically underestimate Mg²⁺ ions bound to DNA, and much of this deviation is attributed to dielectric saturation associated with high valency ions.     

Improved model of hydrated calcium ion for molecular dynamics simulations using classical biomolecular force fields

Calcium ions (Ca²⁺) play key roles in various fundamental biological processes such as cell signaling and brain function. Molecular dynamics (MD) simulations have been used to study such interactions, however, the accuracy of the Ca²⁺ models provided by the standard MD force fields has not been rigorously tested. Here, we assess the performance of the Ca²⁺ models from the most popular classical force fields AMBER and CHARMM by computing the osmotic pressure of model compounds and the free energy of DNA–DNA interactions. In the simulations performed using the two standard models, Ca²⁺ ions are seen to form artificial clusters with chloride, acetate, and phosphate species; the osmotic pressure of CaAc₂ and CaCl₂ solutions is a small fraction of the experimental values for both force fields. Using the standard parameterization of Ca²⁺ ions in the simulations of Ca²⁺‐mediated DNA–DNA interactions leads to qualitatively wrong outcomes: both AMBER and CHARMM simulations suggest strong inter‐DNA attraction whereas, in experiment, DNA molecules repel one another. The artificial attraction of Ca²⁺ to DNA phosphate is strong enough to affect the direction of the electric field‐driven translocation of DNA through a solid‐state nanopore. To address these shortcomings of the standard Ca²⁺ model, we introduce a custom model of a hydrated Ca²⁺ ion and show that using our model brings the results of the above MD simulations in quantitative agreement with experiment. Our improved model of Ca²⁺ can be readily applied to MD simulations of various biomolecular systems, including nucleic acids, proteins and lipid bilayer membranes. © 2016 Wiley Periodicals, Inc. Biopolymers 105: 752–763, 2016.     

Probing the interaction between cHAVc3 peptide and the EC1 domain of E-cadherin using NMR and molecular dynamics simulations

The goal of this work is to probe the interaction between cyclic cHAVc3 peptide and the EC1 domain of human E-cadherin protein. Cyclic cHAVc3 peptide (cyclo(1,6)Ac-CSHAVC-NH₂) binds to the EC1 domain as shown by chemical shift perturbations in the 2D ¹H,-¹⁵N-HSQC NMR spectrum. The molecular dynamics (MD) simulations of the EC1 domain showed folding of the C-terminal tail region into the main head region of the EC1 domain. For cHAVc3 peptide, replica exchange molecular dynamics (REMD) simulations generated five structural clusters of cHAVc3 peptide. Representative structures of cHAVc3 and the EC1 structure from MD simulations were used in molecular docking experiments with NMR constraints to determine the binding site of the peptide on EC1. The results suggest that cHAVc3 binds to EC1 around residues Y36, S37, I38, I53, F77, S78, H79, and I94. The dissociation constants (K_d values) of cHAVc3 peptide to EC1 were estimated using the NMR chemical shifts data and the estimated K_ds are in the range of .5 × 10^?5–7.0 × 10^?5 M.     

Atomistic Simulations of Rare Gas Transport through Breathable Single-wall Nanotubes with Constrictions and Knees

We present the results of molecular dynamics (MD) computer simulations of rare gas diffusion through breathable nanotubes with pentagon–heptagon pair defects resulting in constrictions and knees. Diffusion involves interrupted high speed “choppy” motion with intermittent reversal in velocity direction. Single atoms exhibit a spiral-like path, in contrast to atoms traveling in groups. Considerable resistance to flow appears to reside in the upstream section of the nanotube where density gradients are small, prior to the constriction. Subsequently, considerable density gradients are present and speeds increase, becoming greatest at the tube exit. For the nanotubes examined, Kr and Xe diffusion was too hindered to provide reliable results. Diffusion of He through the nanotubes with knees occurs in a single-file fashion nearly along the center of the tube and the knee has no detectable effect on the diffusion kinetics. Transport diffusion coefficients are in the order of 10^-4–10^-2?cm²/s.     

The new program OPAL for molecular dynamics simulations and energy refinements of biological macromolecules

Summary A new program for molecular dynamics (MD) simulation and energy refinement of biological macromolecules, OPAL, is introduced. Combined with the supporting program TRAJEC for the analysis of MD trajectories, OPAL affords high efficiency and flexibility for work with diferent force fields, and offers a user-friendly interface and extensive trajectory analysis capabilities. Salient features are computational speeds of up to 1.5 GFlops on vector supercomputers such as the NEC SX-3, ellipsoidal boundaries to reduce the system size for studies in explicit solvents, and natural treatment of the hydrostatic pressure. Practical applications of OPAL are illustrated with MD simulations of pure water, energy minimization of the NMR structure of the mixed disulfide of a mutant E. coli glutaredoxin with glutathione in different solvent models, and MD simulations of a small protein, pheromone Er-2, using either instantaneous or time-averaged NMR restraints, or no restraints.Abbreviations D diffusion constant in cm²/s - Er-2 pheromone 2 from Euplotes raikovi - GFlop one billion floating point operations per second - Grx(C14S)-SG mixed disulfide between a mutant E. coli glutaredoxin, with Cys¹⁴ replaced by Ser, and glutathione - MD molecular dynamics - NOE nuclear Overhauser enhancement - rmsd root-mean-square deviation - density in g/cm³     

Global fold of human cannabinoid type 2 receptor probed by solid‐state 13C‐, 15N‐MAS NMR and molecular dynamics simulations

The global fold of human cannabinoid type 2 (CB₂) receptor in the agonist‐bound active state in lipid bilayers was investigated by solid‐state ¹³C‐ and ¹⁵N magic‐angle spinning (MAS) NMR, in combination with chemical‐shift prediction from a structural model of the receptor obtained by microsecond‐long molecular dynamics (MD) simulations. Uniformly ¹³C‐ and ¹⁵N‐labeled CB₂ receptor was expressed in milligram quantities by bacterial fermentation, purified, and functionally reconstituted into liposomes. ¹³C MAS NMR spectra were recorded without sensitivity enhancement for direct comparison of C_α, C_β, and C?O bands of superimposed resonances with predictions from protein structures generated by MD. The experimental NMR spectra matched the calculated spectra reasonably well indicating agreement of the global fold of the protein between experiment and simulations. In particular, the ¹³C chemical shift distribution of C_α resonances was shown to be very sensitive to both the primary amino acid sequence and the secondary structure of CB₂. Thus the shape of the C_α band can be used as an indicator of CB₂ global fold. The prediction from MD simulations indicated that upon receptor activation a rather limited number of amino acid residues, mainly located in the extracellular Loop 2 and the second half of intracellular Loop 3, change their chemical shifts significantly (≥1.5 ppm for carbons and ≥5.0 ppm for nitrogens). Simulated two‐dimensional ¹³C_α(i)? ¹³C?O(i) and ¹³C?O(i)? ¹⁵NH(i + 1) dipolar‐interaction correlation spectra provide guidance for selective amino acid labeling and signal assignment schemes to study the molecular mechanism of activation of CB₂ by solid‐state MAS NMR. Proteins 2014; 82:452–465. © 2013 Wiley Periodicals, Inc.     

Molecular dynamics simulations of palmitoyloleoylphosphatidylglycerol bilayers

Phosphatidylglycerol (PG) is one of the important components of biological membranes, but there is a paucity of experimental data to test the accuracy of molecular dynamics (MD) simulations. This work consists of testing the accuracy of the CHARMM36 (C36) lipid force field on 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoglycerol (POPG) lipid bilayers. MD simulations of POPG lipid bilayers are compared to recently available X-ray and neutron scattering and deuterium NMR measurements. Overall, the C36 lipid force field accurately represents the X-ray and neutron form factors, bilayer and hydrocarbon thicknesses and chain deuterium order parameters. The surface area per lipid from MD simulations with C36 (67.7 ± 0.2 Å²) is in excellent agreement with the experimentally determined value of 66.0 ± 1.3 Å². C36 outperforms the lipid force field developed by Berger et al. [15] and suggests that past studies with this force field may result in lateral areas that are too small. Moreover, our studies give some insight into the structural model used in experiments and suggest that the functional form for the head group may not be Gaussian-like. Based on our simulations, the POPG lipid in the C36 lipid force field is well parameterised and can be used for other PG lipids and membrane models with mixed lipids.     

Is soil carbon disappearing? The dynamics of soil organic carbon in Java

Research in the soil of the tropics mostly has demonstrated the decline of soil organic carbon (SOC) after conversion of primary forest to plantation and cultivated lands. This paper illustrates the dynamics of SOC on the island of Java, Indonesia, from 1930 to 2010. We used 2002 soil profile observations containing organic carbon (C) analysis in the topsoil, which were collected by the Indonesian Center for Agricultural Land Resources Research & Development from 1923 to 2007. Results show the obvious decline of SOC values from around 2% in 1930–1940 to 0.8% in 1960–1970. However, there has been an increase of SOC content since 1970, with a median level of 1.1% in the year 2000. Our analysis suggests that the human influence and agricultural practices on SOC in Java have been a stronger influence than the environmental factors. SOC for the top 10 cm has shown a net accumulation rate of 0.2–0.3 Mg C ha^?1 yr^?1 during the period 1990–2000. These findings give rise to optimism for increased soil C sequestration in the tropics.     

Molecular Dynamics Simulations of Complexes Between Wild-Type and Mutant Anthrax Protective Antigen Variants and a Model Anthrax Toxin Receptor

Abstract

Bacillus anthracis, a spore-forming infectious bacterium, produces a toxin consisting of three proteins: lethal factor (LF), edema factor (EF), and protective antigen (PA). LF and EF possess intracellular enzymatic functions, the net effect of which is to severely compromise host innate immunity. During an anthrax infection PA plays the critical role of facilitating entry of both EF and LF toxins into host cell cytoplasm. Crystal structures of all three of the anthrax toxins have been determined, as well as the crystal structure of the (human) von Willebrand factor A (integrin VWA/I domain)—an anthrax toxin receptor. A theoretical structure of the complex between VWA/I and PA has also been reported. Here we report on the results of 1,000 psec molecular dynamics (MD) simulations carried out on complexes between the Anthrax Protective Antigen Domain 4 (PA-D4) and the von Willebrand Factor A (VWA/I). MD simulations (using Insight II software) were carried out for complexes containing wildtype (WT) PA-D4, as well as for complexes containing three different mutants of PA-D4, one containing three substitutions in the PA-D4 “small loop” (residues 679–693) (D683A/L685E/Y688C), one containing a single substitution at a key site at the PA-D4—receptor interface (K679A) and another containing a deletion of eleven residues at the C-terminus of PA (A724–735). All three sets of PA mutations have been shown experimentally to result in serious deficiencies in PA function. Our MD results are consistent with these findings. Major disruptions in interactions were observed between the mutant PA-D4 domains and the anthrax receptor during the MD simulations. Many secondary structural features in PA-D4 are also severely compromised when VWA complexes with mutant variants of PA-D4 are subjected to MD simulations. These MD simulation results clearly indicate the importance of the mutated PA-D4 residues in both the “small loop” and at the carboxyl terminus in maintaining a PA conformation that is capable of effective interaction with the anthrax toxin receptor.     

Accelerating multiple replica molecular dynamics simulations using the Intel® Xeon Phi™ coprocessor

We investigate the performance increase provided by the Intel^® Xeon Phi? coprocessor in multiple replica molecular dynamics applications using a novel parallelisation scheme. The benefits of the proposed parallelisation scheme are demonstrated by glycine in water, a system of significant interest in the crystallisation simulation community. The molecular dynamics (MD) engine consists of initially serial LAMMPS and NAMD subroutines, and is subsequently modified and parallelised using a heterogeneous programming model, where each MPI rank is paired with a unique Intel^® Xeon Phi? coprocessor and CPU socket. The MD engine is parallelised using an OpenMP atom domain decomposition algorithm on the Intel^® Xeon Phi? coprocessor and OpenMP task parallelism on the host CPU socket. Using nodes with two Intel^® Xeon Phi? coprocessors, one per socket, we demonstrate that a factor of five reduction in the required computational resources is achieved per replica with the coprocessor, when compared against employing the standard spatial domain decomposition algorithm with no accelerator. Furthermore, the proposed parallelisation scheme achieves ideal weak scaling with respect to the number of employed MPI ranks (replicas). The Intel^® Xeon Phi? coprocessor not only allows us to the increase performance output per socket by a factor of five, when compared against no accelerators, but also significantly reduces the parallelisation complexity necessary to achieve this performance, as the Intel^® Xeon Phi? coprocessor operates using the simple OpenMP programming language.     

"SCINT," a Fortran IV program for the analysis of data for doubly labeled and variably quenched liquid scintillation samples

For some time now, we have been using double label techniques in our laboratory. In our work we generate large numbers of samples containing ¹⁴C and ³H and varying amounts of quenching agents. The analysis of this data led us to develop SCINT, a simple and versatile program written in FORTRAN IV for use on our Control Data 1700 computer facility. We feel that such a program may be of use to the many workers now involved with dual isotope procedures who also have access to similar computer facilities.     

Solid‐state NMR and simulation studies of equinatoxin II N‐terminus interaction with lipid bilayers

The interaction with model membranes of a peptide, EqtII_1–32, corresponding to the N‐terminal region of the pore‐forming toxin equinatoxin II (EqtII) has been studied using solid‐state NMR and molecular dynamics (MD) simulations. The distances between specifically labeled nuclei in [¹⁹F‐para]Phe16‐[1‐¹³C]Leu19 and [¹⁹F‐para]Phe16‐[¹⁵N]Leu23 analogs of EqtII_1–32 measured by REDOR in lyophilized peptide were in agreement with published crystal and solution structures. However, in both DMPC and mixed DMPC:SM membrane environments, significant changes in the distances between the labeled amino acid pairs were observed, suggesting changes in helical content around the experimentally studied region, 16–23, in the presence of bilayers. ¹⁹F‐³¹P REDOR experiments indicated that the aromatic ring of Phe16 is in contact with lipid headgroups in both membrane environments. For the DMPC:SM mixed bilayers, a closer interaction between Phe16 side chains and lipid headgroups was observed, but an increase in distances was observed for both labeled amino acid pairs compared with those measured for EqtII_1–32 in pure DMPC bilayers. The observed differences between DMPC and DMPC:SM bilayers may be due to the greater affinity of EqtII for the latter. MD simulations of EqtII_1–32 in water, on a pure DMPC bilayer, and on a mixed DMPC:SM bilayer indicate significant peptide secondary structural differences in the different environments, with the DMPC‐bound peptide adopting helical formations at residues 16–24, whereas the DMPC:SM‐bound peptide exhibits a longer helical stretch, which may contribute to its enhanced activity against PC:SM compared with pure PC bilayers. Proteins 2010. © 2009 Wiley‐Liss, Inc.