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.

Investigating the high affinity and low sequence specificity of calmodulin binding to its targets.

Calmodulin (CaM) is a calcium binding protein that regulates a wide range of enzymes. Recently the structures of a number of complexes between CaM and synthetic target peptides have been determined. The peptides correspond to the CaM-binding domain of skeletal and smooth muscle myosin light-chain kinase (MLCK) and calmodulin-dependent protein kinase II alpha. Comparison of the peptide-free and peptide-bound structures reveals that CaM undergoes a large conformational change when forming a complex, resulting in the formation of a binding surface that provides for an optimal interaction with its target. In this work, the available co-ordinates of the NMR solution structure of CaM-skeletal MLCK peptide are used as a basis upon which several molecular models of binding are built. The detailed features of the protein's peptide binding surface are revealed through two-dimensional topographical projections. Negatively charged margins at the binding surface extremities interact strongly with basic peptide residues separated by nine or ten positions. The binding surface core is hydrophobic and displays a groove with four deep pockets, which can accommodate bulky peptide residues at relative positions 4 and 8 (pocket A), 11 (pocket B), 13 (pocket C), 14 and 17 (pocket D). Therefore, both electrostatic and van der Waals' features contribute to the high affinity binding. A search for alternative peptide placements in the binding tunnel reveals the dominant role of specific electrostatic interactions in the binding energy. Apolar interactions are more permissive, such that the hydrophobic side-chains that line the binding tunnel adapt in order to maintain favourable van der Waals' contacts. The model suggests that the structure can accommodate large peptide translations (up to 5 A) and a reversed peptide binding mode, with a little loss in binding interaction energy. These calculations are compared with available experimental data, providing a structural rationale for the low sequence specificity of the CaM target recognition.

Domain motions in dihydrofolate reductase: a molecular dynamics study.

Molecular dynamics simulations have been carried out on the enzyme dihydrofolate reductase from Lactobacillus casei complexed with methotrexate, NADPH and 264 crystallographic water molecules. Analysis of correlations in atomic fluctuations reveal the presence of highly correlated motion (correlation coefficient > 0.6) in the region between residues 30 to 35 and 85 to 90 leading to the identification of two domains, an "adenosine-binding domain" and a "large domain", which rotate by 3 to 4 degrees with respect to each other. The strongest correlation (> 0.6) within the large domain involves a coupling between the motions of the "teen-loop", and the spatially contiguous loops linking beta 6-beta 7 and beta 7-beta 8. Moreover, there is a significant correlation (approximately 0.5) between the adenosine fragment of NADPH and the pteridine and p-aminobenzoyl fragments of methotrexate, which are separated by approximately 17 A, and is lost on removal of "rigid-body" motion from the original trajectory. This provides support for the idea that the relative motion of the two domains is a means by which the occupation of the binding site for the adenosine end of the coenzyme can affect methotrexate binding and vice versa. Quasiharmonic vibrational analysis of the trajectory reveals that the overall dynamics of the system are governed by domain motions whose contributions are dominant at low frequencies. In addition, different low-frequency modes are responsible for separately coupling the adenosine-binding site and parts of methotrexate.

Understanding the mechanism of ice binding by type III antifreeze proteins.

Type III antifreeze proteins (AFPs) are present in the body fluids of some polar fishes where they inhibit ice growth at subzero temperatures. Previous studies of the structure of type III AFP by NMR and X-ray identified a remarkably flat surface on the protein containing amino acids that were demonstrated to be important for interaction with ice by mutational studies. It was proposed that this protein surface binds onto the (1 0 [\bar 1] 0) plane of ice with the key amino acids interacting directly with the water molecules in the ice crystal. Here, we show that the mechanism of type III AFP interaction with ice crystals is more complex than that proposed previously. We report a high-resolution X-ray structure of type III AFP refined at 1.15 A resolution with individual anisotropic temperature factors. We report the results of ice-etching experiments that show a broad surface coverage, suggesting that type III AFP binds to a set of planes that are parallel with or inclined at a small angle to the crystallographic c-axis of the ice crystal. Our modelling studies, performed with the refined structure, confirm that type III AFP can make energetically favourable interactions with several ice surfaces.

Locally accessible conformations of proteins: multiple molecular dynamics simulations of crambin.

Multiple molecular dynamics (MD) simulations of crambin with different initial atomic velocities are used to sample conformations in the vicinity of the native structure. Individual trajectories of length up to 5 ns sample only a fraction of the conformational distribution generated by ten independent 120 ps trajectories at 300 K. The backbone atom conformational space distribution is analyzed using principal components analysis (PCA). Four different major conformational regions are found. In general, a trajectory samples only one region and few transitions between the regions are observed. Consequently, the averages of structural and dynamic properties over the ten trajectories differ significantly from those obtained from individual trajectories. The nature of the conformational sampling has important consequences for the utilization of MD simulations for a wide range of problems, such as comparisons with X-ray or NMR data. The overall average structure is significantly closer to the X-ray structure than any of the individual trajectory average structures. The high frequency (less than 10 ps) atomic fluctuations from the ten trajectories tend to be similar, but the lower frequency (100 ps) motions are different. To improve conformational sampling in molecular dynamics simulations of proteins, as in nucleic acids, multiple trajectories with different initial conditions should be used rather than a single long trajectory.

Conformational change in the activation of lipase: an analysis in terms of low-frequency normal modes.

The interfacial activation of Rhizomucor miehei lipase (RmL) involves the motion of an alpha-helical region (residues 82-96) which acts as a "lid" over the active site of the enzyme, undergoing a displacement from a "closed" to an "open" conformation upon binding of substrate. Normal mode analyses performed in both low and high dielectric media reveal that low-frequency vibrational modes contribute significantly to the conformational transition between the closed and open conformations. In these modes, the lid displacement is coupled to local motions of active site loops as well as global breathing motions. Atomic fluctuations of the first hinge of the lid (residues 83-84) are substantially larger in the low dielectric medium than in the high dielectric medium. Our results also suggest that electrostatic interactions of Arg86 play an important role in terms of both the intrinsic stability of the lid and its displacement, through enhancement of hinge mobility in a high dielectric medium. Additional calculations demonstrate that the observed patterns of atomic fluctuations are an intrinsic feature of the protein structure and not dependent on the nature of specific energy minima.

Structural consequences of the B5 histidine --> tyrosine mutation in human insulin characterized by X-ray crystallography and conformational analysis.

The addition of phenols to hexameric insulin solutions produces a particularly stable hexamer, resulting from a rearrangement in which residues B1-B8 change from an extended conformation (T-state) to form an alpha-helix (R-state). The R-state is, in part, stabilized by nonpolar interactions between the phenolic molecule and residue B5 His at the dimer-dimer interface. The B5 His --> Tyr mutant human insulin was constructed to see if the tyrosine side chain would mimic the effect of phenol binding in the hexamer and induce the R-state. In partial support of this hypothesis, the molecule crystallized as a half-helical hexamer (T(3)R(3)) in conditions that conventionally promote the fully nonhelical (T6) form. As expected, in the presence of phenol or resorcinol, the B5 Tyr hexamers adopt the fully helical (R6) conformation. Molecular modeling calculations were performed to investigate the conformational preference of the T-state B5 Tyr side chain in the T(3)R(3) form, this side chain being associated with structural perturbations of the A7-A10 loop in an adjacent hexamer. For an isolated dimer, several different orientations of the side chain were found, which were close in energy and readily interconvertible. In the crystal environment only one of these conformations remains low in energy; this conformation corresponds to that observed in the crystal structure. This suggests that packing constraints around residue B5 Tyr result in the observed structural rearrangements. Thus, rather than promoting the R-state in a manner analogous to phenol, the mutation appears to destabilize the T-state. These studies highlight the role of B5 His in determining hexamer conformation and in mediating crystal packing interactions, properties that are likely be important in vivo.

Locating interaction sites on proteins: the crystal structure of thermolysin soaked in 2% to 100% isopropanol.

Multiple-solvent crystal structure determination (MSCS) allows the position and orientation of bound solvent fragments to be identified by determining the structure of protein crystals soaked in organic solvents. We have extended this technique by the determination of high-resolution crystal structures of thermolysin (TLN), generated from crystals soaked in 2% to 100% isopropanol. The procedure causes only minor changes to the conformation of the protein, and an increasing number of isopropanol interaction sites could be identified as the solvent concentration is increased. Isopropanol occupies all four of the main subsites in the active site, although this was only observed at very high concentrations of isopropanol for three of the four subsites. Analysis of the isopropanol positions shows little correlation with interaction energy computed using a molecular mechanics force field, but the experimentally determined positions of isopropanol are consistent with the structures of known protein-ligand complexes of TLN.

Demonstration that 1-trans-epoxysuccinyl-L-leucylamido-(4-guanidino) butane (E-64) is one of the most effective low Mr inhibitors of trypsin-catalysed hydrolysis. Characterization by kinetic analysis and by energy minimization and molecular dynamics simulation of the E-64-beta-trypsin complex.

1-trans-Epoxysuccinyl-L-leucylamido(4-guanidino)butane (E-64) was shown to inhibit beta-trypsin by a reversible competitive mechanism; this contrasts with the widely held view that E-64 is a class-specific inhibitor of the cysteine proteinases and reports in the literature that it does not inhibit a number of other enzymes including, notably, trypsin. The K1, value (3 x 10(-5) M) determined by kinetic analysis of the hydrolysis of N alpha-benzoyl-L-arginine 4-nitroanilide in Tris/HCl buffer, pH 7.4, at 25 degrees C, I = 0.1, catalysed by beta-trypsin is comparable with those for the inhibition of trypsin by benzamidine and 4-aminobenzamidine, which are widely regarded as the most effective low Mr inhibitors of this enzyme. Computer modelling of the beta-trypsin-E64 adsorptive complex, by energy minimization, molecular dynamics simulation and Poisson-Boltzmann electrostatic-potential calculations, was used to define the probable binding mode of E-64; the ligand lies parallel to the active-centre cleft, anchored principally by the dominant electrostatic interaction of the guanidinium cation at one end of the E-64 molecule with the carboxylate anion of Asp-171 (beta-trypsin numbering from Ile-1) in the S1-subsite, and by the interaction of the carboxylate substituent on C-2 of the epoxide ring at the other end of the molecule with Lys-43; the epoxide ring of E-64 is remote from the catalytic site serine hydroxy group. The possibility that E-64 might bind to the cysteine proteinases clostripain (from Clostridium histolyticum) and alpha-gingivain (one of the extracellular enzymes from phyromonas gingivalis) in a manner analogous to that deduced for the beta-trypsin-E-64 complex is discussed.