CHARMM-GUI Nanodisc Builder for modeling and simulation of various nanodisc systems

Nanodiscs are discoidal protein–lipid complexes that have wide applications in membrane protein studies. Modeling and simulation of nanodiscs are challenging due to the absence of structures of many membrane scaffold proteins (MSPs) that wrap around the membrane bilayer. We have developed CHARMM-GUI Nanodisc Builder ( http://www.charmm-gui.org/input/nanodisc ) to facilitate the setup of nanodisc simulation systems by modeling the MSPs with defined size and known structural features. A total of 11 different nanodiscs with a diameter from 80 to 180 Å are made available in both the all-atom CHARMM and two coarse-grained (PACE and Martini) force fields. The usage of the Nanodisc Builder is demonstrated with various simulation systems. The structures and dynamics of proteins and lipids in these systems were analyzed, showing similar behaviors to those from previous all-atom and coarse-grained nanodisc simulations. We expect the Nanodisc Builder to be a convenient and reliable tool for modeling and simulation of nanodisc systems. © 2019 Wiley Periodicals, Inc.     

Combination of the CHARMM27 force field with united-atom lipid force fields

Computer simulations offer a valuable way to study membrane systems, from simple lipid bilayers to large transmembrane protein complexes and lipid-nucleic acid complexes for drug delivery. Their accuracy depends on the quality of the force field parameters used to describe the components of a particular system. We have implemented the widely used CHARMM22 and CHARMM27 force fields in the GROMACS simulation package to (i) combine the CHARMM22 protein force field with two sets of united-atom lipids parameters; (ii) allow comparisons of the lipid CHARMM27 force field with other lipid force fields or lipid-protein force field combinations. Our tests do not show any particular issue with the combination of the all-atom CHARMM22 force field with united-atoms lipid parameters, although pertinent experimental data are lacking to assess the quality of the lipid-protein interactions. The conversion utilities allow automatic generation of GROMACS simulation files with CHARMM nucleic acids and protein parameters and topologies, starting from pdb files using the standard GROMACS pdb2gmx method. CMAP is currently not implemented.     

Force field dependence of phospholipid headgroup and acyl chain properties: Comparative molecular dynamics simulations of DMPC bilayers

The reliability of molecular simulations largely depends on the quality of the empirical force field parameters. Force fields used in lipid simulations continue to be improved to enhance the agreement with experiments for a number of different properties. In this work, we have carried out molecular dynamics simulations of neat DMPC bilayers using united‐atom Berger force field and three versions of all‐atom CHARMM force fields. Three different systems consisting of 48, 72, and 96 lipids were studied. Both particle mesh Ewald (PME) and spherical cut‐off schemes were used to evaluate the long‐range electrostatic interactions. In total, 21 simulations were carried out and analyzed to find out the dependence of lipid properties on force fields, system size, and schemes to calculate long‐range interactions. The acyl chain order parameters calculated from Berger and the recent versions of CHARMM simulations have shown generally good agreement with the experimental results. However, both sets of force fields deviate significantly from the experimentally observed P‐C dipolar coupling values for the carbon atoms that link the choline and glycerol groups with the phosphate groups. Significant differences are also observed in several headgroup parameters between CHARMM and Berger simulations. Our results demonstrate that when changes were introduced to improve CHARMM force field using PME scheme, all the headgroup parameters have not been reoptimized. The headgroup properties are likely to play a significant role in lipid–lipid, protein–lipid, and ligand–lipid interactions and hence headgroup parameters in phospholipids require refinement for both Berger and CHARMM force fields. © 2009 Wiley Periodicals, Inc.J Comput Chem, 2010     

Exploring the inter-molecular interactions in amyloid-β protofibril with molecular dynamics simulations and molecular mechanics Poisson-Boltzmann surface area free energy calculations

Aggregation of amyloid-β (Aβ) peptides correlates with the pathology of Alzheimer's disease. However, the inter-molecular interactions between Aβ protofibril remain elusive. Herein, molecular mechanics Poisson-Boltzmann surface area analysis based on all-atom molecular dynamics simulations was performed to study the inter-molecular interactions in Aβ(17-42) protofibril. It is found that the nonpolar interactions are the important forces to stabilize the Aβ(17-42) protofibril, while electrostatic interactions play a minor role. Through free energy decomposition, 18 residues of the Aβ(17-42) are identified to provide interaction energy lower than -2.5 kcal/mol. The nonpolar interactions are mainly provided by the main chain of the peptide and the side chains of nine hydrophobic residues (Leu17, Phe19, Phe20, Leu32, Leu34, Met35, Val36, Val40, and Ile41). However, the electrostatic interactions are mainly supplied by the main chains of six hydrophobic residues (Phe19, Phe20, Val24, Met35, Val36, and Val40) and the side chains of the charged residues (Glu22, Asp23, and Lys28). In the electrostatic interactions, the overwhelming majority of hydrogen bonds involve the main chains of Aβ as well as the guanidinium group of the charged side chain of Lys28. The work has thus elucidated the molecular mechanism of the inter-molecular interactions between Aβ monomers in Aβ(17-42) protofibril, and the findings are considered critical for exploring effective agents for the inhibition of Aβ aggregation.     

Biomolecular simulations of membranes: physical properties from different force fields

Phospholipid force fields are of ample importance for the simulation of artificial bilayers, membranes, and also for the simulation of integral membrane proteins. Here, we compare the two most applied atomic force fields for phospholipids, the all-atom CHARMM27 and the united atom Berger force field, with a newly developed all-atom generalized AMBER force field (GAFF) for dioleoylphosphatidylcholine molecules. Only the latter displays the experimentally observed difference in the order of the C2 atom between the two acyl chains. The interfacial water dynamics is smoothly increased between the lipid carbonyl region and the bulk water phase for all force fields; however, the water order and with it the electrostatic potential across the bilayer showed distinct differences between the force fields. Both Berger and GAFF underestimate the lipid self-diffusion. GAFF offers a consistent force field for the atomic scale simulation of biomembranes.     

Classical and quantum mechanical/molecular mechanical molecular dynamics simulations of alanine dipeptide in water: comparisons with IR and vibrational circular dichroism spectra

We have implemented the combined quantum mechanical (QM)/molecular mechanical (MM) molecular dynamics (MD) simulations of alanine dipeptide in water along with the polarizable and nonpolarizable classical MD simulations with different models of water. For the QM/MM MD simulation, the alanine dipeptide is treated with the AM1 or PM3 approximations and the fluctuating solute dipole moment is calculated by the Mulliken population analysis. For the classical MD simulations, the solute is treated with the polarizable or nonpolarizable AMBER and polarizable CHARMM force fields and water is treated with the TIP3P, TIP4P, or TIP5P model. It is found that the relative populations of right-handed alpha-helix and extended beta and P(II) conformations in the simulation trajectory strongly depend on the simulation method. For the QM/MM MD simulations, the PM3/MM shows that the P(II) conformation is dominant, whereas the AM1/MM predicts that the dominant conformation is alpha(R). Polarizable CHARMM force field gives almost exclusively P(II) conformation and other force fields predict that both alpha-helical and extended (beta and P(II)) conformations are populated with varying extents. Solvation environment around the dipeptide is investigated by examining the radial distribution functions and numbers and lifetimes of hydrogen bonds. Comparing the simulated IR and vibrational circular dichroism spectra with experimental results, we concluded that the dipeptide adopts the P(II) conformation and PM3/MM, AMBER03 with TIP4P water, and AMBER polarizable force fields are acceptable for structure determination of the dipeptide considered in this paper.     

A nonadditive methanol force field: bulk liquid and liquid-vapor interfacial properties via molecular dynamics simulations using a fluctuating charge model

We study the bulk and interfacial properties of methanol via molecular dynamics simulations using a CHARMM (Chemistry at HARvard Molecular Mechanics) fluctuating charge force field. We discuss the parametrization of the electrostatic model as part of the ongoing CHARMM development for polarizable protein force fields. The bulk liquid properties are in agreement with available experimental data and competitive with existing fixed-charge and polarizable force fields. The liquid density and vaporization enthalpy are determined to be 0.809 g/cm3 and 8.9 kcal/mol compared to the experimental values of 0.787 g/cm3 and 8.94 kcal/mol, respectively. The liquid structure as indicated by radial distribution functions is in keeping with the most recent neutron diffraction results; the force field shows a slightly more ordered liquid, necessarily arising from the enhanced condensed phase electrostatics (as evidenced by an induced liquid phase dipole moment of 0.7 D), although the average coordination with two neighboring molecules is consistent with the experimental diffraction study as well as with recent density functional molecular dynamics calculations. The predicted surface tension of 19.66+/-1.03 dyn/cm is slightly lower than the experimental value of 22.6 dyn/cm, but still competitive with classical force fields. The interface demonstrates the preferential molecular orientation of molecules as observed via nonlinear optical spectroscopic methods. Finally, via canonical molecular dynamics simulations, we assess the model's ability to reproduce the vapor-liquid equilibrium from 298 to 423 K, the simulation data then used to obtain estimates of the model's critical temperature and density. The model predicts a critical temperature of 470.1 K and critical density of 0.312 g/cm3 compared to the experimental values of 512.65 K and 0.279 g/cm3, respectively. The model underestimates the critical temperature by 8% and overestimates the critical density by 10%, and in this sense is roughly equivalent to the underlying fixed-charge CHARMM22 force field.     

A semi-implicit solvent model for the simulation of peptides and proteins

We present a new model of biomolecules hydration based on macroscopic electrostatic theory, that can both describe the microscopic details of solvent-solute interactions and allow for an efficient evaluation of the electrostatic hydration free energy. This semi-implicit model considers the solvent as an ensemble of polarizable pseudoparticles whose induced dipole describe both the electronic and orientational solvent polarization. In the presented version of the model, there is no mutual dipolar interaction between the particles, and they only interact through short-ranged Lennard-Jones interactions. The model has been integrated into a molecular dynamics code, and offers the possibility to simulate efficiently the conformational evolution of biomolecules. It is able to provide estimations of the electrostatic solvation free energy within short time windows during the simulation. It has been applied to the study of two small peptides, the octaalanine and the N-terminal helix of ribonuclease A, and two proteins, the bovine pancreatic trypsin inhibitor and the B1 immunoglobin-binding domain of streptococcal protein G. Molecular dynamics simulations of these biomolecules, using a slightly modified Amber force field, provide stable and meaningful trajectories in overall agreement with experiments and all-atom simulations. Correlations with respect to Poisson-Boltzmann electrostatic solvation free energies are also presented to discuss the parameterization of the model and its consequences.     

CHARMM fluctuating charge force field for proteins: II protein/solvent properties from molecular dynamics simulations using a nonadditive electrostatic model

A fluctuating charge (FQ) force field is applied to molecular dynamics simulations for six small proteins in explicit polarizable solvent represented by the TIP4P-FQ potential. The proteins include 1FSV, 1ENH, 1PGB, 1VII, 1H8K, and 1CRN, representing both helical and beta-sheet secondary structural elements. Constant pressure and temperature (NPT) molecular dynamics simulations are performed on time scales of several nanoseconds, the longest simulations yet reported using explicitly polarizable all-atom empirical potentials (for both solvent and protein) in the condensed phase. In terms of structure, the FQ force field allows deviations from native structure up to 2.5 A (with a range of 1.0 to 2.5 A). This is commensurate to the performance of the CHARMM22 nonpolarizable model and other currently existing polarizable models. Importantly, secondary structural elements maintain native structure in general to within 1 A (both helix and beta-strands), again in good agreement with the nonpolarizable case. In qualitative agreement with QM/MM ab initio dynamics on crambin (Liu et al. Proteins 2001, 44, 484), there is a sequence dependence of average condensed phase atomic charge for all proteins, a dependence one would anticipate considering the differing chemical environments around individual atoms; this is a subtle quantum mechanical feature captured in the FQ model but absent in current state-of-the-art nonpolarizable models. Furthermore, there is a mutual polarization of solvent and protein in the condensed phase. Solvent dipole moment distributions within the first and second solvation shells around the protein display a shift towards higher dipole moments (increases on the order of 0.2-0.3 Debye) relative to the bulk; protein polarization is manifested via the enhanced condensed phase charges of typical polar atoms such as backbone carbonyl oxygens, amide nitrogens, and amide hydrogens. Finally, to enlarge the sample set of proteins, gas-phase minimizations and 1 ps constant temperature simulations are performed on various-sized proteins to compare to earlier work by Kaminsky et al. (J Comp Chem 2002, 23, 1515). The present work establishes the feasibility of applying a fully polarizable force field for protein simulations and demonstrates the approach employed in extending the CHARMM force field to include these effects.     

Polarizable empirical force field for sulfur‐containing compounds based on the classical Drude oscillator model

Condensed‐phase computational studies of molecules using molecular mechanics approaches require the use of force fields to describe the energetics of the systems as a function of structure. The advantage of polarizable force fields over nonpolarizable (or additive) models lies in their ability to vary their electronic distribution as a function of the environment. Toward development of a polarizable force field for biological molecules, parameters for a series of sulfur‐containing molecules are presented. Parameter optimization was performed to reproduce quantum mechanical and experimental data for gas phase properties including geometries, conformational energies, vibrational spectra, and dipole moments as well as for condensed phase properties such as heats of vaporization, molecular volumes, and free energies of hydration. Compounds in the training set include methanethiol, ethanethiol, propanethiol, ethyl methyl sulfide, and dimethyl disulfide. The molecular volumes and heats of vaporization are in good accordance with experimental values, with the polarizable model performing better than the CHARMM22 nonpolarizable force field. Improvements with the polarizable model were also obtained for molecular dipole moments and in the treatment of intermolecular interactions as a function of orientation, in part due to the presence of lone pairs and anisotropic atomic polarizability on the sulfur atoms. Significant advantage of the polarizable model was reflected in calculation of the dielectric constants, a property that CHARMM22 systematically underestimates. The ability of this polarizable model to accurately describe a range of gas and condensed phase properties paves the way for more accurate simulation studies of sulfur‐containing molecules including cysteine and methionine residues in proteins. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

Simulation of multiphase systems utilizing independent force fields to control intraphase and interphase behavior

Fixed‐charge empirical force fields have been developed and widely used over the past three decades for all‐atom molecular simulations. Most simulation programs providing these methods enable only one set of force field parameters to be used for the entire system. Whereas this is generally suitable for single‐phase systems, the molecular environment at the interface between two phases may be sufficiently different from the individual phases to require a different set of parameters to be used to accurately represent the system. Recently published simulations of peptide adsorption to material surfaces using the CHARMM force field have clearly demonstrated this issue by revealing that calculated values of adsorption free energy substantially differ from experimental results. Whereas nonbonded parameters could be adjusted to correct this problem, this cannot be done without also altering the conformational behavior of the peptide in solution, for which CHARMM has been carefully tuned. We have developed a dual‐force‐field approach (Dual‐FF) to address this problem and implemented it in the CHARMM simulation package. This Dual‐FF method provides the capability to use two separate sets of nonbonded force field parameters within the same simulation: one set to represent intraphase interactions and a separate set to represent interphase interactions. Using this approach, we show that interfacial parameters can be adjusted to correct errors in peptide adsorption free energy without altering peptide conformational behavior in solution. This program thus provides the capability to enable both intraphase and interphase molecular behavior to be accurately and efficiently modeled in the same simulation. © 2012 Wiley Periodicals, Inc.     

Peptide loop-closure kinetics from microsecond molecular dynamics simulations in explicit solvent

End-to-end contact formation rates of several peptides were recently measured by tryptophan triplet quenching (Lapidus et al. Proc. Natl. Acad. Sci. U.S.A. 2000, 97, 7220). Motivated by these experiments, we study loop-closure kinetics for two peptides of different lengths, Cys-(Ala-Gly-Gln)n-Trp (n = 1, 2), in multiple all-atom explicit-solvent molecular dynamics simulations with different initial conditions and force fields. In 150 simulations of approximately 20 ns each, we collect data covering 1.0 and 0.8 micros for the penta-peptide simulated with the AMBER and CHARMM force fields, respectively, and about 0.5 micros each with the two force fields for the octa-peptide. These extensive simulations allow us to analyze the dynamics of peptides in the unfolded state with atomic resolution, thus probing early events in protein folding, and to compare molecular dynamics simulations directly with experiment. The calculated lifetimes of the tryptophan triplet state are in the range of 50-100 ns, in agreement with experimental measurements. However, end-to-end contacts form more rapidly, with characteristic times less than 10 ns. The contact formation rates for the two force fields are similar despite differences in the respective ensembles of peptide conformations.     

The CHARMM–TURBOMOLE interface for efficient and accurate QM/MM molecular dynamics,free energies,and excited state properties

The quantum mechanical (QM)/molecular mechanical (MM) interface between Chemistry at HARvard Molecular Mechanics (CHARMM) and TURBOMOLE is described. CHARMM provides an extensive set of simulation algorithms, like molecular dynamics (MD) and free energy perturbation, and support for mature nonpolarizable and Drude polarizable force fields. TURBOMOLE provides fast QM calculations using density functional theory or wave function methods and excited state properties. CHARMM–TURBOMOLE is well‐suited for extended QM/MM MD simulations using first principles methods with large (triple‐ζ) basis sets. We demonstrate these capabilities with a QM/MM simulation of Mg²⁺(aq), where the MM outer sphere water molecules are represented using the SWM4‐NDP Drude polarizable force field and the ion and inner coordination sphere are represented using QM PBE, PBE0, and MP2 methods. The relative solvation free energies of Mg²⁺ and Zn²⁺ were calculated using thermodynamic integration. We also demonstrate the features for excited state properties. We calculate the time‐averaged solution absorption spectrum of indole, the emission spectrum of the indole excited state, and the electronic circular dichroism spectrum of an oxacepham. © 2014 Wiley Periodicals, Inc.     

A transferable coarse-grained model for hydrogen-bonding liquids

We present here a recent development of a generalized coarse-grained model for use in molecular simulations. In this model, interactions between coarse-grained particles consist of both van der Waals and explicit electrostatic components. As a result, the coarse-grained model offers the transferability that is lacked by most current effective-potential based approaches. The previous center-of-mass framework (P. A. Golubkov and P. Ren, J. Chem. Phys., 2006, 125, 64103) is generalized here to include arbitrary off-center interaction sites for both Gay-Berne and multipoles. The new model has been applied to molecular dynamic simulations of neat methanol liquid. By placing a single point multipole at the oxygen atom rather than at the center of mass of methanol, there is a significant improvement in the ability to capture hydrogen-bonding. The critical issue of transferability of the coarse-grained model is verified on methanol-water mixtures, using parameters derived from neat liquids without any modification. The mixture density and internal energy from coarse-grained molecular dynamics simulations show good agreement with experimental measurements, on a par with what has been obtained from more detailed atomic models. By mapping the dynamics trajectory from the coarse-grained simulation into the all-atom counterpart, we are able to investigate atomic-level structure and interaction. Atomic radial distribution functions of neat methanol, neat water and mixtures compare favorably to experimental measurements. Furthermore, hydrogen-bonded 6- and 7-molecule chains of water and methanol observed in the mixture are in agreement with previous atomic simulations.     

Ion solvation thermodynamics from simulation with a polarizable force field

Thermodynamic measurements of the solvation of salts and electrolytes are relatively straightforward, but it is not possible to separate total solvation free energies into distinct cation and anion contributions without reference to an additional extrathermodynamic assumption. The present work attempts to resolve this difficulty using molecular dynamics simulations with the AMOEBA polarizable force field and perturbation techniques to directly compute absolute solvation free energies for potassium, sodium, and chloride ions in liquid water and formamide. Corresponding calculations are also performed with two widely used nonpolarizable force fields. The simulations with the polarizable force field accurately reproduce in vacuo quantum mechanical results, experimental ion-cluster solvation enthalpies, and experimental solvation free energies for whole salts, while the other force fields do not. The results indicate that calculations with a polarizable force field can capture the thermodynamics of ion solvation and that the solvation free energies of the individual ions differ by several kilocalories from commonly cited values.     

Static and dynamic properties of the interface between a polymer brush and a melt of identical chains

Molecular-dynamics simulations of a short-chain polymer melt between two brush-covered surfaces under shear have been performed. The end-grafted polymers which constitute the brush have the same chemical properties as the free chains in the melt and provide a soft deformable substrate. Polymer chains are described by a coarse-grained bead-spring model, which includes excluded volume and backbone connectivity of the chains. The grafting density of the brush layer offers a way of controlling the behavior of the surface without altering the molecular interactions. We perform equilibrium and nonequilibrium molecular-dynamics simulations at constant temperature and volume using the dissipative particle dynamics thermostat. The equilibrium density profiles and the behavior under shear are studied as well as the interdigitation of the melt into the brush, the orientation on different length scales (bond vectors, radius of gyration, and end-to-end vector) of free and grafted chains, and velocity profiles. The obtained boundary conditions and slip length show a rich behavior as a function of grafting density and shear velocity.