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.     

A simplified model of local structure in aqueous proline amino acid revealed by first-principles molecular dynamics simulations

Aqueous proline solutions are deceptively simple as they can take on complex roles such as protein chaperones, cryoprotectants, and hydrotropic agents in biological processes. Here, a molecular level picture of proline/water mixtures is developed. Car-Parrinello ab initio molecular dynamics (CPAIMD) simulations of aqueous proline amino acid at the B-LYP level of theory, performed using IBM's Blue Gene/L supercomputer and massively parallel software, reveal hydrogen-bonding propensities that are at odds with the predictions of the CHARMM22 empirical force field but are in better agreement with results of recent neutron diffraction experiments. In general, the CPAIMD (B-LYP) simulations predict a simplified structural model of proline/water mixtures consisting of fewer distinct local motifs. Comparisons of simulation results to experiment are made by direct evaluation of the neutron static structure factor S(Q) from CPAIMD (B-LYP) trajectories as well as to the results of the empirical potential structure refinement reverse Monte Carlo procedure applied to the neutron data.     

Polarizability rescaling and atom-based Thole scaling in the CHARMM Drude polarizable force field for ethers

Within the CHARMM polarizable force field based on the classical Drude oscillator, atomic polarizabilities are derived via fitting to ab initio calculated data on isolated gas phase molecules, with an empirical scaling factor applied to account for differences between the gas and condensed phases. In the development of polarizable models for the ethers, a polarizability scaling factor of 0.7 was previously applied [Vorobyov et al. J Comput Chem 3:1120–1133, 2007]. While the resulting force field models gave good agreement with a variety of experimental data, they systematically underestimated the liquid phase dielectric constants. Here, a new CHARMM polarizable model is developed for the ethers, employing a polarizability scaling factor of 0.85 and including atom-based Thole scale factors recently introduced into the CHARMM Drude polarizable force field [Harder et al. J Phys Chem B 112:3509-3521, 2008]. The new model offers a significant improvement in the reproduction of liquid phase dielectric constants, while maintaining the good agreement of the previous model with all other experimental and quantum mechanical data, highlighting the sensitivity of liquid phase properties to the choice of atomic polarizability parameters.     

Reparameterization of all-atom dipalmitoylphosphatidylcholine lipid parameters enables simulation of fluid bilayers at zero tension

Molecular dynamics simulations of dipalmitoylphosphatidylcholine (DPPC) lipid bilayers using the CHARMM27 force field in the tensionless isothermal-isobaric (NPT) ensemble give highly ordered, gel-like bilayers with an area per lipid of approximately 48 A(2). To obtain fluid (L(alpha)) phase properties of DPPC bilayers represented by the CHARMM energy function in this ensemble, we reparameterized the atomic partial charges in the lipid headgroup and upper parts of the acyl chains. The new charges were determined from the electron structure using both the Mulliken method and the restricted electrostatic potential fitting method. We tested the derived charges in molecular dynamics simulations of a fully hydrated DPPC bilayer. Only the simulation with the new restricted electrostatic potential charges shows significant improvements compared with simulations using the original CHARMM27 force field resulting in an area per lipid of 60.4 +/- 0.1 A(2). Compared to the 48 A(2), the new value of 60.4 A(2) is in fair agreement with the experimental value of 64 A(2). In addition, the simulated order parameter profile and electron density profile are in satisfactory agreement with experimental data. Thus, the biologically more interesting fluid phase of DPPC bilayers can now be simulated in all-atom simulations in the NPT ensemble by employing our modified CHARMM27 force field.     

Structural equilibrium of DNA represented with different force fields.

We have recently indicated preliminary evidence of different equilibrium average structures with the CHARMM and AMBER force fields in explicit solvent molecular dynamics simulations on the DNA duplex d(C5T5) . d(A5G5) (Feig, M. and B.M. Pettitt, 1997, Experiment vs. Force Fields: DNA conformation from molecular dynamics simulations. J. Phys. Chem. B. (101:7361-7363). This paper presents a detailed comparison of DNA structure and dynamics for both force fields from extended simulation times of 10 ns each. Average structures display an A-DNA base geometry with the CHARMM force field and a base geometry that is intermediate between A- and B-DNA with the AMBER force field. The backbone assumes B form on both strands with the AMBER force field, while the CHARMM force field produces heterogeneous structures with the purine strand in A form and the pyrimidine strand in dynamical equilibrium between A and B conformations. The results compare well with experimental data for the cytosine/guanine part but fail to fully reproduce an overall B conformation in the thymine/adenine tract expected from crystallographic data, particularly with the CHARMM force field. Fluctuations between A and B conformations are observed on the nanosecond time scale in both simulations, particularly with the AMBER force field. Different dynamical behavior during the first 4 ns indicates that convergence times of several nanoseconds are necessary to fully establish a dynamical equilibrium in all structural quantities on the time scale of the simulations presented here.     

DNA polymorphism: a comparison of force fields for nucleic acids

The improvements of the force fields and the more accurate treatment of long-range interactions are providing more reliable molecular dynamics simulations of nucleic acids. The abilities of certain nucleic acid force fields to represent the structural and conformational properties of nucleic acids in solution are compared. The force fields are AMBER 4.1, BMS, CHARMM22, and CHARMM27; the comparison of the latter two is the primary focus of this paper. The performance of each force field is evaluated first on its ability to reproduce the B-DNA decamer d(CGATTAATCG)(2) in solution with simulations in which the long-range electrostatics were treated by the particle mesh Ewald method; the crystal structure determined by Quintana et al. (1992) is used as the starting point for all simulations. A detailed analysis of the structural and solvation properties shows how well the different force fields can reproduce sequence-specific features. The results are compared with data from experimental and previous theoretical studies.     

Effects of force fields on the conformational and dynamic properties of amyloid β(1‐40) dimer explored by replica exchange molecular dynamics simulations

The conformational space and structural ensembles of amyloid beta (Aβ) peptides and their oligomers in solution are inherently disordered and proven to be challenging to study. Optimum force field selection for molecular dynamics (MD) simulations and the biophysical relevance of results are still unknown. We compared the conformational space of the Aβ(1‐40) dimers by 300 ns replica exchange MD simulations at physiological temperature (310 K) using: the AMBER‐ff99sb‐ILDN, AMBER‐ff99sb*‐ILDN, AMBER‐ff99sb‐NMR, and CHARMM22* force fields. Statistical comparisons of simulation results to experimental data and previously published simulations utilizing the CHARMM22* and CHARMM36 force fields were performed. All force fields yield sampled ensembles of conformations with collision cross sectional areas for the dimer that are statistically significantly larger than experimental results. All force fields, with the exception of AMBER‐ff99sb‐ILDN (8.8 ± 6.4%) and CHARMM36 (2.7 ± 4.2%), tend to overestimate the α‐helical content compared to experimental CD (5.3 ± 5.2%). Using the AMBER‐ff99sb‐NMR force field resulted in the greatest degree of variance (41.3 ± 12.9%). Except for the AMBER‐ff99sb‐NMR force field, the others tended to under estimate the expected amount of β‐sheet and over estimate the amount of turn/bend/random coil conformations. All force fields, with the exception AMBER‐ff99sb‐NMR, reproduce a theoretically expected β‐sheet‐turn‐β‐sheet conformational motif, however, only the CHARMM22* and CHARMM36 force fields yield results compatible with collapse of the central and C‐terminal hydrophobic cores from residues 17‐21 and 30‐36. Although analyses of essential subspace sampling showed only minor variations between force fields, secondary structures of lowest energy conformers are different.     

Refinement of homology-based protein structures by molecular dynamics simulation techniques

The use of classical molecular dynamics simulations, performed in explicit water, for the refinement of structural models of proteins generated ab initio or based on homology has been investigated. The study involved a test set of 15 proteins that were previously used by Baker and coworkers to assess the efficiency of the ROSETTA method for ab initio protein structure prediction. For each protein, four models generated using the ROSETTA procedure were simulated for periods of between 5 and 400 nsec in explicit solvent, under identical conditions. In addition, the experimentally determined structure and the experimentally derived structure in which the side chains of all residues had been deleted and then regenerated using the WHATIF program were simulated and used as controls. A significant improvement in the deviation of the model structures from the experimentally determined structures was observed in several cases. In addition, it was found that in certain cases in which the experimental structure deviated rapidly from the initial structure in the simulations, indicating internal strain, the structures were more stable after regenerating the side-chain positions. Overall, the results indicate that molecular dynamics simulations on a tens to hundreds of nanoseconds time scale are useful for the refinement of homology or ab initio models of small to medium-size proteins.     

Study of the betulin molecule in a water environment; ab initio and molecular simulation calculations

Ab initio and molecular simulation methods were used in calculations of the neutral individual betulin molecule, and molecular simulations were used to optimize the betulin molecule immersed in various amounts of water. Individual betulin was optimized in different force fields to find the one exhibiting best agreement with ab initio calculations obtained in the Gaussian03 program. Dihedral torsions of active groups of betulin were determined for both procedures, and related calculated structures were compared successfully. The selected force field was used for subsequent optimization of betulin in a water environment, and a conformational search was performed using quench molecular dynamics. The total energies of betulin and its interactions in water bulk were calculated, and the influence of water on betulin structure was investigated.     

Comparison of force fields for Alzheimer's A : A case study for intrinsically disordered proteins

Intrinsically disordered proteins are essential for biological processes such as cell signalling, but are also associated to devastating diseases including Alzheimer's disease, Parkinson's disease or type II diabetes. Because of their lack of a stable three‐dimensional structure, molecular dynamics simulations are often used to obtain atomistic details that cannot be observed experimentally. The applicability of molecular dynamics simulations depends on the accuracy of the force field chosen to represent the underlying free energy surface of the system. Here, we use replica exchange molecular dynamics simulations to test five modern force fields, OPLS, AMBER99SB, AMBER99SB*ILDN, AMBER99SBILDN‐NMR and CHARMM22*, in their ability to model Aβ₄₂, an intrinsically disordered peptide associated with Alzheimer's disease, and compare our results to nuclear magnetic resonance (NMR) experimental data. We observe that all force fields except AMBER99SBILDN‐NMR successfully reproduce local NMR observables, with CHARMM22* being slightly better than the other force fields.     

A modified version of the Cornell et al. force field with improved sugar pucker phases and helical repeat

We have examined some subtle parameter modifications to the Cornell et al. force field, which has proven quite successful in reproducing nucleic acid properties, but whose C2'-endo sugar pucker phase and helical repeat for B DNA appear to be somewhat underestimated. Encouragingly, the addition of a single V2 term involving the atoms C(sp3)-O-(sp3)-C(sp3)-N(sp2), which can be nicely rationalized because of the anomeric effect (lone pairs on oxygen are preferentially oriented relative to the electron withdrawing N), brings the sugar pucker phase of C2'-endo sugars to near perfect agreement with ab initio calculations (W near 162 degrees). Secondly, the use of high level ab initio calculations on entire nucleosides (in contrast to smaller model systems necessitated in 1994-95 by computer limitations) lets one improve the chi torsional potential for nucleic acids. Finally, the O(sp3)-C(sp3)- C(sp3)-O(sp3) V2 torsional potential has been empirically adjusted to reproduce the ab initio calculated relative energy of C2'-endo and C3'-endo nucleosides. These modifications are tested in molecular dynamics simulations of mononucleosides (to assess sugar pucker percentages) and double helices of DNA and RNA (to assess helical and sequence specific structural properties). In both areas, the modified force field leads to improved agreement with experimental data.     

Influence of non-bonded parameters on the quality of NMR structures: A new force field for NMR structure calculation

The effects of different non-bonded parameters of force fields for NMR structure calculation on the quality of the resulting NMR solution structures were investigated using Interleukin 4 as a model system. NMR structure ensembles were calculated with an ab initio protocol using torsion angle dynamics. The calculations were repeated with five different non-bonded energy functions and parameters. The resulting ensembles were compared with the available X-ray structures, and their quality was assessed with common structure validation programs. In addition, the impact of torsion angle restraints and dihedral energy terms for the sidechains and the backbone was studied. The further improvement of the quality by refinement in explicit solvent was demonstrated. The optimal parameters, including those necessary for water refinement, are available in the new version of the PARALLHDG force field.     

Molecular dynamics simulations of the unfolding of beta(2)-microglobulin and its variants

In this study, we examined the unfolding processes of native beta(2)-microglobulin and two related variants, one with an N-terminal hexapeptide deletion DeltaN6 and another with Lys57-Asp58 cleavage, by high-temperature molecular dynamics simulations. Three simulation models were used, molecular dynamics (MD) simulations with explicit water solvation, MD simulations with the CHARMM EEF1 force field and Langevin dynamics with the CHARMM EEF1 force field. Our simulations reproduce many of the experimentally observed structural changes. The most striking agreement is in the beta-strands to alpha-helix transition. In our simulations, strands beta(3), beta(4) and beta(5) consistently change to alpha-helix, whereas beta(8) changes to an alpha-helix only briefly. Through comparisons of the conformational behavior of the native, the DeltaN6 and the Lys-cut beta(2)-m, using the three simulation methods, we identified the consensus conformational changes that differentiate between the native beta(2)-m and its two variants. We found that the main effect of the removal of the N-terminal hexapeptide is to increase the separation between strands beta(2) and beta(6) and to facilitate the beta to alpha transition. On the other hand, the lysine cleavage only increases the flexibility of strand beta(5) and does not affect the interactions between strands beta(2) and beta(6). These conformational changes may relate to polymerization tendencies of these variants.     

Experimental validation of molecular dynamics simulations of lipid bilayers: a new approach

A novel protocol has been developed for comparing the structural properties of lipid bilayers determined by simulation with those determined by diffraction experiments, which makes it possible to test critically the ability of molecular dynamics simulations to reproduce experimental data. This model-independent method consists of analyzing data from molecular dynamics bilayer simulations in the same way as experimental data by determining the structure factors of the system and, via Fourier reconstruction, the overall transbilayer scattering-density profiles. Multi-nanosecond molecular dynamics simulations of a dioleoylphosphatidylcholine bilayer at 66% RH (5.4 waters/lipid) were performed in the constant pressure and temperature ensemble using the united-atom GROMACS and the all-atom CHARMM22/27 force fields with the GROMACS and NAMD software packages, respectively. The quality of the simulated bilayer structures was evaluated by comparing simulation with experimental results for bilayer thickness, area/lipid, individual molecular-component distributions, continuous and discrete structure factors, and overall scattering-density profiles. Neither the GROMACS nor the CHARMM22/27 simulations reproduced experimental data within experimental error. The widths of the simulated terminal methyl distributions showed a particularly strong disagreement with the experimentally observed distributions. A comparison of the older CHARMM22 with the newer CHARMM27 force fields shows that significant progress is being made in the development of atomic force fields for describing lipid bilayer systems empirically.     

Harmonic force field for nitro compounds

Molecular simulations leading to sensors for the detection of explosive compounds require force field parameters that can reproduce the mechanical and vibrational properties of energetic materials. We developed precise harmonic force fields for alanine polypeptides and glycine oligopeptides using the FUERZA procedure that uses the Hessian tensor (obtained from ab initio calculations) to calculate precise parameters. In this work, we used the same procedure to calculate generalized force field parameters of several nitro compounds. We found a linear relationship between force constant and bond distance. The average angle in the nitro compounds was 116°, excluding the 90° angle of the carbon atoms in the octanitrocubane. The calculated parameters permitted the accurate molecular modeling of nitro compounds containing many functional groups. Results were acceptable when compared with others obtained using methods that are specific for one type of molecule, and much better than others obtained using methods that are too general (these ignore the chemical effects of surrounding atoms on the bonding and therefore the bond strength, which affects the mechanical and vibrational properties of the whole molecule).     

Interchain hydrophobic clustering promotes rigidity in HIV-1 protease flap dynamics: new insights from Molecular Dynamics

The dynamics of HIV-1 protease (HIV-pr), a drug target for HIV infection, has been studied extensively by both computational and experimental methods. The flap dynamics of HIV-pr is considered to be more important for better ligand binding and enzymatic actions. Moreover, it has been demonstrated that the drug-induced mutations can change the flap dynamics of HIV-pr affecting the binding affinity of the ligands. Therefore, detailed understanding of flap dynamics is essential for designing better inhibitors. Previous computational investigations observed significant variation in the flap opening in nanosecond time scale indicating that the dynamics is highly sensitive to the simulation protocols. To understand the sensitivity of the flap dynamics on the force field and simulation protocol, molecular dynamics simulations of HIV-pr have been performed with two different AMBER force fields, ff99 and ff02. Two different trajectories (20?ns each) were obtained using the ff99 and ff02 force field. The results showed polarizable force field (ff02) make the flap tighter than the nonpolarizable force field (ff99). Some polar interactions and hydrogen bonds involving flap residues were found to be stronger with ff02 force field. The formation of interchain hydrophobic cluster (between flap tip of one chain and active site wall of another chain) was found to be dominant in the semi-open structures obtained from the simulations irrespective of the force field. It is proposed that an inhibitor, which will promote this interchain hydrophobic clustering, may make the flaps more rigid, and presumably the effect of mutation would be small on ligand binding.     

Inclusion of Many-Body Effects in the Additive CHARMM Protein CMAP Potential Results in Enhanced Cooperativity of α-Helix and β-Hairpin Formation

Folding simulations on peptides and proteins using empirical force fields have demonstrated the sensitivity of the results to details of the backbone potential. A recently revised version of the additive CHARMM protein force field, which includes optimization of the backbone CMAP potential to achieve good balance between different types of secondary structure, correcting the α-helical bias present in the former CHARMM22/CMAP energy function, is shown to result in improved cooperativity for the helix-coil transition. This is due to retention of the empirical corrections introduced in the original CMAP to reproduce folded protein structures—corrections that capture many-body effects missing from an energy surface fitted to gas phase calculations on dipeptides. The experimental temperature dependence of helix formation in (AAQAA)₃ and parameters for helix nucleation and elongation are in much better agreement with experiment than those obtained with other recent force fields. In contrast, CMAP parameters derived by fitting to a vacuum quantum mechanical surface for the alanine dipeptide do not reproduce the enhanced cooperativity, showing that the empirical backbone corrections, and not some other feature of the force field, are responsible. We also find that the cooperativity of β-hairpin formation is much improved relative to other force fields we have studied. Comparison with (ϕ,ψ) distributions from the Protein Data Bank further justifies the inclusion of many-body effects in the CMAP. These results suggest that the revised energy function will be suitable for both simulations of unfolded or intrinsically disordered proteins and for investigating protein-folding mechanisms.     

The range and shielding of dipole-dipole interactions in phospholipid bilayers

In molecular dynamics simulations of lipid bilayers, the structure is sensitive to the precise treatment of electrostatics. The dipole-dipole interactions between headgroup dipoles are not long-ranged, but the area per lipid and, through it, other properties of the bilayer are very sensitive to the detailed balance between the perpendicular and in-plane components of the headgroup dipoles. This is affected by the detailed properties of the cutoff scheme or if long-range interactions are included by Ewald or particle-mesh Ewald techniques. Interaction between the in-plane components of the headgroup dipoles is attractive and decays as the inverse sixth power of distance. The interaction is screened by the square of a dielectric permittivity close to the value for water. Interaction between the components perpendicular to the membrane plane is repulsive and decays as the inverse third power of distance. These interactions are screened by a dielectric permittivity of the order 10. Thus, despite the perpendicular components being much smaller in magnitude than the in-plane components, they will dominate the interaction energies at large distances.