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.

Configuration of influenza hemagglutinin fusion peptide monomers and oligomers in membranes.

The 20 N-terminal residues of the HA2 subunit of influenza hemagglutinin (HA), known as the fusion peptide, play a crucial role in membrane fusion. Molecular dynamics simulations with implicit solvation are employed here to study the structure and orientation of the fusion peptide in membranes. As a monomer the alpha-helical peptide adopts a shallow, slightly tilted orientation along the lipid tail-head group interface. The average angle of the peptide with respect to membrane plane is 12.4 degrees . We find that the kinked structure proposed on the basis of NMR data is not stable in our model because of the high energy cost related to the membrane insertion of polar groups. Because hemagglutinin-mediated membrane fusion is promoted by low pH, we examined the effect of protonation of the Glu and Asp residues. The configurations of the protonated peptides were slightly deeper in the membrane but at similar angles. Finally, because HA is a trimer, we modeled helical fusion peptide trimers. We find that oligomerization affects the insertion depth of the peptide and its orientation with respect to the membrane: a trimer exhibits equally favorable configurations in which some or all of the helices in the bundle insert obliquely deep into the membrane.

Solvent size vs cohesive energy as the origin of hydrophobicity.

Two main physical explanations of hydrophobicity seem to be currently competing. The classical, intuitive view attributes it to the fact that interactions between water molecules are much stronger than those between water and nonpolar groups. The second, "heretical" view attributes it to the small size of the water molecule which increases the entropic cost of opening up a cavity to accommodate the solute. Here we examine the solvation of methane in water and in model liquids that lack one or more of water's properties and report a detailed decomposition of the solvation free energy, enthalpy, entropy, and heat capacity in these solvents. The results fully support the classical view. It is found that fluids with strong intermolecular interactions favor expulsion of methane to its pure phase or to CCl(4), whereas fluids with weak intermolecular interactions do not. However, the specific thermodynamic signature of the hydrophobic effect (entropy driven at room temperature with a large heat capacity change) is a result of the hydrogen-bonding structure of water.

On the unfolding of alpha-lytic protease and the role of the pro region.

Molecular dynamics simulations of alpha-lytic protease (alphaLP) alone and complexed with its pro region (PRO) are performed to understand the origin of its high unfolding (and folding) barrier when it is alone and how the pro region lowers this barrier. At room temperature, alphaLP exhibits lower dynamic fluctuations than alpha-chymotrypsin. Simulation of PRO alone led to reorientation of its N terminal helix and collapse to a more compact state. A model for the uncleaved proenzyme was built and found to be stable in the time scale of the simulations. Energetic analysis suggests that the origin of strain in the uncleaved proenzyme compared with the cleaved complex is in the intramolecular backbone electrostatic interactions of the cleaved strand. In high temperature simulations, the interaction of the long beta hairpin of the enzyme with the C terminal beta sheet of PRO is among the most stable in the complex and a likely "nucleation site" for folding. In the course of unfolding, the C terminal tail of PRO is sometimes observed to intervene between the long hairpin and the aspartate loop of the enzyme, perhaps thereby lowering the energy barrier for separation of the two hairpins. Tighter interactions at the interface between the enzyme and its pro region are also occasionally observed, providing an additional mechanism for unfolding catalysis. Simulations of a mutant enzyme where the buried ion pair residues R102 and D142 were replaced by W and L, respectively, did not display any distinguishable behavior compared with the wild type.

Effective energy functions for protein structure prediction.

Protein structure prediction, fold recognition, homology modeling and design rely mainly on statistical effective energy functions. Although the theoretical foundation of such functions is not clear, their usefulness has been demonstrated in many applications. Molecular mechanics force fields, particularly when augmented by implicit solvation models, provide physical effective energy functions that are beginning to play a role in this area.

Understanding beta-hairpin formation.

The kinetics of formation of protein structural motifs (e.g., alpha-helices and beta-hairpins) can provide information about the early events in protein folding. A recent study has used fluorescence measurements to monitor the folding thermodynamics and kinetics of a 16-residue beta-hairpin. In the present paper, we obtain the free energy surface and conformations involved in the folding of an atomistic model for the beta-hairpin from multicanonical Monte Carlo simulations. The results suggest that folding proceeds by a collapse that is downhill in free energy, followed by rearrangement to form a structure with part of the hydrophobic cluster; the hairpin hydrogen bonds propagate outwards in both directions from the partial cluster. Such a folding mechanism differs from the published interpretation of the experimental results, which is based on a helix-coil-type phenomenological model.

Effective energy function for proteins in solution.

A Gaussian solvent-exclusion model for the solvation free energy is developed. It is based on theoretical considerations and parametrized with experimental data. When combined with the CHARMM 19 polar hydrogen energy function, it provides an effective energy function (EEF1) for proteins in solution. The solvation model assumes that the solvation free energy of a protein molecule is a sum of group contributions, which are determined from values for small model compounds. For charged groups, the self-energy contribution is accounted for primarily by the exclusion model. Ionic side-chains are neutralized, and a distance-dependent dielectric constant is used to approximate the charge-charge interactions in solution. The resulting EEF1 is subjected to a number of tests. Molecular dynamics simulations at room temperature of several proteins in their native conformation are performed, and stable trajectories are obtained. The deviations from the experimental structures are similar to those observed in explicit water simulations. The calculated enthalpy of unfolding of a polyalanine helix is found to be in good agreement with experimental data. Results reported elsewhere show that EEF1 clearly distinguishes correctly from incorrectly folded proteins, both in static energy evaluations and in molecular dynamics simulations and that unfolding pathways obtained by high-temperature molecular dynamics simulations agree with those obtained by explicit water simulations. Thus, this energy function appears to provide a realistic first approximation to the effective energy hypersurface of proteins.

Discrimination of the native from misfolded protein models with an energy function including implicit solvation.

An essential requirement for theoretical protein structure prediction is an energy function that can discriminate the native from non-native protein conformations. To date most of the energy functions used for this purpose have been extracted from a statistical analysis of the protein structure database, without explicit reference to the physical interactions responsible for protein stability. The use of the statistical functions has been supported by the widespread belief that they are superior for such discrimination to physics-based energy functions. An effective energy function which combined the CHARMM vacuum potential with a Gaussian model for the solvation free energy is tested for its ability to discriminate the native structure of a protein from misfolded conformations; the results are compared with those obtained with the vacuum CHARMM potential. The test is performed on several sets of misfolded structures prepared by others, including sets of about 650 good decoys for six proteins, as well as on misfolded structures of chymotrypsin inhibitor 2. The vacuum CHARMM potential is successful in most cases when energy minimized conformations are considered, but fails when applied to structures relaxed by molecular dynamics. With the effective energy function the native state is always more stable than grossly misfolded conformations both in energy minimized and molecular dynamics-relaxed structures. The present results suggest that molecular mechanics (physics-based) energy functions, complemented by a simple model for the solvation free energy, should be tested for use in the inverse folding problem, and supports their use in studies of the effective energy surface of proteins in solution. Moreover, the study suggests that the belief in the superiority of statistical functions for these purposes may be ill founded.

Heat capacity and compactness of denatured proteins.

One of the striking results of protein thermodynamics is that the heat capacity change upon denaturation is large and positive. This change is generally ascribed to the exposure of non-polar groups to water on denaturation, in analogy to the large heat capacity change for the transfer of small non-polar molecules from hydrocarbons to water. Calculations of the heat capacity based on the exposed surface area of the completely unfolded denatured state give good agreement with experimental data. This result is difficult to reconcile with evidence that the heat denatured state in the absence of denaturants is reasonably compact. In this work, sample conformations for the denatured state of truncated CI2 are obtained by use of an effective energy function for proteins in solution. The energy function gives denatured conformations that are compact with radii of gyration that are slightly larger than that of the native state. The model is used to estimate the heat capacity, as well as that of the native state, at 300 and 350 K via finite enthalpy differences. The calculations show that the heat capacity of denaturation can have large positive contributions from non-covalent intraprotein interactions because these interactions change more with temperature in non-native conformations than in the native state. Including this contribution, which has been neglected in empirical surface area models, leads to heat capacities of unfolding for compact denatured states that are consistent with the experimental heat capacity data. Estimates of the stability curve of CI2 made with the effective energy function support the present model.

Protein sidechain conformer prediction: a test of the energy function.

BACKGROUND: Homology modeling is an important technique for making use of the rapidly increasing number of protein sequences in the absence of structural information. The major problems in such modeling, once the alignment has been made, concern the positions of loops and the orientations of sidechains. Although progress has been made in recent years for sidechain prediction, current methods appear to have a limit on the order of 70% in their accuracy. It is important to have an understanding of this limitation, which for energy-based methods could arise from inaccuracies of the potential function. RESULTS: A test of the CHARMM function for sidechain prediction was performed. To eliminate the multiple-residue search problem, the minimum energy positions of individual sidechains in ten proteins were calculated in the presence of all other sidechains in their crystal orientations. This test provides a necessary condition that any energy function useful for sidechain placement must satisfy. For chi1 x chi2 rotations, the accuracies were 77.4% and 89.5%, respectively, and in the presence of crystal waters were 86.5% and 94.9%, respectively. If there was an error, the crystal structure usually corresponded to an alternative local minimum on the calculated energy map. Prediction accuracy correlated with the size of the energy gap between primary and secondary minima. CONCLUSIONS: The results indicate that the errors in current sidechain prediction schemes cannot be attributed to the potential energy function per se. The test used here establishes a necessary condition that any proposed energy-based sidechain prediction method, as well as many statistically based methods, must satisfy.

"New view" of protein folding reconciled with the old through multiple unfolding simulations.

Twenty-four molecular dynamics trajectories of chymotrypsin inhibitor 2 provide a direct demonstration of the diversity of unfolding pathways. Comparison with experiments suggests that the transition state region for folding and unfolding occurs early with only 25 percent of the native contacts and that the root-mean-square deviations between contributing structures can be as large as 15 angstroms. Nevertheless, a statistically preferred unfolding pathway emerges from the simulations; disruption of tertiary interactions between the helix and a two-stranded portion of the beta sheet is the primary unfolding event. The results suggest a synthesis of the "new" and the classical view of protein folding with a preferred pathway on a funnel-like average energy surface.

Dynamics and unfolding pathways of a hyperthermophilic and a mesophilic rubredoxin.

Molecular dynamics simulations in solution are performed for a rubredoxin from the hyperthermophilic archaeon Pyrococcus furiosus (RdPf) and one from the mesophilic organism Desulfovibrio vulgaris (RdDv). The two proteins are simulated at four temperatures: 300 K, 373 K, 473 K (two sets), and 500 K; the various simulations extended from 200 ps to 1,020 ps. At room temperature, the two proteins are stable, remain close to the crystal structure, and exhibit similar dynamic behavior; the RMS residue fluctuations are slightly smaller in the hyperthermophilic protein. An analysis of the average energy contributions in the two proteins is made; the results suggest that the intraprotein energy stabilizes RdPf relative to RdDv. At 373 K, the mesophilic protein unfolds rapidly (it begins to unfold at 300 ps), whereas the hyperthermophilic does not unfold over the simulation of 600 ps. This is in accord with the expected stability of the two proteins. At 473 K, where both proteins are expected to be unstable, unfolding behavior is observed within 200 ps and the mesophilic protein unfolds faster than the hyperthermophilic one. At 500 K, both proteins unfold; the hyperthermophilic protein does so faster than the mesophilic protein. The unfolding behavior for the two proteins is found to be very similar. Although the exact order of events differs from one trajectory to another, both proteins unfold first by opening of the loop region to expose the hydrophobic core. This is followed by unzipping of the beta-sheet. The results obtained in the simulation are discussed in terms of the factors involved in flexibility and thermostability.

The effect of loratadine on activated cells of the nasal mucosa in patients with allergic rhinitis.

The effect of loratadine on the numbers of activated cells--cells expressing interleukin-2 receptors(IL-2R), HLA-DR antigens and proliferating cell nuclear antigen (PCNA)--in the nasal mucosa was studied in 48 patients with allergic rhinitis. Patients were treated with either loratadine (10 mg) or placebo for 1 month. At the end of treatment, a significant decrease in the symptom scores was noted in both groups of patients. However, the clinical score was significantly lower in the loratadine group compared to the placebo group. At the end of treatment, the numbers of IL-2R+, HLA-DR+ and PCNA+ cells were significantly decreased only in the group on loratadine. An almost significant correlation was also observed between the numbers of IL-2R+ cells and symptoms in the loratadine group. Our results show that loratadine exerts its beneficial effect possibly by inhibiting both the action of histamine and immune activation.