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.

Understanding protein folding via free-energy surfaces from theory and experiment.

The ability of protein molecules to fold into their highly structured functional states is one of the most remarkable evolutionary achievements of biology. In recent years, our understanding of the way in which this complex self-assembly process takes place has increased dramatically. Much of the reason for this advance has been the development of energy surfaces (landscapes), which allow the folding reaction to be described and visualized in a meaningful manner. Analysis of these surfaces, derived from the constructive interplay between theory and experiment, has led to the development of a unified mechanism for folding and a recognition of the underlying factors that control the rates and products of the folding process.

Is protein unfolding the reverse of protein folding? A lattice simulation analysis.

Simulations and experiments that monitor protein unfolding under denaturing conditions are commonly employed to study the mechanism by which a protein folds to its native state in a physiological environment. Due to the differences in conditions and the complexity of the reaction, unfolding is not necessarily the reverse of folding. To assess the relevance of temperature initiated unfolding studies to the folding problem, we compare the folding and unfolding of a 125-residue protein model by Monte Carlo dynamics at two temperatures; the lower one corresponds to the range used in T -jump experiments and the higher one to the range used in unfolding simulations of all-atom models. The trajectories that lead from the native state to the denatured state at these elevated temperatures are less diverse than those observed in the folding simulations. At the lower temperature, the system unfolds through a mandatory intermediate that corresponds to a local free energy minimum. At the higher temperature, no such intermediate is observed, but a similar pathway is followed. The structures contributing to the unfolding pathways resemble most closely those that make up the "fast track" of folding. The transition state for unfolding at the lower temperature (above Tm) is determined and is found to be more structured than the transition state for folding below the melting temperature. This shift towards the native state is consistent with the Hammond postulate. The implications for unfolding simulations of higher resolution models and for unfolding experiments of proteins are discussed.

Defining cooperativity in gene regulation locally through intrinsic noise.

Regulatory networks in cells may comprise a variety of types of molecular interactions. The most basic are pairwise interactions, in which one species controls the behaviour of another (e.g. a transcription factor activates or represses a gene). Higher-order interactions, while more subtle, may be important for determining the function of networks. Here, the authors systematically expand a simple master equation model for a gene to derive an approach for robustly assessing the cooperativity (effective copy number) with which a transcription factor acts. The essential idea is that moments of a joint distribution of protein copy numbers determine the Hill coefficient of a cis-regulatory input function without non-linear fitting. The authors show that this method prescribes a definition of cooperativity that is meaningful even in highly complex situations in which the regulation does not conform to a simple Hill function. To illustrate the utility of the method, the authors measure the cooperativity of the transcription factor CI in simulations of phage- and show how the cooperativity accurately reflects the behaviour of the system. The authors numerically assess the effects of deviations from ideality, as well as possible sources of error. The relationship to other definitions of cooperativity and issues for experimentally realising the procedure are discussed.

A metastable state in folding simulations of a protein model.

The native state of a protein is generally believed to be the global free energy minimum. However, there is increasing evidence that kinetically selected states play a role in the biological function of some proteins. In a recent folding study of a 125-residue heteropolymer model, one of 200 sequences was found to fold repeatedly to a particular local minimum that did not interconvert to the global minimum. The kinetic preference for this 'metastable' state is shown to derive from an entropic barrier associated with inserting a tail segment into the protein interior of the serpin-like global minimum structure. The relation of the present results to the role of metastable states in functioning and pathogenic proteins is discussed.

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.

Factors that affect the folding ability of proteins.

The folding ability of a heteropolymer model for proteins subject to Monte Carlo dynamics on a simple cubic lattice is shown to be strongly correlated with the stability of the native state. We consider a number of estimates of the stability that can be determined without simulation, including the energy gap between the native state and the structurally dissimilar part of the spectrum (Z score) and, for sequences with fully compact native states, the gap in energy between the native and first excited fully compact states. These estimates are found to be more robust predictors of folding ability than a parameter sigma that requires simulation for its evaluation: sigma = 1 - Tf/Ttheta, where Tf is the temperature at which the fluctuation of an order parameter is at its maximum and Ttheta is the temperature at which the specific heat is at its maximum. We show that the interpretation of Ttheta as the collapse transition temperature is not correct in general and that the correlation between sigma and the folding ability arises from the fact that sigma is related to the energy gap (Z score).

Use of quantitative structure-property relationships to predict the folding ability of model proteins.

We investigate the folding of a 125-bead heteropolymer model for proteins subject to Monte Carlo dynamics on a simple cubic lattice. Detailed study of a few sequences revealed a folding mechanism consisting of a rapid collapse followed by a slow search for a stable core that served as the transition state for folding to a near-native intermediate. Rearrangement from the intermediate to the native state slowed folding further because it required breaking native-like local structure between surface monomers so that those residues could condense onto the core. We demonstrate here the generality of this mechanism by a statistical analysis of a 200 sequence database using a method that employs a genetic algorithm to pick the sequence attributes that are most important for folding and an artificial neural network to derive the corresponding functional dependence of folding ability on the chosen sequence attributes [quantitative structure-property relationships (QSPRs)]. QSPRs that use three sequence attributes yielded substantially more accurate predictions than those that use only one. The results suggest that efficient search for the core is dependent on both the native state's overall stability and its amount of kinetically accessible, cooperative structure, whereas rearrangement from the intermediate is facilitated by destabilization of contacts between surface monomers. Implications for folding and design are discussed.

The folding mechanism of larger model proteins: role of native structure.

The folding mechanism of a 125-bead heteropolymer model for proteins is investigated with Monte Carlo simulations on a cubic lattice. Sequences that do and do not fold in a reasonable time are compared. The overall folding behavior is found to be more complex than that of models for smaller proteins. Folding begins with a rapid collapse followed by a slow search through the semi-compact globule for a sequence-dependent stable core with about 30 out of 176 native contacts which serves as the transition state for folding to a near-native structure. Efficient search for the core is dependent on structural features of the native state. Sequences that fold have large amounts of stable, cooperative structure that is accessible through short-range initiation sites, such as those in anti-parallel sheets connected by turns. Before folding is completed, the system can encounter a second bottleneck, involving the condensation and rearrangement of surface residues. Overly stable local structure of the surface residues slows this stage of the folding process. The relation of the results from the 125-mer model studies to the folding of real proteins is discussed.

Use of a quantitative structure-property relationship to design larger model proteins that fold rapidly.

A quantitative structure-property relationship (QSPR) was used to design model protein sequences that fold repeatedly and relatively rapidly to stable target structures. The specific model was a 125-residue heteropolymer chain subject to Monte Carlo dynamics on a simple cubic lattice. The QSPR was derived from an analysis of a database of 200 sequences by a statistical method that uses a genetic algorithm to select the sequence attributes that are most important for folding and a neural network to determine the corresponding functional dependence of folding ability on the chosen attributes. The QSPR depends on the number of anti-parallel sheet contacts, the energy gap between the native state and quasi-continuous part of the spectrum and the total energy of the contacts between surface residues. Two Monte Carlo procedures were used in series to optimize both the target structures and the sequences. We generated 20 fully optimized sequences and 60 partially optimized control sequences and tested each for its ability to fold in dynamic MC simulations. Although sequences in which either the number of anti-parallel sheet contacts or the energy of the surface residues is non-optimal are capable of folding almost as well as fully optimized ones, sequences in which only the energy gap is optimized fold markedly more slowly. Implications of the results for the design of proteins are discussed.

Uracil-DNA glycosylase acts by substrate autocatalysis.

In humans, uracil appears in DNA at the rate of several hundred bases per cell each day as a result of misincorporation of deoxyuridine (dU) or deamination of cytosine. Four enzymes that catalyse the hydrolysis of the glycosylic bond of dU in DNA to yield an apyridiminic site as the first step in base excision repair have been identified in the human genome. The most efficient and well characterized of these uracil-DNA glycosylases is UDG (also known as UNG and present in almost all known organisms), which excises U from single- or double-stranded DNA and is associated with DNA replication forks. We used a hybrid quantum-mechanical/molecular-mechanical (QM/MM) approach to determine the mechanism of catalysis by UDG. In contrast to the concerted associative mechanism proposed initially, we show here that the reaction proceeds in a stepwise dissociative manner. Cleavage of the glycosylic bond yields an intermediate comprising an oxocarbenium cation and a uracilate anion. Subsequent attack by a water molecule and transfer of a proton to D145 result in the products. Surprisingly, the primary contribution to lowering the activation energy comes from the substrate, rather than from the enzyme. This 'autocatalysis' derives from the burial and positioning of four phosphate groups that stabilize the rate-determining transition state. The importance of these phosphates explains the residual activity observed for mutants that lack key residues. A corresponding catalytic mechanism could apply to the DNA glycosylases TDG and SMUG1, which belong to the same structural superfamily as UDG.

One of two unstructured domains of Ii becomes ordered in complexes with MHC class II molecules.

We studied the role of the invariant chain (Ii) protein's structure in its ability to form complexes with major histocompatibility complex class II molecules. Multidimensional nuclear magnetic resonance experiments demonstrated that Ii contains two unstructured, flexible domains: a 39 residue sequence that contains a region (CLIP) critical for Ii/class II complex formation and becomes rapidly ordered when Ii/class II complexes are assembled, and a 30 residue sequence that contains the insertion point for a protease inhibitor domain included in an alternative splice form of Ii. Mobility of these domains guarantees accessibility to CLIP and the inhibitor insert, and ordering of the CLIP-containing domain may provide protection against proteolysis and contribute, along with Ii's compact 118-192 domain, to allotype-independent class II binding.