Richard M. Venable

CHARMM: The biomolecular simulation program

CHARMM (Chemistry at HARvard Molecular Mechanics) is a highly versatile and widely used molecular simulation program. It has been developed over the last three decades with a primary focus on molecules of biological interest, including proteins, peptides, lipids, nucleic acids, carbohydrates, and small molecule ligands, as they occur in solution, crystals, and membrane environments. For the study of such systems, the program provides a large suite of computational tools that include numerous conformational and path sampling methods, free energy estimators, molecular minimization, dynamics, and analysis techniques, and model‐building capabilities. The CHARMM program is applicable to problems involving a much broader class of many‐particle systems. Calculations with CHARMM can be performed using a number of different energy functions and models, from mixed quantum mechanical‐molecular mechanical force fields, to all‐atom classical potential energy functions with explicit solvent and various boundary conditions, to implicit solvent and membrane models. The program has been ported to numerous platforms in both serial and parallel architectures. This article provides an overview of the program as it exists today with an emphasis on developments since the publication of the original CHARMM article in 1983.

Update of the CHARMM all-atom additive force field for lipids: Validation on six lipid types

A significant modification to the additive all-atom CHARMM lipid force field (FF) is developed and applied to phospholipid bilayers with both choline and ethanolamine containing head groups and with both saturated and unsaturated aliphatic chains. Motivated by the current CHARMM lipid FF (C27 and C27r) systematically yielding values of the surface area per lipid that are smaller than experimental estimates and gel-like structures of bilayers well above the gel transition temperature, selected torsional, Lennard-Jones and partial atomic charge parameters were modified by targeting both quantum mechanical (QM) and experimental data. QM calculations ranging from high-level ab initio calculations on small molecules to semiempirical QM studies on a 1,2-dipalmitoyl-sn-phosphatidylcholine (DPPC) bilayer in combination with experimental thermodynamic data were used as target data for parameter optimization. These changes were tested with simulations of pure bilayers at high hydration of the following six lipids: DPPC, 1,2-dimyristoyl-sn-phosphatidylcholine (DMPC), 1,2-dilauroyl-sn-phosphatidylcholine (DLPC), 1-palmitoyl-2-oleoyl-sn-phosphatidylcholine (POPC), 1,2-dioleoyl-sn-phosphatidylcholine (DOPC), and 1-palmitoyl-2-oleoyl-sn-phosphatidylethanolamine (POPE); simulations of a low hydration DOPC bilayer were also performed. Agreement with experimental surface area is on average within 2%, and the density profiles agree well with neutron and X-ray diffraction experiments. NMR deuterium order parameters (S(CD)) are well predicted with the new FF, including proper splitting of the S(CD) for the aliphatic carbon adjacent to the carbonyl for DPPC, POPE, and POPC bilayers. The area compressibility modulus and frequency dependence of (13)C NMR relaxation rates of DPPC and the water distribution of low hydration DOPC bilayers also agree well with experiment. Accordingly, the presented lipid FF, referred to as C36, allows for molecular dynamics simulations to be run in the tensionless ensemble (NPT), and is anticipated to be of utility for simulations of pure lipid systems as well as heterogeneous systems including membrane proteins.

CHARMM-GUI Membrane Builder toward realistic biological membrane simulations.

CHARMM‐GUI Membrane Builder, http://www.charmm‐gui.org/input/membrane, is a web‐based user interface designed to interactively build all‐atom protein/membrane or membrane‐only systems for molecular dynamics simulations through an automated optimized process. In this work, we describe the new features and major improvements in Membrane Builder that allow users to robustly build realistic biological membrane systems, including (1) addition of new lipid types, such as phosphoinositides, cardiolipin (CL), sphingolipids, bacterial lipids, and ergosterol, yielding more than 180 lipid types, (2) enhanced building procedure for lipid packing around protein, (3) reliable algorithm to detect lipid tail penetration to ring structures and protein surface, (4) distance‐based algorithm for faster initial ion displacement, (5) CHARMM inputs for P21 image transformation, and (6) NAMD equilibration and production inputs. The robustness of these new features is illustrated by building and simulating a membrane model of the polar and septal regions of E. coli membrane, which contains five lipid types: CL lipids with two types of acyl chains and phosphatidylethanolamine lipids with three types of acyl chains. It is our hope that CHARMM‐GUI Membrane Builder becomes a useful tool for simulation studies to better understand the structure and dynamics of proteins and lipids in realistic biological membrane environments.

Additive empirical force field for hexopyranose monosaccharides

We present an all‐atom additive empirical force field for the hexopyranose monosaccharide form of glucose and its diastereomers allose, altrose, galactose, gulose, idose, mannose, and talose. The model is developed to be consistent with the CHARMM all‐atom biomolecular force fields, and the same parameters are used for all diastereomers, including both the α‐ and β‐anomers of each monosaccharide. The force field is developed in a hierarchical manner and reproduces the gas‐phase and condensed‐phase properties of small‐molecule model compounds corresponding to fragments of pyranose monosaccharides. The resultant parameters are transferred to the full pyranose monosaccharides, and additional parameter development is done to achieve a complete hexopyranose monosaccharide force field. Parametrization target data include vibrational frequencies, crystal geometries, solute–water interaction energies, molecular volumes, heats of vaporization, and conformational energies, including those for over 1800 monosaccharide conformations at the MP2/cc‐pVTZ//MP2/6‐31G(d) level of theory. Although not targeted during parametrization, free energies of aqueous solvation for the model compounds compare favorably with experimental values. Also well‐reproduced are monosaccharide crystal unit cell dimensions and ring pucker, densities of concentrated aqueous glucose systems, and the thermodynamic and dynamic properties of the exocyclic torsion in dilute aqueous systems. The new parameter set expands the CHARMM additive force field to allow for simulation of heterogeneous systems that include hexopyranose monosaccharides in addition to proteins, nucleic acids, and lipids.

Molecular Dynamics Studies of Polyethylene Oxide and Polyethylene Glycol: Hydrodynamic Radius and Shape Anisotropy

A revision (C35r) to the CHARMM ether force field is shown to reproduce experimentally observed conformational populations of dimethoxyethane. Molecular dynamics simulations of 9, 18, 27, and 36-mers of polyethylene oxide (PEO) and 27-mers of polyethylene glycol (PEG) in water based on C35r yield a persistence length lambda = 3.7 A, in quantitative agreement with experimentally obtained values of 3.7 A for PEO and 3.8 A for PEG; agreement with experimental values for hydrodynamic radii of comparably sized PEG is also excellent. The exponent upsilon relating the radius of gyration and molecular weight (R(g) proportional, variantM(w)(upsilon)) of PEO from the simulations equals 0.515 +/- 0.023, consistent with experimental observations that low molecular weight PEG behaves as an ideal chain. The shape anisotropy of hydrated PEO is 2.59:1.44:1.00. The dimension of the middle length for each of the polymers nearly equals the hydrodynamic radius R(h)obtained from diffusion measurements in solution. This explains the correspondence of R(h) and R(p), the pore radius of membrane channels: a polymer such as PEG diffuses with its long axis parallel to the membrane channel, and passes through the channel without substantial distortion.

Additive and classical drude polarizable force fields for linear and cyclic ethers

Empirical force field parameters consistent with the CHARMM additive and classical Drude based polarizable force fields are presented for linear and cyclic ethers. Initiation of the optimization process involved validation of the aliphatic parameters based on linear alkanes and cyclic alkanes. Results showed the transfer to cyclohexane to yield satisfactory agreement with target data; however, in the case of cyclopentane direct transfer of the Lennard-Jones parameters was not sufficient due to ring strain, requiring additional optimization of these parameters for this molecule. Parameters for the ethers were then developed starting with the available aliphatic parameters, with the nonbond parameters for the oxygens optimized to reproduce both gas- and condensed-phase properties. Nonbond parameters for the polarizable model include the use of an anisotropic electrostatic model on the oxygens. Parameter optimization emphasized the development of transferable parameters between the ethers of a given class. The ether models are shown to be in satisfactory agreement with both pure solvent and aqueous solvation properties, and the resulting parameters are transferable to test molecules. The presented force field will allow for simulation studies of ethers in condensed phase and provides a basis for ongoing developments in both additive and polarizable force fields for biological molecules.

Brownian dynamics simulation of a lipid chain in a membrane bilayer

A Brownian dynamics simulation of a lipid chain is used to model the motional properties of a dipalmitoyl phosphatidylcholine bilayer. The effects of the bilayer environment on the chain are represented by a mean field derived from an extension of the Marcelja model. The simulation was run 44 million steps, the equivalent of approximately 0.66 μs for a viscosity of 2.2 cp. The results are compared with those of a 30 million step simulation of the chain in the absence of the mean field. Deuterium order parameters for the methylene groups along the chain and the average chain length calculated from the mean field trajectory are shown to converge to the experimentally determined values for DPPC with an appropriate choice of parameters. An analysis of the torsional dynamics of the chain, including transition rates and kink probabilities, is carried out. It is demonstrated that kink formation is sometimes, though not always, concerted. A comparison of the membrane and free chain simulations implies that the in...

Structural Basis of Peroxide-mediated Changes in Human Hemoglobin A NOVEL OXIDATIVE PATHWAY

Hydrogen peroxide (H2O2) triggers a redox cycle between ferric and ferryl hemoglobin (Hb) leading to the formation of a transient protein radical and a covalent hemeprotein cross-link. Addition of H2O2 to highly purified human hemoglobin (HbA0) induced structural changes that primarily resided within β subunits followed by the internalization of the heme moiety within α subunits. These modifications were observed when an equal molar concentration of H2O2 was added to HbA0 yet became more abundant with greater concentrations of H2O2. Mass spectrometric and amino acid analysis revealed for the first time that βCys-93 and βCys-112 were oxidized extensively and irreversibly to cysteic acid when HbA0 was treated with H2O2. Oxidation of further amino acids in HbA0 exclusive to the β-globin chain included modification of βTrp-15 to oxyindolyl and kynureninyl products as well as βMet-55 to methionine sulfoxide. These findings may therefore explain the premature collapse of the β subunits as a result of the H2O2 attack. Analysis of a tryptic digest of the main reversed phase-high pressure liquid chromatography fraction revealed two α-peptide fragments (α128 - α139) and a heme moiety with the loss of iron, cross-linked between αSer-138 and the porphyrin ring. The novel oxidative pathway of HbA0 modification detailed here may explain the diverse oxidative, toxic, and potentially immunogenic effects associated with the release of hemoglobin from red blood cells during hemolytic diseases and/or when cell-free Hb is used as a blood substitute.

A simulation based model of NMR T1 relaxation in lipid bilayer vesicles

The results of the Brownian dynamics simulation of a hydrocarbon chain in a membrane bilayer described in the preceding paper are used to analyze the 13C NMR T1 relaxation in lipid bilayer vesicles. The analysis shows that the frequency dependence of the relaxation does not arise from gauche–trans isomerization or from axial rotation of the entire lipid molecule. However, a model in which fast axial rotation (D∥≊2×1010 s−1) and slow noncollective diffusive director fluctuations (D⊥≊1–2×108 s−1) are superimposed on the internal motions quantitatively accounts for both the magnitude and frequency dependence of the T1 data. An effective viscosity for the interior of the bilayer in the range of 1 cp, and a director order parameter of 0.5–0.7 are required to fit the NMR data. Collective effects do not appear necessary for explaining the NMR T1 data in vesicles, although they may be important for multilamellar dispersions.

Molecular dynamics simulations of gel (LβI) phase lipid bilayers in constant pressure and constant surface area ensembles

The results of a series of molecular dynamics simulations of the gel state of a dipalmitoylphosphatidylcholine bilayer at 293 K are described. The simulations, ranging from 40 ps to 2.5 ns, show clearly that: a flexible cell geometry is essential during equilibration; Ewald summation of electrostatics is superior to spherical cutoff methods; water exchange with the carbonyl group of chain 2 takes place on the ns time scale, while there is almost no hydration of chain 1. There is overall good agreement (D-spacing, chain tilt, fraction gauche, and area compressibility modulus) with experiment, though the surface area per lipid is slightly underestimated. The randomization of torsion 1 of chain 2 from exclusively gauche minus (as specified in the initial condition modeled from the crystal structure of a related lipid) to a mixture of g+/g− over the course of approximately 2 ns is a critical feature of the study. The torsional equilibration proceeded steadily when simulating at constant surface tension, but w...