Search

A fundamental appreciation for how biological macromolecules work requires knowledge of structure and dynamics. Molecular dynamics simulations provide powerful tools for the exploration of the conformational energy landscape accessible to these molecules, and the rapid increase in computational power coupled with improvements in methodology makes this an exciting time for the application of simulation to structural biology. In this Perspective we survey two areas, protein folding and enzymatic catalysis, in which simulations have contributed to a general understanding of mechanism. We also describe results for the [F.sub.1] ATPase molecular motor and the Src family of signaling proteins as examples of applications of simulations to specific biological systems.

[F.sub.1][F.sub.o]-ATP synthase is the enzyme responsible for most of the ATP synthesis in living systems. The catalytic domain [F.sub.1] of the [F.sub.1][F.sub.o] complex, [F.sub.1]-ATPase, has the ability to hydrolyze ATP. A fundamental problem in the development of a detailed mechanism for this enzyme is that it has not been possible to determine experimentally the relation between the ligand binding affinities measured in solution and the different conformations of the catalytic [beta] subunits ([[beta].sub.TP], [[beta].sub.DP],[[beta].sub.E]) observed in the crystal structures of the mitochondrial enzyme, M[F.sub.1]. Using free energy difference simulations for the hydrolysis reaction ATP+[H.sub.2]O [right arrow] ADP+[P.sub.i] in the [[beta].sub.TP] and [[beta].sub.DP] sites and unisite hydrolysis data, we are able to identify [[beta].sub.TP] as the 'tight' ([K.sub.D] = [10.sup.-12] M, M[F.sub.1]) binding site for ATP and [[beta].sub.DP] as the 'loose' site. An energy decomposition analysis demonstrates how certain residues, some of which have been shown to be important in catalysis, modulate the free energy of the hydrolysis reaction in the [[beta].sub.TP] and [[beta].sub.DP] sites, even though their structures are very similar. Combined with the recently published simulations of the rotation cycle of [F.sub.1]-ATPase, the present results make possible a consistent description of the binding change mechanism of [F.sub.1]-ATPase at an atomic level of detail.

Molecular dynamics simulations of the essential part of the photosynthetic reaction center of Rhodopseudomonas viridis was done to examine the coupling of electron transfer with protein structure and dynamics of the protein complex. The experiment showed the presence of a coupling reaction between the thermal motion of segments of photosynthetic reaction centers with electron transfer and that this coupling process regulated the rate of electron transfer. A fast and temperature-independent dielectron relaxation of the protein matrix stabilized the dissociated electron pair.

Molecular dynamics simulations of myoglobin after ligand photodissociation reveal that the out-of-plan movement of the heme iron has a quick subsecond stage, after which there is a more gradual nonexponential process having greater global protein relaxation. The different trajectory behaviors exhibited indicate an inhomogeneous element to the relaxation.

F[sub 1]F[sub 0]-ATP synthase is the enzyme responsible for most of the ATP synthesis in living systems. The catalytic domain F[sub 1] of the F[sub 1]F[sub 0] complex, F[sub 1]-ATPase, has the ability to hydrolyze ATP. A fundamental problem in the development of a detailed mechanism for this enzyme is that it has not been possible to determine experimentally the relation between the ligand binding affinities measured in solution and the different conformations of the catalytic β subunits (β[sub TP], β[sub DP], β[sub E]) observed in the crystal structures of the mitochondrial enzyme, MF[sub 1]. Using free energy difference simulations for the hydrolysis reaction ATP+H[sub 2]O → ADP+ Pi in the β[sub TP] and β[sub DP] sites and unisite hydrolysis data, we are able to identify β[sub TP] as the "tight" (K[sub D] = 10[sup -12] M, MF[sub 1]) binding site for ATP and β[sub DP] as the "loose" site. An energy decomposition analysis demonstrates how certain residues, some of which have been shown to be important in catalysis, modulate the free energy of the hydrolysis reaction in the β[sub TP] and β[sub DP] sites, even though their structures are very similar. Combined with the recently published simulations of the rotation cycle of F[sub 1]-ATPase, the present results make possible a consistent description of the binding change mechanism of F[sub 1]-ATPase at an atomic level of detail. [ABSTRACT FROM AUTHOR]

A vaccine which is effective against the HIV virus is considered to be the best solution to the ongoing global HIV/AIDS epidemic. In the past thirty years, numerous attempts to develop an effective vaccine have been made with little or no success, due, in large part, to the high mutability of the virus. More recent studies showed that a vaccine able to elicit broadly neutralizing antibodies (bnAbs), that is, antibodies that can neutralize a high fraction of global virus variants, has promise to protect against HIV. Such a vaccine has been proposed to involve at least three separate stages: First, activate the appropriate precursor B cells; second, shepherd affinity maturation along pathways toward bnAbs; and, third, polish the Ab response to bind with high affinity to diverse HIV envelopes (Env). This final stage may require immunization with a mixture of Envs. In this paper, we set up a framework based on theory and modeling to design optimal panels of antigens to use in such a mixture. The designed antigens are characterized experimentally and are shown to be stable and to be recognized by known HIV antibodies., Competing Interests: The authors declare no competing interest.

Two mechanisms have been proposed for the function of motor proteins: The power stroke and the Brownian ratchet. The former refers to generation of a large downhill free energy gradient over which the motor protein moves nearly irreversibly in making a step, whereas the latter refers to biasing or rectifying the diffusive motion of the motor. Both mechanisms require input of free energy, which generally involves the processing of an ATP (adenosine 5'-triphosphate) molecule. Recent advances in experiments that reveal the details of the stepping motion of motor proteins, together with computer simulations of atomistic structures, have provided greater insights into the mechanisms. Here, we compare the various models of the power stroke and the Brownian ratchet that have been proposed. The 2 mechanisms are not mutually exclusive, and various motor proteins employ them to different extents to perform their biological function. As examples, we discuss linear motor proteins Kinesin-1 and myosin-V, and the rotary motor F 1 -ATPase, all of which involve a power stroke as the essential element of their stepping mechanism., Competing Interests: The authors declare no conflict of interest.

Our understanding of the rotary-coupling mechanism of F 1 -ATPase has been greatly enhanced in the last decade by advances in X-ray crystallography, single-molecular imaging, and theoretical models. Recently, Volkán-Kacsó and Marcus [S. Volkán-Kacsó, R. A. Marcus, Proc. Natl. Acad. Sci. U.S.A. 112, 14230 (2015)] presented an insightful thermodynamic model based on the Marcus reaction theory coupled with an elastic structural deformation term to explain the observed γ-rotation angle dependence of the adenosine triphosphate (ATP)/adenosine diphosphate (ADP) exchange rates of F 1 -ATPase. Although the model is successful in correlating single-molecule data, it is not in agreement with the available theoretical results. We describe a revision of the model, which leads to consistency with the simulation results and other experimental data on the F 1 -ATPase rotor compliance. Although the free energy liberated on ATP hydrolysis by F 1 -ATPase is rapidly dissipated as heat and so cannot contribute directly to the rotation, we show how, nevertheless, F 1 -ATPase functions near the maximum possible efficiency. This surprising result is a consequence of the differential binding of ATP and its hydrolysis products ADP and P i along a well-defined pathway., Competing Interests: The authors declare no conflict of interest.

Small multidrug resistance (SMR) pumps represent a minimal paradigm of proton-coupled membrane transport in bacteria, yet no high-resolution structure of an SMR protein is available. Here, atomic-resolution structures of the Escherichia coli efflux-multidrug resistance E ( EmrE ) multidrug transporter in ligand-bound form are refined using microsecond molecular dynamics simulations biased using low-resolution data from X-ray crystallography. The structures are compatible with existing mutagenesis data as well as NMR and biochemical experiments, including pKas of the catalytic glutamate residues and the dissociation constant ([Formula: see text]) of the tetraphenylphosphonium + cation. The refined structures show the arrangement of residue side chains in the EmrE active site occupied by two different ligands and in the absence of a ligand, illustrating how EmrE can adopt structurally diverse active site configurations. The structures also show a stable, well-packed binding interface between the helices H4 of the two monomers, which is believed to be crucial for EmrE dimerization. Guided by the atomic details of this interface, we design proteolysis-resistant stapled peptides that bind to helix H4 of an EmrE monomer. The peptides are expected to interfere with the dimerization and thereby inhibit drug transport. Optimal positions of the peptide staple were determined using free-energy simulations of peptide binding to monomeric EmrE Three of the four top-scoring peptides selected for experimental testing resulted in significant inhibition of proton-driven ethidium efflux in live cells without nonspecific toxicity. The approach described here is expected to be of general use for the design of peptide therapeutics.

The rotary motor enzyme FoF1-ATP synthase uses the proton-motive force across a membrane to synthesize ATP from ADP and Pi (H2PO4(-)) under cellular conditions that favor the hydrolysis reaction by a factor of 2 × 10(5). This remarkable ability to drive a reaction away from equilibrium by harnessing an external force differentiates it from an ordinary enzyme, which increases the rate of reaction without shifting the equilibrium. Hydrolysis takes place in the neighborhood of one conformation of the catalytic moiety F1-ATPase, whose structure is known from crystallography. By use of molecular dynamics simulations we trap a second structure, which is rotated by 40° from the catalytic dwell conformation and represents the state associated with ATP binding, in accord with single-molecule experiments. Using the two structures, we show why Pi is not released immediately after ATP hydrolysis, but only after a subsequent 120° rotation, in agreement with experiment. A concerted conformational change of the α3β3 crown is shown to induce the 40° rotation of the γ-subunit only when the βE subunit is empty, whereas with Pi bound, βE serves as a latch to prevent the rotation of γ. The present results provide a rationalization of how F1-ATPase achieves the coupling between the small changes in the active site of βDP and the 40° rotation of γ.

Pentameric ligand-gated ion channels (pLGICs) play a central role in intercellular communication in the nervous system and are involved in fundamental processes such as attention, learning, and memory. They are oligomeric protein assemblies that convert a chemical signal into an ion flux through the postsynaptic membrane, but the molecular mechanism of gating ions has remained elusive. Here, we present atomistic molecular dynamics simulations of the prokaryotic channels from Gloeobacter violaceus (GLIC) and Erwinia chrysanthemi (ELIC), whose crystal structures are thought to represent the active and the resting states of pLGICs, respectively, and of the eukaryotic glutamate-gated chloride channel from Caenorhabditis elegans (GluCl), whose open-channel structure was determined complexed with the positive allosteric modulator ivermectin. Structural observables extracted from the trajectories of GLIC and ELIC are used as progress variables to analyze the time evolution of GluCl, which was simulated in the absence of ivermectin starting from the structure with bound ivermectin. The trajectory of GluCl with ivermectin removed shows a sequence of structural events that couple agonist unbinding from the extracellular domain to ion-pore closing in the transmembrane domain. Based on these results, we propose a structural mechanism for the allosteric communication leading to deactivation/activation of the GluCl channel. This model of gating emphasizes the coupling between the quaternary twisting and the opening/closing of the ion pore and is likely to apply to other members of the pLGIC family.

Base excision repair of genotoxic nucleobase lesions in the genome is critically dependent upon the ability of DNA glycosylases to locate rare sites of damage embedded in a vast excess of undamaged DNA, using only thermal energy to fuel the search process. Considerable interest surrounds the question of how DNA glycosylases translocate efficiently along DNA while maintaining their vigilance for target damaged sites. Here, we report the observation of strandwise translocation of 8-oxoguanine DNA glycosylase, MutM, along undamaged DNA. In these complexes, the protein is observed to translocate by one nucleotide on one strand while remaining untranslocated on the complementary strand. We further report that alterations of single base-pairs or a single amino acid substitution (R112A) can induce strandwise translocation. Molecular dynamics simulations confirm that MutM can translocate along DNA in a strandwise fashion. These observations reveal a previously unobserved mode of movement for a DNA-binding protein along the surface of DNA.

Large conformational transitions play an essential role in the function of many proteins, but experiments do not provide the atomic details of the path followed in going from one end structure to the other. For the hemoglobin tetramer, the transition path between the unliganded (T) and tetraoxygenated (R) structures is not known, which limits our understanding of the cooperative mechanism in this classic allosteric system, where both tertiary and quaternary changes are involved. The conjugate peak refinement algorithm is used to compute an unbiased minimum energy path at atomic detail between the two end states. Although the results confirm some of the proposals of Perutz [Perutz MF (1970) Stereochemistry of cooperative effects in haemoglobin. Nature 228:726-734], the subunit motions do not follow the textbook description of a simple rotation of one αβ-dimer relative to the other. Instead, the path consists of two sequential quaternary rotations, each involving different subdomains and axes. The quaternary transitions are preceded and followed by phases of tertiary structural changes. The results explain the recent photodissociation measurements, which suggest that the quaternary transition has a fast (2 μs) as well as a slow (20 μs) component and provide a testable model for single molecule FRET experiments.

Molecular dynamics (MD) simulations provide essential information about the thermodynamics and dynamics of proteins. To construct the free-energy surface from equilibrium trajectories, it is necessary to group the individual snapshots in a meaningful way. The inherent structures (IS) are shown to provide an appropriate discretization of the trajectory and to avoid problems that can arise in clustering algorithms that have been employed previously. The IS-based approach is illustrated with a 30-ns room temperature "native" state MD simulation of a 10-residue peptide in a beta-hairpin conformation. The transitions between the IS are used to construct a configuration space network from which a one-dimensional free-energy profile is obtained with the mincut method. The results demonstrate that the IS approach is useful and that even for this simple system, there exists a nontrivial organization of the native state into several valleys separated by barriers as high as 3 kcal/mol. Further, by introducing a coarse-grained network, it is demonstrated that there are multiple pathways connecting the valleys. This scenario is hidden when the snapshots of the trajectory are used directly with rmsd clustering to compute the free-energy profile. Application of the IS approach to the native state of the PDZ2 signaling domain indicates its utility for the study of biologically relevant systems.

PR65 is the two-layered (alpha-alpha solenoid) HEAT-repeat (Huntingtin, elongation factor 3, a subunit of protein phosphatase 2A, PI3 kinase target of rapamycin 1) scaffold of protein phosphatase PP2A. Molecular dynamics simulations predict that, at forces expected in living systems, PR65 undergoes (visco-)elastic deformations in response to pulling/pushing on its ends. At lower forces, smooth global flexural and torsional changes occur via even redistribution of stress along the hydrophobic core of the molecule. At intermediate forces, helix-helix separation along one layer ("fracturing") leads to global relaxation plus loss of contact in the other layer to unstack the affected units. Fracture sites are determined by unusual sequences in contiguous interhelix turns. Normal mode analysis of the heterotrimeric PP2A enzyme reveals that its ambient conformational fluctuations are dominated by elastic deformations of PR65, which introduce a mechanical linkage between the separately bound regulatory and catalytic subunits. PR65-dominated fluctuations of PP2A have the effect of opening and closing the enzyme's substrate binding/catalysis interface, as well as altering the positions of certain catalytic residues. These results suggest that substrate binding/catalysis are sensitive to mechanical force. Force could be imposed from the outside (e.g., in PP2A's response to spindle tension) or arise spontaneously (e.g., in PP2A's interaction with unstructured proteins such as Tau, a microtubule-associated Alzheimer's-implicated protein). The presented example supports the view that conformation and function of protein complexes can be modulated by mechanical energy inputs, as well as by chemical energy inputs from ligand binding. Given that helical-repeat proteins are involved in many cellular processes, the findings also encourage the view that mechanical forces may be of widespread importance.

Each step of the kinesin motor involves a force-generating molecular rearrangement. Although significant progress has been made in elucidating the broad features of the kinesin mechanochemical cycle, molecular details of the force generation mechanism remain a mystery. Recent molecular dynamics simulations have suggested a mechanism in which the forward drive is produced when the N-terminal cover strand forms a beta-sheet with the neck linker to yield the cover-neck bundle. We tested this proposal by comparing optical trapping motility measurements of cover strand mutants with the wild-type. Motility data, as well as kinetic analyses, revealed impairment of the force-generating capacity accompanied by a greater load dependence in the mechanochemical cycle. In particular, a mutant with the cover strand deleted functioned only marginally, despite the fact that the cover strand, the N-terminal "dangling end," unlike the neck linker and nucleotide-binding pocket, is not involved with any previously considered energy transduction pathway. Furthermore, a constant assisting load, likely in lieu of a power stroke, was shown to rescue forward motility in the cover strand deletion mutant. Our results support a stepping mechanism driven by dynamic cover-neck bundle formation. They also suggest a strategy to generate motors with altered mechanical characteristics by targeting the force-generating element.

A fundamental problem in the analysis of protein folding and other complex reactions in which the entropy plays an important role is the determination of the activation free energy from experimental measurements or computer simulations. This article shows how to combine minimum-cut-based free energy profiles (F(C)), obtained from equilibrium molecular dynamics simulations, with conventional histogram-based free energy profiles (F(H)) to extract the coordinate-dependent diffusion coefficient on the F(C) (i.e., the method determines free energies and a diffusive preexponential factor along an appropriate reaction coordinate). The F(C), in contrast to the F(H), is shown to be invariant with respect to arbitrary transformations of the reaction coordinate, which makes possible partition of configuration space into basins in an invariant way. A "natural coordinate," for which F(H) and F(C) differ by a multiplicative constant (constant diffusion coefficient), is introduced. The approach is illustrated by a model one-dimensional system, the alanine dipeptide, and the folding reaction of a double beta-hairpin miniprotein. It is shown how the results can be used to test whether the putative reaction coordinate is a good reaction coordinate.

F(o)F(1)-ATP synthase manufactures the energy "currency," ATP, of living cells. The soluble F(1) portion, called F(1)-ATPase, can act as a rotary motor, with ATP binding, hydrolysis, and product release, inducing a torque on the gamma-subunit. A coarse-grained plastic network model is used to show at a residue level of detail how the conformational changes of the catalytic beta-subunits act on the gamma-subunit through repulsive van der Waals interactions to generate a torque that drives unidirectional rotation, as observed experimentally. The simulations suggest that the calculated 85 degrees substep rotation is driven primarily by ATP binding and that the subsequent 35 degrees substep rotation is produced by product release from one beta-subunit and a concomitant binding pocket expansion of another beta-subunit. The results of the simulation agree with single-molecule experiments [see, for example, Adachi K, et al. (2007) Cell 130:309-321] and support a tri-site rotary mechanism for F(1)-ATPase under physiological condition.

Nicotinic acetylcholine receptors (nAChR) are pentameric ligand-gated ion channels composed of subunits that consist of an extracellular domain that carries the ligand-binding site and a distinct ion-pore domain. Signal transduction results from the allosteric coupling between the two domains: the distance from the binding site to the gate of the pore domain is 50 A. Normal mode analysis with a C(alpha) Gaussian network of a new structural model of the neuronal alpha7 nAChR showed that the lowest mode involves a global quaternary twist motion that opens the ion pore. A molecular probe analysis, in which the network is modified at each individual amino acid residue, demonstrated that the major effect is to change the frequency, but not the form, of the twist mode. The largest effects were observed for the ligand-binding site and the Cys-loop. Most (24/27) of spontaneous mutations known to cause congenital myasthenia and autosomal dominant nocturnal frontal lobe epilepsy are located either at the interface between subunits or, within a given subunit, at the interface between rigid blocks. These interfaces are modified significantly by the twist mode. The present analysis, thus, supports the quaternary twist model of the nAChR allosteric transition and provides a qualitative interpretation of the effect of the mutations responsible for several receptor pathologies.

An understanding of the thermodynamics and kinetics of protein folding requires a knowledge of the free energy surface governing the motion of the polypeptide chain. Because of the many degrees of freedom involved, surfaces projected on only one or two progress variables are generally used in descriptions of the folding reaction. Such projections result in relatively smooth surfaces, but they could mask the complexity of the unprojected surface. Here we introduce an approach to determine the actual (unprojected) free energy surface and apply it to the second beta-hairpin of protein G, which has been used as a model system for protein folding. The surface is represented by a disconnectivity graph calculated from a long equilibrium folding-unfolding trajectory. The denatured state is found to have multiple low free energy basins. Nevertheless, the peptide shows exponential kinetics in folding to the native basin. Projected surfaces obtained from the present analysis have a simple form in agreement with other studies of the beta-hairpin. The hidden complexity found for the beta-hairpin surface suggests that the standard funnel picture of protein folding should be revisited.

Amyloid fibril formation involves nonfibrillar oligomeric intermediates, which are important as possible cytotoxic species in neurodegenerative diseases. However, their transient nature and polydispersity have made it difficult to identify their formation mechanism or structure. We have investigated the dimerization process, the first step in aggregate formation, by multiple molecular dynamics simulations of five beta-sheet-forming peptides. Contrary to the regular beta-sheet structure of the amyloid fibril, the dimers exhibit all possible combinations of beta-sheets, with an overall preference for antiparallel arrangements. Through statistical analysis of 1,000 dimerization trajectories, each 1 ns in length, we have demonstrated that the observed distribution of dimer configurations is kinetically determined; hydrophobic interactions orient the peptides so as to minimize the solvent accessible surface area, and the dimer structures become trapped in energetically unfavorable conformations. Once the hydrophobic contacts are present, the backbone hydrogen bonds form rapidly by a zipper-like mechanism. The initial nonequilibrium structures formed are stable during the 1-ns simulation time for all five peptides at room temperature. In contrast, at higher temperatures, where rapid equilibration among different configurations occurs, the distribution follows the global energies. The relaxation time of dimers at room temperature was estimated to be longer than the time for diffusional encounters with other oligomers at typical concentrations. These results suggest that kinetic trapping could play a role in the structural evolution of early aggregates in amyloid fibrillogenesis., (Copyright 2004 The National Academy of Sciencs of the USA)

The ability of proteins to fold to well defined compact structures is one of the most remarkable examples of the effect of natural selection on biological molecules. To understand their properties, including the stability, the mechanism of folding, and the possibilities of misfolding and association, it is necessary to know the protein free energy landscape. We use NMR data as restraints in a Monte Carlo sampling procedure to determine the ensemble of structures populated by human alpha-lactalbumin in the presence of increasing concentrations of urea. The ensembles of structures that represent the partially folded states of the protein show that two structural cores, corresponding to portions of the alpha and beta domains of the native protein, are preserved even when the native-like interactions that define their existence are substantially weakened. Analysis of the network of residual contacts reveals the presence of a complex interface region between the two structural cores and indicates that the development of specific interactions within this interface is the key step in achieving the native structure. The relative probabilities of the conformations determined from the NMR data are used to construct a coarse-grained free energy landscape for alpha-lactalbumin in the absence of urea. The form of the landscape, together with the existence of distinct cores, supports the concept that robustness and modularity are the properties that make possible the folding of complex proteins.

Although the binding change mechanism of rotary catalysis by which F1-ATPase hydrolyzes ATP has been supported by equilibrium, kinetic, and structural observations, many questions concerning the function remain unanswered. Because of the importance of this enzyme, the search for a full understanding of its mechanism is a key problem in structural biology. Making use of the results of free energy simulations and experimental binding constant measurements, a model is developed for the free energy change during the hydrolysis cycle. This model makes possible the development of a kinetic scheme for ATP hydrolysis by F1-ATPase, in which the rate constants are associated with specific configurations of the beta subunits. An essential new element is that the strong binding site for ADP,Pi is shown to be the betaDP site, in contrast to the strong binding site for ATP, which is betaTP. This result provides a rationale for the rotation of the gamma subunit, which induces the cooperativity required for a tri-site binding change mechanism. The model explains a series of experimental data, including the ATP concentration dependence of the rate of hydrolysis and catalytic site occupation for both the Escherichia coli F1-ATPase (EcF1) and Thermophilic Bacillus PS3 F1-ATPase (TF1), which have different behavior.

We present a general approach in which theory and experiments are combined in an iterative manner to provide a detailed description of the transition state ensemble (TSE) for folding. The method is illustrated by applying it to TNfn3, a fibronectin type III domain protein. In the first iteration, a coarse-grained determination of the TSE is carried out by using a limited set of experimental phi values as constraints in a molecular dynamics sampling simulation. The resulting model of the TSE is used to determine the additional residues whose phi value measurement would provide the most information for refining the TSE. Successive iterations with an increasing number of phi value measurements are carried out until no further changes in the properties of the TSE are detected or there are no additional residues whose phi values can be measured. In the study of TNfn3 three iterations were necessary to achieve self-consistency. A retrospective application of the method can be used to determine the accuracy of the TSE results and to find "key residues" for folding, i.e., those that are most important for the formation of the TSE. The approach reported here is an efficient method for finding the structures that make up the TSEs for protein folding. Its use will improve future efforts for their experimental determination and refinement.

Chorismate mutase acts at the first branch-point of aromatic amino acid biosynthesis and catalyzes the conversion of chorismate to prephenate. The results of molecular dynamics simulations of the substrate in solution and in the active site of chorismate mutase are reported. Two nonreactive conformers of chorismate are found to be more stable than the reactive pseudodiaxial chair conformer in solution. It is shown by QM/MM molecular dynamics simulations, which take into account the motions of the enzyme, that when these inactive conformers are bound to the active site, they are rapidly converted to the reactive chair conformer. This result suggests that one contribution of the enzyme is to bind the more prevalent nonreactive conformers and transform them into the active form in a step before the chemical reaction. The motion of the reactive chair conformer in the active site calculated by using the QM/MM potential generates transient structures that are closer to the transition state than is the stable CHAIR conformer.

Unfolding of proteins by forced stretching with atomic force microscopy or laser tweezer experiments complements more classical techniques using chemical denaturants or temperature. Forced unfolding is of particular interest for proteins that are under mechanical stress in their biological function. For beta-sandwich proteins (a fibronectin type III and an immunoglobulin domain), both of which appear in the muscle protein titin, the results of stretching simulations show important differences from temperature-induced unfolding, but there are common features that point to the existence of folding cores. Intermediates detected by comparing unfolding with a biasing perturbation and a constant pulling force are not evident in temperature-induced unfolding. For an alpha-helical domain (alpha-spectrin), which forms part of the cytoskeleton, there is little commonality in the pathways from unfolding induced by stretching and temperature. Comparison of the forced unfolding of the two beta-sandwich proteins and two alpha-helical proteins (the alpha-spectrin domain and an acyl-coenzyme A-binding protein) highlights important differences within and between protein classes that are related to the folding topologies and the relative stability of the various structural elements.

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.

The isomerization of chorismate to prephenate by chorismate mutase in the biosynthetic pathway that forms Tyr and Phe involves C5---O (ether) bond cleavage and C1---C9 bond formation in a Claisen rearrangement. Development of negative charge on the ether oxygen, stabilized by Lys-168 and Glu-246, is inferred from the structure of a complex with a transition state analogue (TSA) and from the pH-rate profile of the enzyme and the E246Q mutant. These studies imply a protonated Glu-246 well above pH 7. Here, several 500-ps molecular dynamics simulations test the stability of enzyme-TSA complexes by using a solvated system with stochastic boundary conditions. The simulated systems are (i) protonated Glu-246 (stable), (ii) deprotonated Glu-246 (unstable), (iii) deprotonated Glu-246 plus one H2O between Glu-246 and the ether oxygen (unstable), (iv) the E246Q mutant (stable), and (v) addition of OH- between protonated Glu-246 and the ether oxygen. In (v), a local conformational change of Lys-168 displaced the OH- into the solvent region, suggesting a possible rate-determining step that precedes the catalytic step. In a 500-ps simulation of the enzyme complexed with the reactant chorismate or the product prephenate, no water molecule remained near the oxygen of the ligand. Calculations using the linearized Poisson-Boltzmann equation show that the effective pKa of Glu-246 is shifted from 5.8 to 8.1 as the negative charge on the ether oxygen of the TSA is changed from -0.56 electron to -0.9 electron. Altogether, these results support retention of a proton on Glu-246 to high pH and the absence of a water molecule in the catalytic steps.

Normal mode calculations on individual subunits and a multisubunit construct are used to analyze the structural transitions that occur during the GroEL cycle. The normal modes demonstrate that the specific displacements of the domains (hinge bending, twisting) observed in the structural studies arise from the intrinsic flexibility of the subunits. The allosteric mechanism (positive cooperativity within a ring, negative cooperativity between rings) is shown to be based on coupled tertiary structural changes, rather than the quaternary transition found in classic allosteric proteins. The results unify static structural data from x-ray crystallography and cryoelectron microscopy with functional measurements of binding and cooperativity.

The calculated folding thermodynamics of a simple off-lattice three-helix-bundle protein model under equilibrium conditions shows the experimentally observed protein transitions: a collapse transition, a disordered-to-ordered globule transition, a globule to native-state transition, and the transition from the active native state to a frozen inactive state. The cooperativity and physical origin of the various transitions are explored with a single "optimization" parameter and characterized with the Lindemann criterion for liquid versus solid-state dynamics. Below the folding temperature, the model has a simple free energy surface with a single basin near the native state; the surface is similar to that calculated from a simulation of the same three-helix-bundle protein with an all-atom representation [Boczko, E. M. & Brooks III, C. L. (1995) Science 269, 393-396].

Conformational changes in ras p21 triggered by the hydrolysis of GTP play an essential role in the signal transduction pathway. The path for the conformational change is determined by molecular dynamics simulation with a holonomic constraint directing the system from the known GTP-bound structure (with the gamma-phosphate removed) to the GDP-bound structure. The simulation is done with a shell of water molecules surrounding the protein. In the switch I region, the side chain of Tyr-32, which undergoes a large displacement, moves through the space between loop 2 and the rest of the protein, rather than on the outside of the protein. As a result, the charged residues Glu-31 and Asp-33, which interact with Raf in the homologous RafRBD-Raps complex, remain exposed during the transition. In the switch II region, the conformational changes of alpha2 and loop 4 are strongly coupled. A transient hydrogen bonding complex between Arg-68 and Tyr-71 in the switch II region and Glu-37 in switch I region stabilizes the intermediate conformation of alpha2 and facilitates the unwinding of a helical turn of alpha2 (residues 66-69), which in turn permits the larger scale motion of loop 4. Hydrogen bond exchange between the protein and solvent molecules is found to be important in the transition. Possible functional implications of the results are discussed.

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.

The transition in barnase from the native state to a compact globule has been studied with high-temperature molecular dynamics simulations. A partial destruction of the alpha-helices and the outer strands of the beta-sheet is observed with water molecules replacing the hydrogen bonds of the secondary structural elements. Simultaneously, the main alpha-helix moves away from the beta-sheet and exposes the principal hydrophobic core, many of whose nonpolar side chains, beginning with the ones near the surface, become solvated by hydrogen-bonded water molecules. This step involves a significant increase in the solvent-exposed surface area; the resulting loss of stability due to the hydrophobic effect may be the major source of the activation barrier in the unfolding reaction. The detailed mechanism described here for the first stage of the denaturation of barnase, including the essential role of water molecules, is likely to be representative of protein denaturation, in general.

Molecular dynamics simulations have been used to compute the difference in the unfolding free energy between wild-type barnase and the mutant in which Ile-96 is replaced by alanine. The simulations yield results (-3.42 and -5.21 kcal/mol) that compare favorably with experimental values (-3.3 and -4.0 kcal/mol). The major contributions to the free energy difference arise from bonding terms involving degrees of freedom of the mutated side chain and from nonbonded interactions of that side chain with its environment in the folded protein. By comparison with simulations of an extended peptide in the absence of solvent, used as a reference state, hydration effects are shown to play a minor role in the overall free energy balance for the Ile----Ala transformation. The implications of these results for our understanding of the hydrophobic effect and its contribution to protein stability are discussed.

Sequences of the principal neutralizing determinant (PND) of the external envelope protein, gp120, from 245 isolates of human immunodeficiency virus type 1 are analyzed. The minimal set of peptides that would elicit antibodies to neutralize a majority of U.S. and European isolates of human immunodeficiency virus type 1 is determined with the assumption that peptides of a given length including the central Gly-Pro-Gly triad are required. In spite of the hypervariability of the PND, 90% of these 245 sequences include peptides from a set of 7 pentapeptides, 13 hexapeptides, or 17 heptapeptides. Tests of these peptide sets on 78 additional PND sequences show that 95% are covered by the 7 pentapeptides, 94% by the 13 hexapeptides, and 86% by the 17 heptapeptides. To anticipate variants not yet observed, single amino acid mutation frequencies from the 245 isolates are used to calculate an expanded set of the 10,000 most probable PND sequences. These sequences cover 86% of the total distribution expected for the central portion of the PND. Peptide lists derived from this expanded set when tested on the 78 additional sequences show that 7 pentapeptides cover 95%, 13 hexapeptides cover 94%, and 17 heptapeptides cover 94%. These results suggest that peptide cocktails of limited size with the potential to cover a large fraction of PND sequence variation may be feasible vaccine candidates.

A detailed model for the structure and dynamics of the interior of the lipid bilayer in the liquid crystal phase is presented. The model includes two classes of motion: (i) the internal dynamics of the chains, determined from Brownian dynamics simulations with a continuous version of the Marcelja mean-field potential, and (ii) noncollective reorientation (axial rotation and wobble) of the entire molecule, introduced by a cone model. The basic unit of the model is a single lipid chain with field parameters adjusted to fit the 2H order parameters and the frequency-dependent 13C NMR T1 relaxation times of dipalmitoyl phosphatidylcholine bilayers. The chain configurations obtained from the trajectory are used to construct a representation of the bilayer. The resulting lipid assembly is consistent with NMR, neutron diffraction, surface area, and density data. It indicates that a high degree of chain disorder and entanglement exists in biological membranes.

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.

Molecular dynamics simulations are used to calculate the incoherent neutron scattering spectra of myoglobin between 80 K and 325 K and compared with experimental data. There is good agreement over the entire temperature range for the elastic, quasi-elastic, and inelastic components of the scattering. This provides support for the accuracy of the simulations of the internal motions that make the dominant contributions to the atomic displacements on a time scale of 0.3-100 ps (100-0.3 cm-1). Analysis of the simulations shows that at low temperatures a harmonic description of the molecule is appropriate and that the molecule is trapped in localized regions of conformational space. At higher temperatures the scattering arises from a combination of vibrations within wells (substates) and transitions between them; the latter contribute to the quasi-elastic scattering.

A normal mode analysis making use of an empirical potential function including local and nonlocal (nonbonded) interactions is performed for the bovine pancreatic trypsin inhibitor in the full conformational space of the molecule (1,740 degrees of freedom); that is, all bond lengths and angles, as well as dihedral angles, are included for the 580-atom system consisting of all heavy atoms and polar hydrogens. The heavy-atom frequency spectrum shows a dense distribution between 3 and 1,800 cm-1, with 350 modes below 216 cm-1. Most of the low-frequency modes, of which many have significant anharmonic character, are found to be delocalized over the protein. The root-mean-square amplitudes of the atomic fluctuations are calculated at 300 K from the normal modes and compared with those obtained from a solution molecular dynamics simulation based on the same potential function; very good agreement is obtained for the variation in the main-chain fluctuations as a function of residue number, though larger differences occur for the side chains. The fluctuations are generally, though not always, dominated by frequencies below 30 cm-1, in accord with the results of the dynamics simulation. The vibrational contributions to the thermodynamic properties of the protein are calculated as a function of temperature; the effects of perturbations on the spectrum, suggested for ligand or substrate binding, are examined. The analysis demonstrates that, in spite of the anharmonic contributions to the potential, a normal mode description can provide useful results concerning the internal motions of proteins.

The rapid fluctuations of protein atoms derived from molecular dynamics simulations can be extrapolated to longer-time motions by effective single-particle stochastic models. This is demonstrated by an analysis of velocity autocorrelation functions for the atoms of lysine side chains in the active site of RNase A. The atomic motions are described by a bounded stochastic model with the friction and noise parameters determined from a molecular dynamics simulation. The low-frequency relaxation behavior is shown to result from collisional damping rather than dephasing. Extrapolation of these results to the quasistochastic motion of the heme group in myoglobin provides an explanation of 57Fe Mössbauer spectroscopic data.

A reaction path is presented by which the effects of oxygen binding in hemoglobin are transmitted from a heme group to the surface of its subunit. Starting from the known deoxy geometry, it is shown by calculations with empirical energy functions and comparisons with available data how the change in heme geometry on ligation introduces a perturbation that leads to the tertiary structural alterations essential for cooperatively. It is found that there is little strain on the unliganded heme; instead, the reduced oxygen affinity of hemoglobin results from the strain on the liganded subunit in a tetramer with the deoxy quarternary structure.

The solution structure and dynamics of the 17-base-pair synthetic operator site OL1, which is recognized by the cI and Cro repressors of bacteriophage lambda, is studied by two-dimensional NMR methods. A sequential assignment strategy for nucleic acids is proposed and illustrated by the assignments of the base and sugar protons of OL1.

Model calculations are performed to evaluate the utility of molecular dynamics with NMR interproton distance restraints for determining the three-dimensional structure of proteins. The system used for testing the method is the 1.5-A resolution crystal structure of crambin (a protein of 46 residues) from which a set of 240 approximate interproton distances of less than 4 A are derived. The convergence properties of the method are examined by using different dynamics protocols and starting from two initial structures; one is a completely extended beta-strand, and the other has residues 7-19 and 23-30 in the form of alpha-helices (as in the crystal structure) with the remaining residues in the form of extended beta-strands. In both cases global and local convergence to the correct final structure is achieved with rms atomic differences between the restrained dynamics structures and the crystal structure of 1.5-2.1 A for the backbone atoms and 2.1-2.8 A for all atoms; the averaged structure has backbone and all atom rms deviations of 1.3 and 1.9 A, respectively. Further, it is shown that a restrained dynamics structure with significantly larger deviations (i.e., 5.7 A for the backbone atoms) can be characterized as incorrect, independent of a knowledge of the crystal structure.

A method is presented for protein secondary structure prediction based on a neural network. A training phase was used to teach the network to recognize the relation between secondary structure and amino acid sequences on a sample set of 48 proteins of known structure. On a separate test set of 14 proteins of known structure, the method achieved a maximum overall predictive accuracy of 63% for three states: helix, sheet, and coil. A numerical measure of helix and sheet tendency for each residue was obtained from the calculations. When predictions were filtered to include only the strongest 31% of predictions, the predictive accuracy rose to 79%.

A large number of mutant proteins with single amino acid substitutions are now being produced. The ability to predict the structural changes expected from such mutations would aid greatly in the efficient utilization of the mutagenic techniques and in the interpretation of the changes in stability and function that result. A minimum perturbation approach is suggested as a first step in such structural predictions and is tested by application to a recently isolated variant of the hemagglutinin glycoprotein. The agreement between the predicted structure and that inferred from the x-ray refinement is encouraging and provides support for the proposed modeling procedure.

A method is presented for finding particular normal modes for large molecules such as proteins and nucleic acids. The method is based on an iterative approach that extracts eigenvectors of interest from the full second-derivative matrix. Application of the method to the interdomain (hinge-bending) motion of lysozyme yields a frequency of 3.6 cm-1. This is similar to the value obtained from earlier adiabatic-energy-minimization studies. Analysis of the mode shows that the relaxation associated with the hinge bending is highly delocalized; that is, the dihedral angle and energy changes are distributed over many residues, including some (e.g., Trp-28) that are distant from the cleft and hinge region.

The mean square amplitudes of atomic fluctuations for a polypeptide (decaglycine) alpha-helix evaluated from molecular dynamics simulations at seven temperatures between 5 and 300 K are compared with analytic harmonic results and with experimental values. Above 100 K the harmonic approximation significantly underestimates the amplitudes of the displacements. Analysis of the time dependence of the fluctuations shows that low-frequency modes (<75 cm(-1)) dominate the atomic fluctuations and that there is a contribution with a very long relaxation time (>10 ps). Quantum corrections to the amplitude of the fluctuations are found to be small above 50 K. The mean square amplitudes obtained from the molecular dynamics simulations are compared with the values derived from x-ray temperature (Debye-Waller) factors for metmyoglobin (80, 250, and 300 K) and ferrocytochrome c (300 K).

A method is described for calculating the reaction rate in globular proteins of activated processes such as ligand binding or enzymatic catalysis. The method is based on the determination of the probability that the system is in the transition state and of the magnitude of the reactive flux for transition-state systems. An "umbrella sampling" simulation procedure is outlined for evaluating the transition-state probability. The reactive flux is obtained from an approach described previously for calculating the dynamics of transition-state trajectories. An application to the rotational isomerization of an aromatic ring in the bovine pancreatic trypsin inhibitor is presented. The results demonstrate the feasibility of calculating rate constants for reactions in proteins and point to the importance of solvent effects for reactions that occur near the protein surface.

The stochastic boundary molecular dynamics method is used to study the structure, dynamics, and energetics of the solvated active site of bovine pancreatic ribonuclease A. Simulations of the native enzyme and of the enzyme complexed with the dinucleotide substrate CpA and the transition-state analog uridine vanadate are compared. Structural features and dynamical couplings for ribonuclease residues found in the simulation are consistent with experimental data. Water molecules, most of which are not observed in crystallographic studies, are shown to play an important role in the active site. Hydrogen bonding of residues with water molecules in the free enzyme is found to mimic the substrate-enzyme interactions of residues involved in binding. Networks of water stabilize the cluster of positively charged active site residues. Correlated fluctuations between the uridine vanadate complex and the distant lysine residues are mediated through water and may indicate a possible role for these residues in stabilizing the transition state.