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.

First-principles calculation of the folding free energy of a three-helix bundle protein.

The folding and unfolding of a three-helix bundle protein were explored with molecular-dynamics simulations, cluster analysis, and weighted-histogram techniques. The folding-unfolding process occurs by means of a "folding funnel," in which a uniform and broad distribution of conformational states is accessible outside of the native manifold. This distribution narrows near a transition region and becomes compact within the native manifold. Key thermodynamic steps in folding include initial interactions around the amino-terminal helix-turn-helix motif, interactions between helices I and II, and, finally, the docking of helix III onto the helix I-II subdomain. A metastable minimum in the calculated free-energy surface is observed at approximately 1.5 times the native volume. Folding-unfolding thermodynamics are dominated by the opposing influences of protein-solvent energy, which favors unfolding, and the overall entropy, which favors folding by means of the hydrophobic effect.

The p53--Mdm2--HAUSP complex is involved in p53 stabilization by HAUSP.

The ubiquitin-specific protease HAUSP is a critical component of the p53-Mdm2 pathway by acting as a specific deubiquitinase for both p53 and Mdm2. Recent structural studies have indicated that p53 and Mdm2 bind to the N-terminal TRAF-like domain of HAUSP in a mutually exclusive manner. To understand the mechanism of HAUSP-mediated effects, we have created a p53 mutant that lacks HAUSP binding based on the crystal structure analysis. Indeed, this mutant p53 protein can be degraded by Mdm2 but fails to interact with HAUSP both in vitro and in vivo. Surprisingly, however, we have found that direct interaction between HAUSP and p53 is not absolutely required for it to antagonize efficiently Mdm2-mediated ubiquitination of p53 and that HAUSP is capable of enzymatically functioning in trans on p53 by using Mdm2 as a bridge. Further, we show that a trimeric protein complex containing p53, Mdm2 and HAUSP can exist in vivo, despite mutually exclusive binding, with Mdm2 serving as a binding mediator for p53 and HAUSP. These findings reveal the complication of HAUSP-mediated effects in the p53-Mdm2 interplay. It also has important implications for the development of novel chemotherapeutic compounds aimed at blocking this protein-protein interaction.

A method for predicting protein structure from sequence.

BACKGROUND: The ability to predict the native conformation of a globular protein from its amino-acid sequence is an important unsolved problem of molecular biology. We have previously reported a method in which reduced representations of proteins are folded on a lattice by Monte Carlo simulation, using statistically-derived potentials. When applied to sequences designed to fold into four-helix bundles, this method generated predicted conformations closely resembling the real ones. RESULTS: We now report a hierarchical approach to protein-structure prediction, in which two cycles of the above-mentioned lattice method (the second on a finer lattice) are followed by a full-atom molecular dynamics simulation. The end product of the simulations is thus a full-atom representation of the predicted structure. The application of this procedure to the 60 residue, B domain of staphylococcal protein A predicts a three-helix bundle with a backbone root mean square (rms) deviation of 2.25-3 A from the experimentally determined structure. Further application to a designed, 120 residue monomeric protein, mROP, based on the dimeric ROP protein of Escherichia coli, predicts a left turning, four-helix bundle native state. Although the ultimate assessment of the quality of this prediction awaits the experimental determination of the mROP structure, a comparison of this structure with the set of equivalent residues in the ROP dime- crystal structure indicates that they have a rms deviation of approximately 3.6-4.2 A. CONCLUSION: Thus, for a set of helical proteins that have simple native topologies, the native folds of the proteins can be predicted with reasonable accuracy from their sequences alone. Our approach suggest a direction for future work addressing the protein-folding problem.

Methodological advances in molecular dynamics simulations of biological systems.

New advances in the techniques used to simulate specific statistical ensembles provide molecular dynamics algorithms that permit rigorous connections to be made between thermodynamic observables and calculated quantities in simulations of biological molecules. The complete inclusion of electrostatic forces in simulations also improves the comparison between calculations of simple structural measures in crystals and X-ray crystallographic results. These advances coupled with the ongoing development of more accurate inter/intramolecular forcefields with enhanced accuracy provide guidelines and benchmarks for comparison as we move to study more complicated biological problems.

Simulations of protein folding and unfolding.

The 'new view' of protein folding is based on a statistical analysis of the landscape for protein folding. This perspective leads to the investigation of the statistical distributions of protein conformations as the folding protein approaches the native state from unfolded states. Molecular dynamics simulations of both the thermodynamics of protein folding and the kinetics of unfolding are beginning to explore the statistical nature of these distributions. They also provide connections between the theory of protein folding landscapes and the experimental observations of the properties of key ensembles of the conformations populated as folding progresses.

A novel agent SL-401 induces anti-myeloma activity by targeting plasmacytoid dendritic cells, osteoclastogenesis and cancer stem-like cells.

Novel therapies for multiple myeloma (MM) can target mechanism(s) in the host-MM bone marrow (BM) microenvironment mediating MM progression and chemoresistance. Our studies showed increased numbers of tumor-promoting, immunosuppressive and drug-resistant plasmacytoid dendritic cells (pDCs) in the MM BM microenvironment. pDC-MM cell interactions upregulate interleukin-3 (IL-3), which stimulates both pDC survival and MM cell growth. Since IL-3 R is highly expressed on pDCs in the MM BM milieu, we here targeted pDCs using a novel IL-3 R-targeted therapeutic SL-401. In both in vitro and in vivo models of MM in its BM milieu, SL-401 decreases viability of pDCs, blocks pDC-induced MM cell growth, and synergistically enhances anti-MM activity of bortezomib and pomalidomide. Besides promoting pDC survival and MM cell growth, IL-3 also mediates progression of osteolytic bone disease in MM. Osteoclast (OCL) progenitor cells express IL-3 R, and we show that SL-401 abrogates monocyte-derived OCL formation and bone resorption. Finally, we show that SL-401 also decreases the viability of IL-3 R-expressing cancer stem-like cells in MM. Overall, our study provides the preclinical basis for clinical trials of SL-401 to block pDC-induced MM cell growth, inhibit osteoclastogenesis and target MM stem-like cell subpopulations to improve patient outcome in MM.

Effects of an S84E mutation of bovine growth hormone in transgenic mice.

The ability of mutant bovine growth hormones (bGH) to serve as either agonist or antagonist has been demonstrated in transgenic mice. We have prepared two transgenic strains of FVB/N mice, one expressing wild-type bGH and a second with a glutamic acid mutation at serine 84 in helix 2. Comparison of their phenotypes to those of nontransgenic littermates indicates that wild-type bGH induces a previously described phenotype for hyper-somatotrophic mice. In contrast, the replacement of the side chain hydroxyl at serine 84 with acetic acid produced a phenotype that expressed bGH at appreciable concentrations, but failed to elicit the phenotype observed with either an agonist or an antagonist of bGH. These results indicate that serine 84 is crucial for the activity of bGH despite this site being distal to the receptor binding surfaces.

Inhibition of mammary tumorigenesis in carcinogen-treated Lewis rats by suppression of prolactin secretion.

One hundred sixty-eight female Lewis rats were treated intragastrically with 10 mg 7,12-dimethylbenz[a]anthracene at 56 and 63 days of age. Pituitary prolactin secretion was suppressed in one-half of these rats by daily sc administrations of 2-bromoergocryptine mesylate (CB-154; 0.4 mg/100 g body wt) from 29 to 90 days of age (series 1) and from 90 to 140 days of age (series 2). Treatment with CB-154 was initiated prior to the onset of palpable mammary carcinomas. Control rats were given injections of saline. Inguinal mammary glands were excised from 10 control and 10 CB-154-treated rats at the cessation of saline and CB-154 treatments and examined for hyperplastic nodules (HN). The remaining rats were palpated weekly for mammary carcinomas (MC) and killed at 200 days of age. Mean number of HN per rat, mean number of MC per rat, and percent of rats with MC were, respectively: series 1--controls, 0.6, 1.5, and 68; CB-154 treatment, 0.5, 1.1, and 62; series 2--controls, 10.4, 2.0, and 94; CB-154 treatment, 5.1, 1.1, and 56. The number of HN and MC was only slightly reduced in rats when prolactin was suppressed during carcinogen treatment (series 1) but markedly reduced when prolactin was suppressed after carcinogen treatment (series 2). These results provide evidence that prolactin is involved in the early development of mammary dysplasias in the carcinogen-treated female Lewis rat.

Bone and kidney adenylate cyclase-stimulating activity produced by a hypercalcemic canine adenocarcinoma line (CAC-8) maintained in nude mice.

The tumor line CAC-8, is a serially transplantable adenocarcinoma maintained in nude mice which originated from a hypercalcemic dog. Nude mice with CAC-8 developed a syndrome of humoral hypercalcemia of malignancy. CAC-8 contained a protein factor which stimulated adenylate cyclase of bone and kidney cells in vitro. The adenylate cyclase (AC) of rat osteosarcoma cell lines, ROS 17/2.8 (ROS) and UMR-106, was stimulated by the tumor extract and potentiated by forskolin (0.1 microM). The ROS cells responded to the lowest concentration of CAC-8 extract, but UMR cells responded with a greater increase in AC activity compared to controls following exposure to CAC-8 extract. Pretreatment of ROS 17/2.8 cells with dexamethasone enhanced the response to CAC-8 extract. The opossum kidney cell line (OK) was less sensitive to the AC-stimulating activity of CAC-8 extract, but AC stimulation was increased in the presence of forskolin. Bovine (1-34) parathyroid hormone (BPTH) (10 nM) stimulated AC equally in ROS, UMR, and OK cells. Isoproterenol (1.0 microM) stimulated AC activity in ROS and UMR cells but not in OK cells. The AC-stimulating activity of CAC-8 appeared to bind to the parathyroid hormone receptor of ROS, UMR, and OK cells since addition of the parathyroid hormone receptor antagonist, [8,18norleucine, 34tyrosine]BPTH (3-34) amide, inhibited CAC-8-mediated cyclic adenosine 5'-monophosphate production and alone did not stimulate AC activity. The AC-stimulating activity of CAC-8 was acid and heat stable. Trypsin digestion reduced BPTH and CAC-8 stimulation of AC to near basal levels and treatment of CAC-8 extract with dithiothreitol reduced AC stimulation in UMR cells by approximately 50%. Extracts of the hypercalcemic tumor line (CAC-8) contained bone and kidney AC-stimulating activity which was enhanced by forskolin and dexamethasone, inhibited by [8,18Nle, 34Tyr]BPTH (3-34) amide, heat stable, trypsin sensitive, inactivated by reduction, and had a relative molecular weight of 34,000 by gel exclusion chromatography. Isolation and characterization of the factor(s) produced by CAC-8 that stimulate AC activity will be useful in further understanding the pathogenesis of humoral hypercalemia of malignancy in animal and human patients.

Specificity of inhibition of calcium- and calmodulin-dependent protein kinase by alloxan.

Studies were undertaken to determine whether the effect of alloxan to inactivate a membrane-bound calcium- and calmodulin-dependent protein kinase was specific for the pancreatic islets and whether inactivation of the kinase occurred also after injection of a diabetogenic dose of alloxan into rats. The effect of alloxan was also examined on similar particulate calcium- and calmodulin-dependent kinases present in two other secretory tissues, mammary acini and forebrain. Exposure of alloxan to cell-free preparations of all secretory tissues examined inhibited the calcium- and calmodulin-dependent kinase activities, suggesting that the specificity of alloxan action was not due to the presence in islets of a kinase uniquely sensitive to alloxan. To determine whether the selective effect of alloxan action was mediated at the cellular level, experiments were performed with alloxan presented to intact cells. Whereas alloxan exposure to viable cell preparations of islets and brain decreased the subsequently measured calcium- and calmodulin-dependent protein kinase activity, the activity measured in mammary acini exposed to these alloxan concentrations was unaffected. Injection (i.v.) of a diabetogenic dose of alloxan (50 mg/kg) produced an immediate (10 min) and selective inactivation of the calcium- and calmodulin-dependent protein kinase in pancreatic islets but had no effect on the similar kinases measured in brain and mammary acini. These results indicate that the unique sensitivity of islets to alloxan may result from the ability of alloxan to rapidly gain intracellular access and then inactivate this kinase activity. The selective effect of alloxan injection on this islet protein kinase is consistent with the hypothesis that inactivation of the kinase by alloxan is related to its diabetogenic effect in vivo.

Characterization of "native" apomyoglobin by molecular dynamics simulation.

We have used molecular dynamics simulation methods to study the structure and fluctuations of "native" apomyoglobin in aqueous solution for a period of greater than 0.5 nanosecond. This work was motivated by the recent attempts of Hughson et al. to characterize the structure and motion of both this molecule and the less compact, acid stabilized I stage, using methods of pulsed H/2H exchange. The study of these systems provides new insights into protein folding intermediates and our simulation has yielded a detailed model for structure and fluctuations in apomyoglobin which complements the experimental studies. We find that local (short-time) fluctuations agree well with fluctuations observed for the holoprotein in aqueous solution, as well as results from the crystallographic B-factors. In addition, the structural features we observe for native apomyoglobin are very similar to the holoprotein, in basic agreement with the findings of Hughson et al. By examining larger-scale motions, developing only over timescales in excess of a 100 picoseconds, we are able to identify conformationally "labile" and "non-labile" regions within native apomyoglobin. These regions correspond extremely well to those identified in the nuclear magnetic resonance experiments as unstable and stable "folding subdomains" in the I state of apomyoglobin. Overall we find that helices A, B, E, G and H show the least amount of motion and helices C, D and F move substantially over the timescales examined. The major motions, and the primary difference between the holo and apo structures as we have observed them, are due to the shifting motion of helices C, D and F into the vacant heme cavity. We also find that motions at the interface of helical segments can be large, with one important exception being the chain segment connecting helices G and H. This segment of chain interacts with the conformationally "non-labile" helix A to form a relatively rigid subdomain composed of helices A, G and H. We believe that these findings provide direct support for the suggestion of Hughson et al. that helices A, G and H constitute a compact subdomain that remains in a native-like conformation as the protein begins to unfold in environments of decreasing pH.

A microscopic view of helix propagation: N and C-terminal helix growth in alanine helices.

Molecular dynamics simulations with umbrella sampling are used to perform free energy simulations of C-terminal and N-terminal helix propagation in small helices of Ace-(Ala)n-NMe, with n= (4,5,10,15), in water. From the resulting free energy surfaces, computed as a function of the terminal psi dihedral angle, the roles of length and end effects in helix propagation are explored. An energetic analysis of the helices, both formed and partially formed, is used to develop a molecular rationalization for the observed trends in helix stability. We find that the microscopic helix propagation parameters vary significantly depending on the end and length of the helix in which the terminal hydrogen bond is forming. A model which considers propagation of the helices from either end as statistically independent yields Zimm-Bragg s parameters in the range of 0.5 to 1.5, depending on helical length. Analysis of the mechanism of helix propagation suggests that 3(10)-helix plays a role in helix formation but its population should be low in the helical state of these model peptides.

Calculations on folding of segment B1 of streptococcal protein G.

We present an investigation of the folding thermodynamics and mechanism of segment B1 of streptococcal protein G. Molecular dynamics simulations of the fully solvated protein are used to probe thermodynamically significant states at different stages of folding. We performed several unfolding simulations to generate a database of initial conditions. The database is analyzed and clustered. The cluster centers extracted from this database were then used as starting points for umbrella sampling of the folding free energy landscape under folding conditions. The resulting sampling was combined with the weighted histogram analysis method. One and two-dimensional free energy surfaces were constructed along several order parameters and used to analyze the folding process. Our findings indicate that an initial collapse precedes the formation of significant native structure. Elements of local structure originate in the regions of the protein shown to have higher H/2H exchange protection factors in early stages of folding. A non-native contact, observed experimentally at the N terminus of the alpha-helix in a peptide excised from the protein, is seen to pre-organize the chain in early stages of folding. Collapse and early structure formation yields a compact globule with a significant number of water molecules present. Desolvation of the protein core is coincident with the final stages of folding from the compact state.

Solvent effects on protein motion and protein effects on solvent motion. Dynamics of the active site region of lysozyme.

The stochastic boundary molecular dynamics methodology is applied to the active site of the enzyme lysozyme. A comparison is made of in vacuo dynamics results from the stochastic boundary method and a full conventional molecular dynamics simulation of lysozyme. Excellent agreement between the two approaches is obtained. The influence of solvent on the residues in the active site region is explored and it is shown that both the structure and dynamics are affected. Of particular importance for the structure of the protein is the solvation of polar residues and the stabilization of like-charged ion pairs. The magnitude of the fluctuations is only slightly altered by the solvent; the overall increase in the root-mean-square fluctuations, relative to the vacuum run, is 11%. The solvent effect on dynamical properties is found not to be simply related to the solvent viscosity. Both the solvent exposure and dynamic aspects of protein-solvent interactions, including the relative time scales of the motions, are shown to play a role. The effects of the protein on solvent dynamics and structure are also observed to be significant. The solvent molecules around atoms in charged, polar and apolar side-chains show markedly different diffusion coefficients as well as exhibiting different solvation structures. One key example is the water around apolar groups, which is much less mobile than bulk water, or water solvating polar groups.

Helicity, circular dichroism and molecular dynamics of proteins.

In protein unfolding studies, reduction in circular dichroism (CD) at 222 nm has been interpreted as loss of helicity. Estimates of the helicity of a protein from its CD spectrum are calibrated by reference to X-ray crystal structures, based on the assumption that the mean residue ellipticity at 222 nm is directly proportional to the number of residues in a helix. We have examined various influences on the CD at 222 nm, using molecular dynamics simulations to provide the structural detail required. We have found that the fragmentation of long helices, without a reduction in the number of helical residues, significantly reduces the mean residue ellipticity at 222 nm. The dynamical motion of the protein and the precise conformation of helical residues also play an important role. We discuss the implications of these factors to the interpretation of CD for the partial unfolding of apomyoglobin.

Thermodynamics of amide hydrogen bond formation in polar and apolar solvents.

We present the initial findings of a theoretical study of hydrogen bond formation between two formamide molecules in water and in carbon tetrachloride. These systems were chosen as the simplest models for secondary structure formation in the polar environment near the protein surface and the apolar environment of the protein interior. We have employed thermodynamic simulation methods to obtain absolute binding free energies and free energy profiles for the formation of peptide hydrogen bonds in the two solvents. We find that the amide hydrogen bond is stable by 8.4 kcal/mol in CCl4, and by 0.3 kcal/mol in water. Our results indicate also that the hydrogen-bonded dimer is 2.2 kcal/mol more stable in water than it is in CCl4. We compare our results with those from experiment, and discuss their use in interpreting mechanisms of protein folding.

Stability of a model beta-sheet in water.

We have used molecular dynamics simulations to determine the stability in water of a model beta-sheet formed by two alanine dipeptide molecules with two intermolecular hydrogen bonds in the closely spaced antiparallel arrangement. In this paper we describe our computations of the binding free energy of the model sheet and a portion of the free energy surface as a function of a reaction co-ordinate for sheet formation. We used the free energy surface to identify stable conformations along the reaction co-ordinate. To determine whether or not the model sheet with two hydrogen bonds is more stable than a single amide hydrogen bond in water, we compared the results of the present calculations to results from our earlier study of linear hydrogen bond formation between two formamide molecules (the formamide "dimer"). The free energy surfaces for the sheet and formamide dimer each have two minima corresponding to locally stable hydrogen-bonded and solvent-separated configurations. The binding free energies of the model sheet and the formamide dimer are -5.5 and -0.34 kcal/mol, respectively. Thus, the model sheet with two hydrogen bonds is quite stable while the simple amide hydrogen bond is only marginally stable. To understand the relative stabilities of the model sheet and formamide dimer in terms of solute-solute and solute-water interactions, we decomposed the free energy differences between hydrogen-bonded and solvent-separated conformations into energetic and entropic contributions. The changes in the peptide-peptide energy and the entropy are roughly twice as large for the sheet as they are for the formamide dimer. The magnitude of the peptide-water energy difference for the sheet is less than twice (by about 3.5 kcal/mol) that for the formamide dimer, and this accounts for the stability of the sheet. The presence of the side-chains and/or blocking groups apparently prevents the amide groups in the sheet from being solvated as favorably in the separated arrangement as in the formamide dimer, where the amide groups are completely exposed to the solvent.

Prediction of quaternary structure of coiled coils. Application to mutants of the GCN4 leucine zipper.

Using a simplified protein model, the equilibrium between different oligomeric species of the wild-type GCN4 leucine zipper and seven of its mutants have been predicted. Over the entire experimental concentration range, agreement with experiment is found in five cases, while in two cases agreement is found over a portion of the concentration range. These studies demonstrate a methodology for predicting coiled coil quaternary structure and allow for the dissection of the interactions responsible for the global fold. In agreement with the conclusion of Harbury et al., the results of the simulations indicate that the pattern of hydrophobic and hydrophilic residues alone is insufficient to define a protein's three-dimensional structure. In addition, these simulations indicate that the degree of chain association is determined by the balance between specific side-chain packing preferences and the entropy reduction associated with side-chain burial in higher-order multimers.

Reverse turns in blocked dipeptides are intrinsically unstable in water.

We have carried out molecular dynamics simulations to study the conformational equilibria of two blocked dipeptides, Ac-Ala-Ala-NHMe and trans-Ac-Pro-Ala-NHMe, in water (Ac, amino-terminal blocking group COCH3; NHMe, carboxy-terminal blocking group NHCH3). Using specialized sampling techniques we computed free-energy surfaces as functions of a conformation co-ordinate that corresponds to hydrogen-bonded reverse turns at small values and to extended conformations at large values. The free-energy difference between hydrogen-bonded reverse turn conformations and extended conformations, determined from the equilibrium constants for reverse turn unfolding, is approximately -5 kcal/mole for Ac-Ala-Ala-NHMe, and -10 kcal/mole for Ac-Pro-Ala-NHMe. These results demonstrate that reverse turns in blocked dipeptides are intrinsically unstable in water. That is, in the absence of strongly stabilizing sequence-specific inter-residue interactions involving side-chains and/or charged terminal groups, the extended conformations of small peptides are highly favored in solution. By thermodynamically decomposing the free-energy differences, we found that the peptide-water entropy is the primary reason for the exceptional stability of the extended conformations of both peptides, and that the differences between the two peptides are primarily due to differences in the peptide-water interactions. In addition, we assessed the "proline effect" on the conformational equilibria by comparing the differences in configurational entropies between the reverse turn and extended conformations of the two peptides. As expected, the extended conformation of the Pro-Ala peptide is destabilized by reduced configurational entropy, but the effect is negligible in the blocked dipeptides. Finally, we compared our results with the results of several other experimental studies to identify some of the specific interactions that may be responsible for stabilizing reverse turns in small peptides in solution.