Thymic stromal lymphopoietin deficiency attenuates experimental autoimmune encephalomyelitis.

In the present study we examined the role of thymic stromal lymphopoietin (TSLP) in experimental autoimmune encephalomyelitis (EAE). Here, we report that TSLP knock-out (KO) mice display a delayed onset of disease and an attenuated form of EAE. This delayed onset was accompanied by a reduced number of encephalitogenic T helper type 1 (Th1) cells in the central nervous system (CNS) of TSLP KO mice. In addition, CD4(+) and CD8(+) T cells from CNS of TSLP KO mice show a reduced activation status in comparison to wild-type mice. It is noteworthy that we could also show that lymph node cells from TSLP KO mice expanded less efficiently and that interleukin (IL)-6-, interferon (IFN)-γ and tumour necrosis factor (TNF)-α levels were reduced. Furthermore, CD3(+) T cells isolated in the preclinical phase from myelin oligodendrocyte glycoprotein peptide 35-55 (MOG(35-55))-immunized TSLP KO mice showed a reduced response after secondary exposure to MOG(35-55), indicating that differentiation of naive T cells into MOG(35-55)-specific effector and memory T cells was impaired in KO mice. The addition of recombinant TSLP enhanced T cell proliferation during MOG(35-55) restimulation, showing that T cells also respond directly to TSLP. In summary, these data demonstrate that expression of, and immune activation by, TSLP contributes significantly to the immunopathology of EAE.

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.

Effects of CMAP and electrostatic cutoffs on the dynamics of an integral membrane protein: the phospholamban study.

In our effort to understand the microscopic structure and dynamics of phospholamban (PLB), a small integral membrane protein, we have performed a series of 5-20 ns molecular dynamics simulations to explore the influence of environment (solution and lipid bi-layer) and force field (CMAP correction and Ewald summation) on the protein behavior. Under all simulation conditions, we have observed the same major features: existence of two well-defined helical domains at the N- and C-termini, and large-amplitude rigid-body motions of these domains. The average inter-helix angle of PLB was sensitive to the environment. In the methanol and water solution trajectories, the two helical domains tended to adopt a closed orientation, with the inter-helix angle below 90 degrees, while in the lipid bi-layer the domains tend to be in an open conformation, with the inter-helix angle above 90 degrees. Within each studied environment, simulations employing different force field models provided qualitatively similar description of PLB structure and dynamics. The only significant discrepancy was the presence of pi-helical hydrogen bonds in trajectories generated with the standard CHARMM22 force field. Simulations with the CMAP correction, with both cutoff and Ewald electrostatics, exhibited predominantly alpha-helical and some 3(10)-helical hydrogen bonding interactions, and no pi-helical hydrogen bonding, in accord with NMR data. Thus, our results indicate that models including CMAP, with both cutoff and Ewald electrostatics, provide the most realistic description of PLB structure and dynamics. Results obtained from these simulations are in a good agreement with the experimental observables. These include helical secondary structure of PLB, the range explored by the inter-helix angle in methanol, as well as the inter-helix distance and C-terminal helix orientation in the DPPC bi-layer. The observed effect of opening up of the PLB inter-helix angle in the lipid environment relative to solution is also qualitatively reproduced in the simulations, as is the more rigid and compact structure of the C-terminal domain in the membrane relative to solution. The populations of conformations with relatively open inter-domain angles, as well as large fluctuations of this coordinate in DPPC bi-layers allow the N-terminal helix to come into contact with the PLB binding site on the calcium ATPase. Additionally, the presence of a twisting motion around the helical axis enables the helix to orient the correct face to the binding site. Another interesting observation is that the phosphorylation sites Ser(16) and Thr(17) are essentially always accessible to solvent, and presumably also to phosphorylation.

Computational characterization of substrate binding and catalysis in S-adenosylhomocysteine hydrolase.

S-Adenosylhomocysteine (AdoHcy) hydrolase catalyzes the reversible hydrolysis of AdoHcy to adenosine (Ado) and homocysteine (Hcy), playing an essential role in modulating the cellular Hcy levels and regulating activities of a host of methyltransferases in eukaryotic cells. This enzyme exists in an open conformation (active site unoccupied) and a closed conformation (active site occupied with substrate or inhibitor) [Turner, M. A., Yang, X., Yin, D., Kuczera, K., Borchardt, R. T., and Howell, P. L. (2000) Cell Biochem. Biophys. 33, 101-125]. To investigate the binding of natural substrates during catalysis, the computational docking program AutoDock (with confirming calculations using CHARMM) was used to predict the binding modes of various substrates or inhibitors with the closed and open forms of AdoHcy hydrolase. The results have revealed that the interaction between a substrate and the open form of the enzyme is nonspecific, whereas the binding of the substrate in the closed form is highly specific with the adenine moiety of a substrate as the main recognition factor. Residues Thr57, Glu59, Glu156, Gln181, Lys186, Asp190, Met351, and His35 are involved in substrate binding, which is consistent with the crystal structure. His55 in the docked model appears to participate in the elimination of water from Ado through the interaction with the 5'-OH group of Ado. In the same reaction, Asp131 removes a proton from the 4' position of the substrate after the oxidation-reduction reaction in the enzyme. To identify the residues that bind the Hcy moiety, AdoHcy was docked to the closed form of AdoHcy hydrolase. The Hcy tail is predicted to interact with His55, Cys79, Asn80, Asp131, Asp134, and Leu344 in a strained conformation, which may lower the reaction barrier and enhance the catalysis rate.

Structure and dynamics of calcium-activated calmodulin in solution.

Two 4-ns molecular dynamics simulations of calcium loaded calmodulin in solution have been performed, using both standard nonbonded cutoffs and Ewald summation to treat electrostatic interactions. Our simulation results are generally consistent with solution experimental studies of calmodulin structure and dynamics, including NMR, cross-linking, fluorescence and x-ray scattering. The most interesting result of the molecular dynamics simulations is the detection of large-scale structural fluctuations of calmodulin in solution. The globular N- and C-terminal domains tend to move approximately like rigid bodies, with fluctuations of interdomain distances within a 7 A range and of interdomain angles by up to 60 deg. Essential dynamics analysis indicates that the three dominant types of motion involve bending of the central helix in two perpendicular planes and a twist in which the domains rotate in opposite directions around the central helix. In the more realistic Ewald trajectory the protein backbone remains mostly within a 2-3 A root-mean-square distance from the crystal structure, the secondary structure within the domains is conserved and middle part of the central helix becomes disordered. The central helix itself exhibits limited fluctuations, with its bend angle exploring the 0-50 degrees range and the end-to-end distance falling in 39-43 A. The results of the two simulations were similar in many respects. However, the cutoff trajectory exhibited a larger deviation from the crystal, loss of several helical hydrogen bonds in the N-terminal domain and lack of structural disorder in the central helix.

Probing the differences between rat liver outer mitochondrial membrane cytochrome b5 and microsomal cytochromes b5.

Two distinct forms of cytochrome b5 exist in the rat hepatocyte. One is associated with the membrane of the endoplasmic reticulum (microsomal, or Mc, cyt b5) while the other is associated with the outer membrane of liver mitochondria (OM cyt b5). Rat OM cyt b5, the only OM cyt b5 identified so far, has a significantly more negative reduction potential and is substantially more stable toward chemical and thermal denaturation than Mc cytochromes b5. In addition, hemin is kinetically trapped in rat OM cyt b5 but not in the Mc proteins. As a result, no transfer of hemin from rat OM cyt b5 to apomyoglobin is observed at pH values as low as 5.2, nor can the thermodyamically favored ratio of hemin orientational isomers be achieved under physiologically relevant conditions. These differences are striking given the similarity of the respective protein folds. A combined theoretical and experimental study has been conducted in order to probe the structural basis behind the remarkably different properties of rat OM and Mc cytochromes b5. Molecular dynamics (MD) simulations starting from the crystal structure of bovine Mc cyt b5 revealed a conformational change that exposes several internal residues to the aqueous environment. The new conformation is equivalent to the "cleft-opened" intermediate observed in a previously reported MD simulation of bovine Mc cyt b5 [Storch, E. M., and Daggett, V. (1995) Biochemistry 34, 9682-9693]. The rat OM protein does not adopt a comparable conformation in MD simulations, thus restricting access of water to the protein interior. Subsequent comparisons of the protein sequences and structures suggested that an extended hydrophobic network encompassing the side chains of Ala-18, Ile-32, Leu-36, and Leu-47 might contribute to the inability of rat OM cyt b5 to adopt the cleft-opened conformation and, hence, stabilize its fold relative to the Mc isoforms. A corresponding network is not present in bovine Mc cyt b5 because positions 18, 32, and 47, are occupied by Ser, Leu, and Arg, respectively. To probe the roles played by Ala-18, Ile-32, and Leu-47 in endowing rat OM cyt b5 with its unusual structural properties, we have replaced them with the corresponding residues in bovine Mc cyt b5. Hence, the I32L (single), A18S/L47R (double), and A18S/L47R/I32L (triple) mutants of rat OM cyt b5 were prepared. The stability of these proteins was found to decrease in the following order: WT rat OM > rat OM I32L > rat OM A18S/L47R > rat OM A18S/L47R/I32L > bovine Mc cyt b5. The decrease in stability of the rat OM protein correlates with the extent to which the hydrophobic cluster involving the side chains of residues 18, 32, 36, and 47 has been disrupted. Complete disruption of the hydrophobic network in the triple mutant is confirmed in a 2.0 A resolution crystal structure of the protein. Disruption of the hydrophobic network also facilitates hemin loss at pH 5.2 for the double and triple mutants, with the less stable triple mutant exhibiting the greater rate of hemin transfer to apomyoglobin. Finally, 1H NMR spectroscopy and side-by-side comparisons of the crystal structures of bovine Mc, rat OM, and rat OM A18S/L47R/I32L cyt b5 allowed us to conclude that the nature of residue 32 plays a key role in controlling the relative stability of hemin orientational isomers A and B in rat OM cyt b5. A similar analysis led to the conclusion that Leu-70 and Ser-71 play a pivotal role in stabilizing isomer A relative to isomer B in Mc cytochromes b5.

Transitions from alpha to pi helix observed in molecular dynamics simulations of synthetic peptides.

Molecular dynamics simulations were carried out for three 13-residue peptides of the form AcNH-A-A-E-X-A-E-A-H-A-A-E-K-A-CONH(2) with X = A, F, and W. All three peptides exhibited unexpected dynamical behavior, undergoing a transition from an alpha-helical to a pi-helical structure in the course of 5-ns trajectories in aqueous solution. Analysis of peptide length, accessible surface, interaction energies, hydrogen bonding, and dihedral angles was consistent with alpha --> pi transitions at 2800, 500, and 800 ps for X = A, F and W, respectively. The transitions occurred sequentially and cooperatively, propagating from the C- to the N-terminus for X = A and W and from the center toward both termini for X = F. The time scale of the overall transition ranged from 300 to 500 ps. For all three peptides the backbone structural transition was accompanied by a concerted rearrangement of the charged side chains, including a 3 A increase in the distance between carboxylate groups of Glu-3 and Glu-6. During the transition the peptide backbone hydrogen-bonding patterns were disrupted at the interface between the alpha-helical and nascent pi-helical regions, with peptide groups forming water-bridged hydrogen bonds. The peptide structures exhibited significant fluidity, with individual residues sampling alpha-, pi-, and 3(10)-helical conformations, as well as a "coil" state, without any intramolecular hydrogen bonds. The studied peptides have been designed to form alpha-helices when incorporated in novel hemoprotein model compounds, peptide-sandwiched mesohemes, which consist of two identical peptides covalently attached to an Fe(III) mesoporphyrin [Liu, D., Williamson, D. A., Kennedy, M. L., Williams, T. D., Morton, M. M., and Benson, D. R. (1999) J. Am. Chem. Soc. 121, 11798-11812]. The possibility of adopting pi-helical structures by the constituent peptides may influence the properties of the hemoprotein models.

O2 activation by nonheme iron complexes: A monomeric Fe(III)-Oxo complex derived from O2.

Iron species with terminal oxo ligands are implicated as key intermediates in several synthetic and biochemical catalytic cycles. However, there is a dearth of structural information regarding these types of complexes because their instability has precluded isolation under ambient conditions. The isolation and structural characterization of an iron(III) complex with a terminal oxo ligand, derived directly from dioxygen (O2), is reported. A stable structure resulted from placing the oxoiron unit within a synthetic cavity lined with hydrogen-bonding groups. The cavity creates a microenvironment around the iron center that aids in regulating O2 activation and stabilizing the oxoiron unit. These cavities share properties with the active sites of metalloproteins, where function is correlated strongly with site structure.

Substrate binding stabilizes S-adenosylhomocysteine hydrolase in a closed conformation.

Comparison of crystal structures of S-adenosylhomocysteine (AdoHcy) hydrolase in the substrate-free, NAD(+) form [Hu, Y., Komoto, J., Huang, Y., Gomi, T., Ogawa, H., Takata, Y., Fujioka, M., and Takusagawa, F. (1999) Biochemistry 38, 8323-8333] and a substrate-bound, NADH form [Turner, M. A., Yuan, C.-S., Borchardt, R. T., Hershfield, M. S., Smith, G. D., and Howell, P. L. (1998) Nat. Struct. Biol. 5, 369-376] indicates large differences in the spatial arrangement of the catalytic and NAD(+) binding domains. The substrate-free, NAD(+) form exists in an "open" form with respect to catalytic and NAD(+) binding domains, whereas the substrate-bound, NADH form exists in a closed form with respect to those domains. To address whether domain closure is induced by substrate binding or its subsequent oxidation, we have measured the rotational dynamics of spectroscopic probes covalently bound to Cys(113) and Cys(421) within the catalytic and carboxyl-terminal domains. An independent domain motion is associated with the catalytic domain prior to substrate binding, suggesting the presence of a flexible hinge element between the catalytic and NAD(+) binding domains. Following binding of substrates (i.e., adenosine or neplanocin A) or a nonsubstrate (i.e., 3'-deoxyadenosine), the independent domain motion associated with the catalytic domain is essentially abolished. Likewise, there is a substantial decrease in the average hydrodynamic volume of the protein that is consistent with a reduction in the overall dimensions of the homotetrameric enzyme following substrate binding and oxidation observed in earlier crystallographic studies. Thus, the catalytic and NAD(+) binding domains are stabilized to form a closed active site through interactions with the substrate prior to substrate oxidation.

Metal-catalyzed oxidation of brain-derived neurotrophic factor (BDNF): selectivity and conformational consequences of histidine modification.

We have studied the metal-catalyzed oxidation (MCO) of brain-derived neurotrophic factor (BDNF) with regard to target sites and potential conformational changes of the protein. The exposure of BDNF to three different levels of ascorbate/Cu(II)/O2 [20 microM Cu(II), 2 mM ascorbate (level 1); 20 microM Cu(II), 4 mM ascorbate (level 2); 40 microM Cu(II), 4 mM ascorbate (level 3)], chosen based on the extent of chemical modification of Met and His, respectively, resulted in the exclusive oxidation of a buried Met residue, Met92, at level 1 but in the predominant oxidation of His at level 3. His modification had a significant impact on the structure of BDNF, as quantified by CD and ANSA fluorescence measurements, while Met oxidation had not, also assessed through complementary oxidation of BDNF through hydrogen peroxide. Our ultimate objective was the correlation of the surface exposure of an oxidized His residue in a protein with potential effects on the conformational integrity of the oxidized protein. In a series of three proteins, human growth hormone (hGH), human relaxin (hR1x), and BDNF, we have now observed that His oxidation is paralleled by significant conformational changes when the target His residue is more surface exposed (hR1x, BDNF) while conformational consequences of His modification are less significant when the target His residues are more buried in the interior of the protein (hGH).

The sensitivity of carboxyl-terminal methionines in calmodulin isoforms to oxidation by H(2)O(2) modulates the ability to activate the plasma membrane Ca-ATPase.

The oxidative modification of methionines within the primary sequence of calmodulin (CaM) results in an inability to activate the PM-Ca-ATPase fully, and may contribute to alterations in calcium homeostasis under conditions of oxidative stress. To identify differences in the sensitivities of CaM isoforms to oxidative modification, we have compared the function and patterns of oxidative modification resulting from the exposure of CaM isolated from bovine testes and wheat germ to H(2)O(2). In comparison to CaM isolated from wheat germ, vertebrate CaM is functionally resistant to oxidant-induced loss of function. The decreased functional sensitivity of vertebrate CaM correlates with a 75 +/- 3% reduction in the rate of oxidative modification of a methionine near the carboxyl terminus (i.e., Met(144) or Met(145)). The extent of oxidative modification to other methionines in these CaM isoforms is similar. These results suggest that the sensitivity of Met(144) or Met(145) to oxidation modulates the ability of CaM to activate the PM-Ca-ATPase. Consistent with this interpretation, a CaM mutant in which glutamines were substituted for Met(144) and Met(145) fully activates the PM-Ca-ATPase irrespective of the oxidative modification of the other seven methionines to their corresponding methionine sulfoxides. The extent of oxidative modification to individual methionines in vertebrate CaM by H(2)O(2) correlates with the time-averaged surface accessibility of individual sulfurs calculated from molecular dynamics simulations. Thus, the sensitivity of individual methionines to oxidative modification is directly related to the solvent accessibility. These results indicate that sequence differences between vertebrate and plant CaM alter the sensitivity of methionines near the carboxyl terminus to oxidative modification because of alterations in their solvent accessibility. We suggest that these sequence differences between CaM isoforms have a regulatory role in modulating the functional sensitivity of CaM to conditions of oxidative stress.

Ligand binding inside the cavities of lacunar and saddle-shaped cyclidene complexes: molecular mechanics and molecular dynamics studies.

Cobalt(II) complexes with tetradentate macrocyclic cyclidene ligands are known to coordinate one additional axial base molecule, leaving the sixth vacant coordination site at the metal ion available for small ligand (e.g., O2) binding. Molecular mechanics and molecular dynamics simulations provide a microscopic view of 1-methylimidazole (MeIm) binding within the cavities of several lacunar (bridged) and saddle-shaped (unbridged) cyclidenes and uncover the roles of the bridges and the walls of the clefts in steric protection of the cobalt(II) coordination site. Short bridges (C3 and C6) prevent inside-the-cavity MeIm binding because of severe ligand distortions leading to high-energy penalties (58 and 25 kcal/mol, respectively), while long bridges (C8 and C12) flip away from the MeIm binding site, allowing for penalty-free MeIm inclusion. In the unbridged saddle-shaped complex, there is no energy difference between inside- and outside-the-cavity MeIm binding. The preferential existence of the coordinatively unsaturated, five-coordinate species Co(unbrCyc)(MeIm)2+ should therefore be explained by electronic, rather than steric, factors. Molecular dynamics and free energy simulations reveal the presence of a weak (ca. 4 kcal/mol in the gas phase and ca. 2 kcal/mol in methanol solution) noncovalent MeIm binding site at the entrance of the cleft of cobalt(II) unbridged cyclidene, at a distance of about 4 A from the metal ion. The macrocycle geometry remains undistorted at such large Co-N(MeIm) separations, while the cavity opens up by 0.9 A upon covalent MeIm binding (Co-N(MeIm) distance of 2 A). An increase in macrocycle strain energy upon MeIm inclusion is compensated by favorable nonbonded interactions between the incoming base and the walls of the unbridged cyclidene.

Structure and function of S-adenosylhomocysteine hydrolase.

In mammals, S-adenosylhomocysteine hydrolase (AdoHcyase) is the only known enzyme to catalyze the breakdown of S-adenosylhomocysteine (AdoHcy) to homocysteine and adenosine. AdoHcy is the product of all adenosylmethionine (AdoMet)-dependent biological transmethylations. These reactions have a wide range of products, and are common in all facets of biometabolism. As a product inhibitor, elevated levels of AdoHcy suppress AdoMet-dependent transmethylations. Thus, AdoHcyase is a regulator of biological transmethylation in general. The three-dimensional structure of AdoHcyase complexed with reduced nicotinamide adenine dinucleotide phosphate (NADH) and the inhibitor (1'R, 2'S, 3'R)-9-(2',3'-dihyroxycyclopenten-1-yl)adenine (DHCeA) was solved by a combination of the crystallographic direct methods program, SnB, to determine the selenium atom substructure and by treating the multiwavelength anomalous diffraction data as a special case of multiple isomorphous replacement. The enzyme architecture resembles that observed for NAD-dependent dehydrogenases, with the catalytic domain and the cofactor-binding domain each containing a modified Rossmann fold. The two domains form a deep active site cleft containing the cofactor and bound inhibitor molecule. A comparison of the inhibitor complex of the human enzyme and the structure of the rat enzyme, solved without inhibitor, suggests that a 17 degrees rigid body movement of the catalytic domain occurs upon inhibitor/substrate binding.

Molecular dynamics study of disulfide bond influence on properties of an RGD peptide.

Three 1 ns length molecular dynamics simulations of an RGD peptide (Ac-Pen-Arg-Gly-Asp-Cys-NH2, with Pen denoting penicillamine) have been performed in aqueous solution, one for the disulfide bridged, and two for the unbridged form. The trajectories were analyzed to identify conformations explored by the two forms and to calculate several properties: NMR vicinal coupling constants, order parameters, dipole moments and diffusion coefficients, in an effort to describe the physical role of the disulfide bond. The cyclic peptide was able to explore several distinct backbone conformations centered around a turn-extended-turn structure. However, its flexibility was limited and it appeared to be 'locked in' into a a family of structures characterized by a high dipole moment and a well-defined conformation of the pharmacophore, which has been previously identified as biologically active. Excellent agreement between the simulated and observed NMR vicinal coupling constants indicates that realistic structures were sampled in the cyclic peptide simulation. The linear form of the peptide was much more flexible than the cyclic one. In the two independent 1 ns simulations of the linear form the explored conformations could be roughly grouped into two classes, of cyclic-like and extended type. Within each simulation the peptide switched between the two classes of structures several times. Exact matches between conformations in the two linear peptide simulations were not found; several conformational regions with backbone rms deviations below 1A were identified, suggesting that representative structures of the linear form have also been identified. In the linear peptide simulations the RGD pharmacophore is able to adopt a wide range of conformations, including the one preferred by the cyclic form. The lower biological activity of the linear peptide compared to the cyclic one may be correlated with the lower population of this structure in the absence of the disulfide bond.

Free-energy simulations of the retinal cis --> trans isomerization in bacteriorhodopsin.

Free-energy profiles for ground-state cis --> trans isomerization of retinal in vacuum, in solution, and in the protein bacteriorhodopsin are calculated using free-energy simulations. The free-energy barriers in the protein were 9 kcal/mol for ionized Asp85 and 14 kcal/mol for neutral Asp85, significantly lower than those found in solution (18 kcal/mol) or vacuum (19 kcal/mol). Therefore, bacteriorhodopsin can be said to act as a catalyst in the isomerization. The barrier in the protein is due mainly to stabilization of the transition state through favorable nonbonded interactions with the protein part of the system, with internal strain and interactions with solvent playing minor roles. The protonated Asp85 simulation models the behavior of the system in the N --> O transition. Our calculated 14 kcal/mol barrier and 4-ms relaxation time for this process are in excellent agreement with experimentally measured values of 12 kcal/mol and 5 ms, respectively. The ionized Asp85 simulation models two hypothetical processes: the N --> O transition with a proton removed from Asp85 and the initial BR568 --> L transition on the ground-state energy surface. The cis-trans isomerization barrier in this system is 9 kcal/mol, the lowest of all the studied cases. The presence of the charged carboxylate group in the ionized Asp85 system leads to strong stabilization of the transition state by interactions with the surroundings and changes the distance between Asp85 and the Schiff base proton compared to the corresponding distance in the neutral Asp85 system. This suggests that the protonation of Asp85 plays an important role in regulating access to the Schiff base proton. For both Asp85 ionization states the calculated cis-trans free-energy difference was close to 0, indicating that the protein can accommodate both retinal isomers equally well. The computed negligible difference between the N and O free-energy levels is in accord with experimental data.

All-atom empirical potential for molecular modeling and dynamics studies of proteins.

New protein parameters are reported for the all-atom empirical energy function in the CHARMM program. The parameter evaluation was based on a self-consistent approach designed to achieve a balance between the internal (bonding) and interaction (nonbonding) terms of the force field and among the solvent-solvent, solvent-solute, and solute-solute interactions. Optimization of the internal parameters used experimental gas-phase geometries, vibrational spectra, and torsional energy surfaces supplemented with ab initio results. The peptide backbone bonding parameters were optimized with respect to data for N-methylacetamide and the alanine dipeptide. The interaction parameters, particularly the atomic charges, were determined by fitting ab initio interaction energies and geometries of complexes between water and model compounds that represented the backbone and the various side chains. In addition, dipole moments, experimental heats and free energies of vaporization, solvation and sublimation, molecular volumes, and crystal pressures and structures were used in the optimization. The resulting protein parameters were tested by applying them to noncyclic tripeptide crystals, cyclic peptide crystals, and the proteins crambin, bovine pancreatic trypsin inhibitor, and carbonmonoxy myoglobin in vacuo and in crystals. A detailed analysis of the relationship between the alanine dipeptide potential energy surface and calculated protein φ, χ angles was made and used in optimizing the peptide group torsional parameters. The results demonstrate that use of ab initio structural and energetic data by themselves are not sufficient to obtain an adequate backbone representation for peptides and proteins in solution and in crystals. Extensive comparisons between molecular dynamics simulations and experimental data for polypeptides and proteins were performed for both structural and dynamic properties. Energy minimization and dynamics simulations for crystals demonstrate that the latter are needed to obtain meaningful comparisons with experimental crystal structures. The presented parameters, in combination with the previously published CHARMM all-atom parameters for nucleic acids and lipids, provide a consistent set for condensed-phase simulations of a wide variety of molecules of biological interest.

Free energy simulations of axial contacts in sickle-cell hemoglobin.

Molecular dynamics simulations have been used to investigate the thermodynamic stability of axial contacts in sickle-cell hemoglobin (HbS). Free energy changes were evaluated for the point mutation beta 121 Glu --> Gln in the axial contact region of HbS crystals. The calculations predict a free energy change of-3.6 kcal/mol per contact for the mutation, which is in qualitative agreement with experimental observations of aggravated sickling found in the double mutant Hb D Los Angeles (beta 6 Glu --> Val. beta 121 Glu --> Gln) relative to HbS (beta 6 Glu --> Val). The beta 121 Glu is sequestered in a salt link with beta 17 Lys located on the same polypeptide chain, making the Glu interactions with its surroundings similar in aggregates and individual hemoglobins. Due to this cancellation of the large electrostatic Glu contributions, the weak nonspecific interactions between the Gln and the neighboring polypeptide chain are the main contributing factor to the enhanced aggregation of Hb D Los Angeles relative to HbS. Together with the previous study of the lateral contact [K. Kuczera et al. (1990) Proceedings of the National Academy of Science USA, Vol. 87, pp, 8481-8485], the present results provide a more complete picture of the forces driving the sickling aggregation. A comparison of different treatments of internal flexibility in free energy simulations and analysis of rate of convergence of the different calculated properties has also been performed.

Nonexponential relaxation after ligand dissociation from myoglobin: a molecular dynamics simulation.

Molecular dynamics simulations of myoglobin after ligand photodissociation show that the out-of-plane motion of the heme iron has a rapid subpicosecond phase followed by a slower nonexponential process involving more global protein relaxation. Individual trajectories show rather different behavior, suggesting there is an inhomogeneous component to the relaxation. The calculated time dependence of the iron motion over 100 ps is in excellent agreement with the frequency shift of band III of the heme group [see Lim, M., Jackson, T. A. & Anfinrud, P. A. (1993) Proc. Natl. Acad. Sci. USA 90, 5801-5804]. If that the barrier to rebinding depends on the out-of-plane iron position, the time dependence obtained from the simulation can explain the nonexponential room-temperature geminate recombination of NO.

Ligand binding and protein relaxation in heme proteins: a room temperature analysis of NO geminate recombination.

Ultrafast absorption spectroscopy is used to study heme-NO recombination at room temperature in aqueous buffer on time scales where the ligand cannot leave its cage environment. While a single barrier is observed for the cage recombination of NO with heme in the absence of globin, recombination in hemoglobin and myoglobin is nonexponential. Examination of hemoglobin with and without inositol hexaphosphate points to proximal constraints as important determinants of the geminate rebinding kinetics. Molecular dynamics simulations of myoglobin and heme-imidazole subsequent to ligand dissociation were used to investigate the transient behavior of the Fe-proximal histidine coordinate and its possible involvement in geminate recombination. The calculations, in the context of the absorption measurements, are used to formulate a distinction between nonexponential rebinding that results from multiple protein conformations (substates) present at equilibrium or from nonequilibrium relaxation of the protein triggered by a perturbation such as ligand dissociation. The importance of these two processes is expected to depend on the time scale of rebinding relative to equilibrium fluctuations and nonequilibrium relaxation. Since NO rebinding occurs on the picosecond time scale of the calculated myoglobin relaxation, a time-dependent barrier is likely to be an important factor in the observed nonexponential kinetics. The general implications of the present results for ligand binding in heme proteins and its time and temperature dependence are discussed. It appears likely that, at low temperatures, inhomogeneous protein populations play an important role and that as the temperature is raised, relaxation effects become significant as well.

Free energy of sickling: A simulation analysis.

Molecular dynamics simulations were performed to calculate the difference between the dimerization free energies of normal human deoxyhemoglobin (HbA) and the mutant sickle-cell deoxyhemoglobin HbS (Glu-beta 6----Val) for one of the lateral contacts in the HbS x-ray structure. The simulations yield a value of--15 kcal/mol. Although there is no quantitative experimental value for comparison, this is in qualitative agreement with the experimental result that HbS self-assembles into multistranded fibers that are responsible for erythrocyte sickling, while HbA does not. The free-energy difference was decomposed into enthalpic and entropic terms, both of which are significant, and the contributions of individual protein residues and of the solvent were examined. Electrostatic effects play the dominant role in favoring dimerization of HbS compared with HbA; van der Waals interactions make a negligible contribution to the difference. Both differential solvation and protein-protein interactions are important. Interactions within the donor tetramer (i.e., that containing the Glu-beta 6 mutation site), as well as those with the acceptor tetramer, contribute to the preferential free energy of dimerization of HbS.