Molecular Dynamics Simulations of Ibuprofen Binding to Aβ Peptides

Using replica exchange molecular dynamics simulations and the implicit solvent model we probed binding of ibuprofen to Aβ_10-40 monomers and amyloid fibrils. We found that the concave (CV) fibril edge has significantly higher binding affinity for ibuprofen than the convex edge. Furthermore, binding of ibuprofen to Aβ monomers, as compared to fibrils, results in a smaller free energy gain. The difference in binding free energies is likely to be related to the presence of the groove on the CV fibril edge, in which ibuprofen tends to accumulate. The confinement effect of the groove promotes the formation of large low-energy ibuprofen clusters, which rarely occur on the surface of Aβ monomers. These observations led us to suggest that the ibuprofen binding mechanism for Aβ fibrils is different from that for monomers. In general, ibuprofen shows a preference to bind to those regions of Aβ monomers (amino terminal) and fibrils (the CV edge) that are also the primary aggregation interfaces. Based on our findings and on available experimental data, we propose a rationale for the ibuprofen antiaggregation effect.     

Effect of Cation Concentration on Molecular Dynamics Simulations of UDP-Glucose

Abstract

Glycosyl esters of nucleoside di or mono-phosphates, generally referred to as “sugar nucleotides”, serve as a sugar donors during the biosynthesis of oligo- and polysaccharides; they are therefore of a primary importance in carbohydrate metabolism in the living world. Molecular dynamics simulations were used to explore the conformational flexibility of one nucleotide sugar, UDP-glucose (UDP-Glc). The AMBER program package was used with some new parameters especially developed for nucleotide sugars. Several simulations on this molecule in aqueous solution, each of 2 ns duration, were carried out for increasing concentrations of monovalent K⁺ and divalent Mg²⁺ ions. For the monovalent ion, it is revealed that its presence and concentration is crucial for the conformational behavior, resulting in the stabilization of the extended conformation. The preferred location of K⁺ is in close proximity to the negatively charged phosphate oxygens, but the ion moves freely and can occupy other sites. Since the size of this cation is close that of the water molecules, the hydration scheme is not perturbed. Completely different results are obtained when the divalent Mg²⁺ cation is introduced in the simulation. A very strong interaction is established between the phosphate group and the cation; as a result the UDP-Glc molecule is locked in a rigid extended geometry. The analyses of the trajectories provide new insight on the role of the metal ion in the catalytic mechanism of glycosyltransferases.     

Molecular Dynamics Simulations of Phospholipid Bilayers

Abstract

Molecular dynamics (MD) simulations at 37°C have been performed on three phospholipid bilayer systems composed of the lipids DLPE, DOPE, and DOPC. The model used included 24 explicit lipid molecules and explicit waters of solvation in the polar head group regions, together with constant-pressure periodic boundary conditions in three dimensions. Using this model, a MD simulation samples part of an infinite planar lipid bilayer. The lipid dynamics and packing behavior were characterized. Furthermore, using the results of the simulations, a number of diverse properties including bilayer structural parameters, hydrocarbon chain order parameters, dihedral conformations, electron density profile, hydration per lipid, and water distribution along the bilayer normal were calculated. Many of these properties are available for the three lipid systems chosen, making them well suited for evaluating the model and protocols used in these simulations by direct comparisons with experimental data. The calculated MD behavior, chain disorder, and lipid packing parameter, i.e. the ratio of the effective areas of hydrocarbon tails and head group per lipid (a_t/a_h), correctly predict the aggregation preferences of the three lipids observed experimentally at 37°C, namely: a gel bilayer for DLPE, a hexagonal tube for DOPE, and a liquid crystalline bilayer for DOPC. In addition, the model and conditions used in the MD simulations led to good agreement of the calculated properties of the bilayers with available experimental results, demonstrating the reliability of the simulations. The effects of the cis unsaturation in the hydrocarbon chains of DOPE and DOPC, compared to the fully saturated one in DLPE, as well as the effects of the different polar head groups of PC and PE with the same unsaturated chains on the lipid packing and bilayer structure have been investigated. The results of these studies indicate the ability of MD methods to provide molecular-level insights into the structure and dynamics of lipid assemblies.     

Molecular Dynamics Simulations for Metallic Nanosystems

Nanotechnology is a crucial field for future scientific development where many different disciplines meet. Computational modelization of nanometer-sized structures is a key issue in this development because (i) it allows a considerable saving of resources and costly experimental setups intended to fabricate nanometric test devices and (ii) nowadays the study of nanometric sized systems is feasible with thoroughly designed computational codes and relatively low cost computational resources. This article describes how molecular dynamics simulations, in combination with potentials obtained in the framework of the embedded atom method, are able to describe the properties of two systems of interest for the development of future nanoelectronic devices: metallic nanowires and metallic nanofilms. Our results show that nanowire stretching results in a series of well-defined geometric structures (shells) and that thin films experiment a crystallographic phase transition for a decreasing number of layers. In both cases, good agreement with experiments is found.     

Lipid Configurations from Molecular Dynamics Simulations

The extent to which current force fields faithfully reproduce conformational properties of lipids in bilayer membranes, and whether these reflect the structural principles established for phospholipids in bilayer crystals, are central to biomembrane simulations. We determine the distribution of dihedral angles in palmitoyl-oleoyl phosphatidylcholine from molecular dynamics simulations of hydrated fluid bilayer membranes. We compare results from the widely used lipid force field of Berger et al. with those from the most recent C36 release of the CHARMM force field for lipids. Only the CHARMM force field produces the chain inequivalence with sn-1 as leading chain that is characteristic of glycerolipid packing in fluid bilayers. The exposure and high partial charge of the backbone carbonyls in Berger lipids leads to artifactual binding of Na⁺ ions reported in the literature. Both force fields predict coupled, near-symmetrical distributions of headgroup dihedral angles, which is compatible with models of interconverting mirror-image conformations used originally to interpret NMR order parameters. The Berger force field produces rotamer populations that correspond to the headgroup conformation found in a phosphatidylcholine lipid bilayer crystal, whereas CHARMM36 rotamer populations are closer to the more relaxed crystal conformations of phosphatidylethanolamine and glycerophosphocholine. CHARMM36 alone predicts the correct relative signs of the time-average headgroup order parameters, and reasonably reproduces the full range of NMR data from the phosphate diester to the choline methyls. There is strong motivation to seek further experimental criteria for verifying predicted conformational distributions in the choline headgroup, including the ³¹P chemical shift anisotropy and ¹⁴N and CD₃ NMR quadrupole splittings.     

Molecular Dynamics Simulations of Lipid Membrane Electroporation

The permeability of cell membranes can be transiently increased following the application of external electric fields. Theoretical approaches such as molecular modeling provide a significant insight into the processes affecting, at the molecular level, the integrity of lipid cell membranes when these are subject to voltage gradients under similar conditions as those used in experiments. This article reports on the progress made so far using such simulations to model membrane—lipid bilayer—electroporation. We first describe the methods devised to perform in silico experiments of membranes subject to nanosecond, megavolt-per-meter pulsed electric fields and of membranes subject to charge imbalance, mimicking therefore the application of low-voltage, long-duration pulses. We show then that, at the molecular level, the two types of pulses produce similar effects: provided the TM voltage these pulses create are higher than a certain threshold, hydrophilic pores stabilized by the membrane lipid headgroups form within the nanosecond time scale across the lipid core. Similarly, when the pulses are switched off, the pores collapse (close) within similar time scales. It is shown that for similar TM voltages applied, both methods induce similar electric field distributions within the membrane core. The cascade of events following the application of the pulses, and taking place at the membrane, is a direct consequence of such an electric field distribution.     

Cryo-Cooling Effect on DHFR Crystal Studied by Replica-Exchange Molecular Dynamics Simulations

Cryo-cooling is routinely performed before x-ray diffraction image collection to reduce the damage to crystals due to ionizing radiation. It has been suggested that although backbone structures are usually very similar between room temperature and cryo-temperature, cryo-cooling may hamper biologically relevant dynamics. In this study, the crystal of Escherichia coli dihydrofolate reductase is studied with replica-exchange molecular dynamics simulation, and the results are compared with the crystal structure determined at cryo-temperature and room temperature with the time-averaged ensemble method. Although temperature dependence of unit cell compaction and root mean-square fluctuation of Cα is found in accord with experiment, it is found that the protein structure at low temperature can be more heterogeneous than the ensemble of structures reported by using the time-averaged ensemble method, encouraging further development of the time-averaged ensemble method and indicating that data should be examined carefully to avoid overinterpretation of one average structure.     

Molecular Dynamics Simulations of DNA-Polycation Complex Formation

Complexes formed from DNA and polycations are of interest because of their potential use in gene therapy; however, there remains a lack of understanding of the structure and formation of DNA-polycation complexes at atomic scale. In this work, molecular dynamics simulations of the DNA duplex d(CGCGAATTCGCG) in the presence of polycation chains are carried out to shed light on the specific atomic interaction that result in complex formation. The structures of complexes formed from DNA with polyethylenimine, which is considered one of the most promising DNA vector candidates, and a second polycation, poly-L-lysine, are compared. After an initial separation of ∼50 Å, the DNA and polycation come together and form a stable complex within 10 ns. The DNA does not undergo any major structural changes on complexation and remains in the B-form. In the formed complex, the charged amine groups of the polycation mainly interact with DNA phosphate groups, with polycation intrusion into the major and minor grooves dependent on the identity and charge state of the polycation. The ability of the polycation to effectively neutralize the charge of the DNA phosphate groups and the resulting influence on the DNA helix interaction are discussed.     

Molecular Dynamics Simulations of some Small Organic Molecules

Abstract

Free energy differences between different conformers of D-ribofuranose, L-malic acid and meso-tartaric acid in solution were calculated using Molecular Dynamics simulations. In case of ribose the α → β transition was studied. For the acids attention was focussed on the transitions between the three possible staggered conformers with respect to the central C-C bond. In all cases a thermodynamic integration method was employed to evaluate the free energy difference. The use of an alternative technique, umbrella sampling, for ribose did not give promising results.

It was shown that one needs a fairly accurate picture of the accessible conformational space in case of flexible molecules like the ones considered here before one can determine meaningful free energy differences. Large hysteresis effects between forward and reverse simulated transitions were observed, but contrary to the general belief they are no direct measure of the accuracy of the calculated ΔG values. In all cases the ΔG values resulting from the simulations and from NMR experiments agree within the, considerable, error limits and for the different forms of D-ribose, L-malic acid and L-tartaric acid the relative order of their populations is also correctly reproduced.     

Molecular Dynamics Simulations of Calcium-Free Calmodulin in Solution

Abstract

A 4-ns molecular dynamics simulation of calcium-free calmodulin in solution has been performed, using Ewald summation to treat electrostatic interactions. Our simulation results were mostly consistent with solution experimental studies, including NMR, fluorescence and x-ray scattering. The secondary structures within the N- and C-terminal domains were conserved in the simulation, with trajectory structures similar to the NMR-derived model structure 1CFD. However, the relative orientations of the domains, for which there are no NMR restraints, differed in details between the simulation and the 1CFD model. The most interesting information provided by the simulations is that the dynamics of calcium-free calmod- ulin in solution is dominated by slow rigid body reorientations of the domains. The interdomain distance fluctuated between 29 and 39 Å, and interdomain orientation angle, defined as the pseudo-dihedral formed by the four calcium binding sites, varied between ?2° and 108°. Similarly, the domain linker region also exhibited significant fluctuations, with its length varying in the 34–45 Å range and its bend angle in the 10–100° range. The simulations are in accord with fluorescence results suggesting that calcium-free calmodulin is more compact and more flexible than the calcium activated form. Surprisingly, quite similar solvent accessibilities of the hydrophobic patches were seen in the calcium-free trajectory described in this work and previously generated calcium-loaded calmodulin simulations. Thus, our simulations suggest a reexamination of the standard model of the structural change of calmodulin upon calcium binding, involving exposure of the hydrophobic patches to solvent.     

Molecular Dynamics Simulations of Molecules in Uniform Flow

Flow at the molecular level induces shear-induced unfolding of single proteins and can drive their assembly, the mechanisms of which are not completely understood. To be able to analyze the role of flow on molecules, we present uniform-flow molecular dynamics simulations at atomic level. The pull module of the GRoningen MAchine for Chemical Simulations package was extended to be able to force-group atoms within a defined layer of the simulation box. Application of this external enforcement to explicit water molecules, together with the coupling to a thermostat, led to a uniform terminal velocity of the solvent water molecules. We monitored the density of the whole system to establish the conditions under which the simulated flow is well-behaved. A maximal velocity of 1.3 m/s can be generated if a pull slice of 8 nm is used, and high velocities would require larger pull slices to still maintain a stable density. As expected, the target velocity increases linearly with the total external force applied. Finally, we suggest an appropriate setup to stretch a protein by uniform flow, in which protein extensions depend on the flow conditions. Our implementation provides an efficient computational tool to investigate the effect of the flow at the molecular level.     

An Iterative Variable-timestep Algorithm for Molecular Dynamics Simulations

A method for performing variable-timestep molecular dynamics integration is described, in which an iterative algorithm is used to select the largest timestep consistent with the desired simulation accuracy. Accuracy in this context is defined in terms of energy conservation, rather than trajectory correctness. Specifically, a timestep-independent measure of the rate of "diffusion" of the total energy is used. This variable timestep approach is compared to fixed-timestep integration for three different hydrocarbon systems (polyethylene, liquid benzene and ethylene), which are modeled with a reactive bond-order potential. These systems represent both equilibrium and highly non-equilibrium systems at temperatures ranging from 298 to 2500 K. The variable-timestep method is found to be approximately twice as computationally efficient as fixed-timestep integration for the non-equilibrium sputtering of polyethylene, and the two methods were competitive for the equilibrium systems. The algorithm requires the specification of two parameters controlling the rates of timestep growth and decay, but it is found that one set of values is appropriate for all three systems studied, and there is reason to believe that the parameters are transferable to other systems. The algorithm was developed specifically for simulations involving disparate timescales, such as are encountered with the reactive bond-order model used here, but it should also prove beneficial for a wide range of molecular dynamics applications.     

Rotamer Dynamics: Analysis of Rotamers in Molecular Dynamics Simulations of Proteins

Given by χ torsional angles, rotamers describe the side-chain conformations of amino acid residues in a protein based on the rotational isomers (hence the word rotamer). Constructed rotamer libraries, based on either protein crystal structures or dynamics studies, are the tools for classifying rotamers (torsional angles) in a way that reflect their frequency in nature. Rotamer libraries are routinely used in structure modeling and evaluation. In this perspective article, we would like to encourage researchers to apply rotamer analyses beyond their traditional use. Molecular dynamics (MD) of proteins highlight the in silico behavior of molecules in solution and thus can identify favorable side-chain conformations. In this article, we used simple computational tools to study rotamer dynamics (RD) in MD simulations. First, we isolated each frame in the MD trajectories in separate Protein Data Bank files via the cpptraj module in AMBER. Then, we extracted torsional angles via the Bio3D module in R language. The classification of torsional angles was also done in R according to the penultimate rotamer library. RD analysis is useful for various applications such as protein folding, study of rotamer-rotamer relationship in protein-protein interaction, real-time correlation between secondary structures and rotamers, study of flexibility of side chains in binding site for molecular docking preparations, use of RD as guide in functional analysis and study of structural changes caused by mutations, providing parameters for improving coarse-grained MD accuracy and speed, and many others. Major challenges facing RD to emerge as a new scientific field involve the validation of results via easy, inexpensive wet-lab methods. This realm is yet to be explored.     

The Role of Molecular Dynamics Simulations for the Study of Slow Dynamics

Abstract

A two step strategy is proposed to study dynamical properties of a physical system much slower than the time scales accessible by molecular dynamics simulations. The strategy is applied to investigate the slow dynamics of supercooled liquids.     

Modelling Diffusion in Zeolites by Molecular Dynamics Simulations

Abstract

The diffusion of molecules sorbed in zeolites is of growing interest for understanding the mechanisms of chemical processes with regard to selectivity and reactivity [1].

MD simulations give insight into physical systems on the molecular level allowing to study and visualize the motion of molecules even beyond the possibilities of experiments [2,3]. Single system parameters can easily be varied to study their influence, also those parameters that are fixed in reality (e.g., the size of particles). We present a cross section of our recent work to illustrate the capabilities of MD: The self diffusion coefficients (D) of a mixture of methane and xenon in silicalite show remarkable deviations from those of the pure species. This is shown and confirmed by PFG NMR experiments [4].

Simulating ethane in zeolite A the mechanism of diffusion has been studied. The effects of rotation on the diffusion lead to cases where D decreases with growing temperature [5].

The independence of self diffusion on lattice vibrations is proven even for zeolites with windows of guest particle size comparing simulations with rigid and vibrating zeolite lattice [6].     

Structural Dynamics of Human Telomeric G-Quadruplex Loops Studied by Molecular Dynamics Simulations

Loops which are linkers connecting G-strands and supporting the G-tetrad core in G-quadruplex are important for biological roles of G-quadruplexes. TTA loop is a common sequence which mainly resides in human telomeric DNA (hTel) G-quadruplex. A series of molecular dynamics (MD) simulations were carried out to investigate the structural dynamics of TTA loops. We found that (1) the TA base pair formed in TTA loops are very stable, the occupied of all hydrogen bonds are more than 0.95. (2) The TA base pair makes the adjacent G-quartet more stable than others. (3) For the edgewise loop and the diagonal loop, most loop bases are stacking with others, only few bases have considerable freedom. (4) The stabilities of these stacking structures are distinct. Part of the loops, especially TA base pairs, and bases stacking with the G-quartet, maintain certain stable conformations in the simulation, but other parts, like TT and TA stacking structures, are not stable enough. For the first time, spontaneous conformational switches of TTA edgewise loops were observed in our long time MD simulations. (5) For double chain reversal loop, it is really hard to maintain a stable conformation in the long time simulation under present force fields (parm99 and parmbsc0), as it has multiple conformations with similar free energies.     

Estimation of Solvent Diffusion Coefficients Using Molecular Dynamics Simulations

Abstract

Molecular dynamics simulations have been used to investigate diffusion in two commonly used industrial solvents, toluene and tetrahydrofuran. Several different models for the solvents are compared (flexible vs. rigid, all-atom vs. united atom), and it is found that united atom and all-atom models of the solvents produce very different diffusion coefficients at the experimental density. This disagreement can be explained by the pressure dependence of the diffusion coefficient, which is found to vary in accord with the Chapman-Enskog result for hard spheres. It is recommended that force fields be parametrized carefully to produce reasonable pressures at the experimental densities, or that simulations be carried out at constant pressure, if they are to be used for the purposes of calculating transport coefficients.     

Molecular Dynamics Simulations of a Calmodulin-peptide Complex in Solution

Abstract

We have performed an 4-ns MD simulation of calmodulin complexed with a target peptide in explicit water, under realistic conditions of constant temperature and pressure, in the presence of a physiological concentration of counterions and using Ewald summation to avoid truncation of long-range electrostatic forces. During the simulation the system tended to perform small fluctuations around a structure similar to, but somewhat looser than the starting crystal structure. The calmodulin-peptide complex was quite rigid and did not exhibit any large amplitude domain motions such as previously seen in apo- and calcium-bound calmodulin. We analyzed the calmodulin-peptide interactions by calculating buried surface areas, CHARMM interaction energies and continuum model interaction free energies. In the trajectory, the protein surface area buried by contact with the peptide is 1373 Ų, approximately evenly divided between the calmodulin N-terminal, C-terminal and central linker regions. A majority of this buried surface, 803 ·A², comes from nonpolar residues, in contrast to the protein as a whole, for which the surface is made up of mostly polar and charged groups. Our continuum calculations indicate that the largest favorable contribution to pep- tide binding comes from burial of molecular surface upon complex formation. Electrostatic contributions are favorable but smaller in the trajectory structures, and actually unfavorable for binding in the crystal structure. Since nonpolar groups make up most of buried surface of the protein, our calculations suggest that the hydrophobic effect is the main driving force for binding the helical peptide to calmodulin, consistent with thermodynamic analysis of experimental data. Besides the burial of nonpolar surface area, secondary contributions to peptide binding come from burial of polar surface and electrostatic interactions. In the nonpolar interactions a crucial role is played by the nine methionines of calmodulin. In the electrostatic interactions the negatively charged protein residues and positively charged peptide residues play a dominant role.