Factors governing fibrillogenesis of polypeptide chains revealed by lattice models.

Using lattice models we explore the factors that determine the tendencies of polypeptide chains to aggregate by exhaustively sampling the sequence and conformational space. The morphologies of the fibril-like structures and the time scales (τ(fib)) for their formation depend on a balance between hydrophobic and Coulomb interactions. The extent of population of an ensemble of N* structures, which are fibril-prone structures in the spectrum of conformations of an isolated protein, is the major determinant of τ(fib). This observation is used to determine the aggregating sequences by exhaustively exploring the sequence space, thus providing a basis for genome wide search of fragments that are aggregation prone.

Probing the mechanisms of fibril formation using lattice models.

Using exhaustive Monte Carlo simulations we study the kinetics and mechanism of fibril formation using lattice models as a function of temperature (T) and the number of chains (M). While these models are, at best, caricatures of peptides, we show that a number of generic features thought to govern fibril assembly are captured by the toy model. The monomer, which contains eight beads made from three letters (hydrophobic, polar, and charged), adopts a compact conformation in the native state. In both the single-layered protofilament (seen for M10) structures, the monomers are arranged in an antiparallel fashion with the "strandlike" conformation that is perpendicular to the fibril axis. Partial unfolding of the folded monomer that populates an aggregation prone conformation (N(*)) is required for ordered assembly. The contacts in the N(*) conformation, which is one of the four structures in the first "excited" state of the monomer, are also present in the native conformation. The time scale for fibril formation is a minimum in the T-range when the conformation N(*) is substantially populated. The kinetics of fibril assembly occurs in three distinct stages. In each stage there is a cascade of events that transforms the monomers and oligomers to ordered structures. In the first "burst" stage, highly mobile oligomers of varying sizes form. The conversion to the N(*) conformation occurs within the oligomers during the second stage in which a vast number of interchain contacts are established. As time progresses, a dominant cluster emerges that contains a majority of the chains. In the final stage, the aggregation of N(*) particles serve as a template onto which smaller oligomers or monomers can dock and undergo conversion to fibril structures. The overall time for growth in the latter stages is well described by the Lifshitz-Slyazov growth kinetics for crystallization from supersaturated solutions. The detailed analysis shows that elements of the three popular models, namely, nucleation and growth, templated assembly, and nucleated conformational conversion are present at various stages of fibril assembly.

Dissecting contact potentials for proteins: relative contributions of individual amino acids.

Knowledge-based contact potentials are routinely used in fold recognition, binding of peptides to proteins, structure prediction, and coarse-grained models to probe protein folding kinetics. The dominant physical forces embodied in the contact potentials are revealed by eigenvalue analysis of the matrices, whose elements describe the strengths of interaction between amino acid side chains. We propose a general method to rank quantitatively the importance of various inter-residue interactions represented in the currently popular pair contact potentials. Eigenvalue analysis and correlation diagrams are used to rank the inter-residue pair interactions with respect to the magnitude of their relative contributions to the contact potentials. The amino acid ranking is shown to be consistent with a mean field approximation that is used to reconstruct the original contact potentials from the most relevant amino acids for several contact potentials. By providing a general, relative ranking score for amino acids, this method permits a detailed, quantitative comparison of various contact interaction schemes. For most contact potentials, between 7 and 9 amino acids of varying chemical character are needed to accurately reconstruct the full matrix. By correlating the identified important amino acid residues in contact potentials and analysis of about 7800 structural domains in the CATH database we predict that it is important to model accurately interactions between small hydrophobic residues. In addition, only potentials that take interactions involving the protein backbone into account can predict dense packing in protein structures.

Continuous anisotropic representation of coarse-grained potentials for proteins by spherical harmonics synthesis.

A new method is presented for extracting statistical potentials dependent on the relative side chain and backbone orientations in proteins. Coarse-grained, anisotropic potentials are constructed for short-, medium-, and long-range interactions using the Boltzmann method and a database of non-homologous protein structures. The new orientation-dependent potentials are analyzed using a spherical harmonics decomposition method with real eigenfunctions. This method permits a more realistic, continuous angular representation of the coarse-grained potentials. Results of tests for discriminating the native protein conformations from large sets of decoy proteins, show that the new continuous distance- and orientation-dependent potentials present significantly improved performance. Novel graphical representations are developed and used to depict the orientational dependence of the interaction potentials. These new continuous anisotropic statistical potentials could be instrumental in developing new computational methods for structure prediction, threading and coarse-grained simulations.

Development of novel statistical potentials for protein fold recognition.

The need to perform large-scale studies of protein fold recognition, structure prediction and protein-protein interactions has led to novel developments of residue-level minimal models of proteins. A minimum requirement for useful protein force-fields is that they be successful in the recognition of native conformations. The balance between the level of detail in describing the specific interactions within proteins and the accuracy obtained using minimal protein models is the focus of many current protein studies. Recent results suggest that the introduction of explicit orientation dependence in a coarse-grained, residue-level model improves the ability of inter-residue potentials to recognize the native state. New statistical and optimization computational algorithms can be used to obtain accurate residue-dependent potentials for use in protein fold recognition and, more importantly, structure prediction.

Enhanced sampling in numerical path integration: an approximation for the quantum statistical density matrix based on the nonextensive thermostatistics.

Here we examine a proposed approximation for the quantum statistical density matrix motivated by the nonextensive thermostatistics of Tsallis and co-workers. The approximation involves replacing the physical potential energy with an effective one, corresponding to a generalized nonextensive statistical ensemble. We examine the convergence properties of averages calculated using the effective potential, and introduce a related method for enhanced sampling in numerical path integration. As a necessary measure, path integral energy estimators are introduced for potentials that involve explicit temperature dependence. This sampling method is found to be effective for path integral simulations involving broken ergodicity.

Probing the origins of increased activity of the E22Q "Dutch" mutant Alzheimer's beta-amyloid peptide.

The amyloid peptide congener A beta(10--35)-NH(2) is simulated in an aqueous environment in both the wild type (WT) and E22Q "Dutch" mutant forms. The origin of the noted increase in deposition activity resulting from the Dutch mutation is investigated. Multiple nanosecond time scale molecular dynamics trajectories were performed and analyzed using a variety of measures of the peptide's average structure, hydration, conformational fluctuations, and dynamics. The results of the study support the conclusions that 1) the E22Q mutant and WT peptide are both stable in "collapsed coil" conformations consistent with the WT structure of, J. Struct. Biol. 130:130--141); 2) the E22Q peptide is more flexible in solution, supporting early claims that its equilibrium structural fluctuations are larger than those of the WT peptide; and 3) the local E22Q mutation leads to a change in the first solvation layer in the region of the peptide's "hydrophobic patch," resulting in a more hydrophobic solvation of the mutant peptide. The simulation results support the view that the noted increase in activity due to the Dutch mutation results from an enhancement of the desolvation process that is an essential step in the aggregation of the peptide.

Simulation study of the structure and dynamics of the Alzheimer's amyloid peptide congener in solution.

The amyloid Abeta(10-35)-NH2 peptide is simulated in an aqueous environment on the nanosecond time scale. One focus of the study is on the validation of the computational model through a direct comparison of simulated statistical averages with experimental observations of the peptide's structure and dynamics. These measures include (1) nuclear magnetic resonance spectroscopy-derived amide bond order parameters and temperature-dependent H(alpha) proton chemical shifts, (2) the peptide's radius of gyration and end-to-end distance, (3) the rates of peptide self-diffusion in water, and (4) the peptide's hydrodynamic radius as measured by quasielastic light scattering experiments. A second focus of the study is the identification of key intrapeptide interactions that stabilize the central structural motif of the peptide. Particular attention is paid to the structure and fluctuation of the central LVFFA hydrophobic cluster (17-21) region and the VGSN turn (24-27) region. There is a strong correlation between preservation of the structure of these elements and interactions between the cluster and turn regions in imposing structure on the peptide monomer. The specific role of these interactions in relation to proposed mechanisms of amyloidosis is discussed.

Energy landscape theory for Alzheimer's amyloid beta-peptide fibril elongation.

Recent experiments on the kinetics of deposition and fibril elongation of the Alzheimer's beta-amyloid peptide on preexisting fibrils are analyzed. A mechanism is developed based on the dock-and-lock scheme recently proposed by Maggio and coworkers to organize their experimental observations of the kinetics of deposition of beta-peptide on preexisting amyloid fibrils and deposits. Our mechanism includes channels for (1) a one-step prion-like direct deposition on fibrils of activated monomeric peptide in solution, and (2) a two-step deposition of unactivated peptide on fibrils and subsequent reorganization of the peptide-fibril complex. In this way, the mechanism and implied "energy landscape" unify a number of schemes proposed to describe the process of fibril elongation. This beta-amyloid landscape mechanism (beta ALM) is found to be in good agreement with existing experimental data. A number of experimental tests of the mechanism are proposed. The mechanism leads to a clear definition of overall equilibrium or rate constants in terms of the energetics of the elementary underlying processes. Analysis of existing experimental data suggests that fibril elongation occurs through a two-step mechanism of nonspecific peptide absorption and reorganization. The mechanism predicts a turnover in the rate of fibril elongation as a function of temperature and denaturant concentration. Proteins 2001;42:217-229.

Vibrational population relaxation of carbon monoxide in the heme pocket of photolyzed carbonmonoxy myoglobin: comparison of time-resolved mid-IR absorbance experiments and molecular dynamics simulations.

The vibrational energy relaxation of carbon monoxide in the heme pocket of sperm whale myoglobin was studied by using molecular dynamics simulation and normal mode analysis methods. Molecular dynamics trajectories of solvated myoglobin were run at 300 K for both the delta- and epsilon-tautomers of the distal His-64. Vibrational population relaxation times of 335 +/- 115 ps for the delta-tautomer and 640 +/- 185 ps for the epsilon-tautomer were estimated by using the Landau-Teller model. Normal mode analysis was used to identify those protein residues that act as the primary "doorway" modes in the vibrational relaxation of the oscillator. Although the CO relaxation rates in both the epsilon- and delta-tautomers are similar in magnitude, the simulations predict that the vibrational relaxation of the CO is faster in the delta-tautomer with the distal His playing an important role in the energy relaxation mechanism. Time-resolved mid-IR absorbance measurements were performed on photolyzed carbonmonoxy hemoglobin (Hb(13)CO). From these measurements, a T(1) time of 600 +/- 150 ps was determined. The simulation and experimental estimates are compared and discussed.

A study of vibrational relaxation of B-state carbon monoxide in the heme pocket of photolyzed carboxymyoglobin.

The vibrational energy relaxation of dissociated carbon monoxide in the heme pocket of sperm whale myoglobin has been studied using equilibrium molecular dynamics simulation and normal mode analysis methods. Molecular dynamics trajectories of solvated myoglobin were run at 300 K for both the delta- and epsilon-tautomers of the distal histidine, His64. Vibrational population relaxation times were estimated using the Landau-Teller model. For carbon monoxide (CO) in the myoglobin epsilon-tautomer, for a frequency of omega0 = 2131 cm-1 corresponding to the B1 state, T1epsilon(B1) = 640 +/- 185 ps, and for a frequency of omega0 = 2119 cm-1 corresponding to the B2 state, T1epsilon(B2) = 590 +/- 175 ps. Although the CO relaxation rates in both the epsilon- and delta-tautomers are similar in magnitude, the simulations predict that the vibrational relaxation of the CO is faster in the delta-tautomer. For CO in the myoglobin delta-tautomer, it was found that the relaxation times were identical within error for the two CO substate frequencies, T1delta(B1) = 335 +/- 115 ps and T1delta(B2) = 330 +/- 145 ps. These simulation results are in reasonable agreement with experimental results of Anfinrud and coworkers (unpublished results). Normal mode calculations were used to identify the dominant coupling between the protein and CO molecules. The calculations suggest that the residues of the myoglobin pocket, acting as a first solvation shell to the CO molecule, contribute the primary "doorway" modes in the vibrational relaxation of the oscillator.

Direct computation of long time processes in peptides and proteins: reaction path study of the coil-to-helix transition in polyalanine.

The MaxFlux reaction path algorithm was used to isolate optimal transition pathways for the coil-to-helix transition in polyalanine. Eighteen transition pathways, each connecting one random coil configuration with an ideal alpha-helical configuration, were computed and analyzed. The transition pathway energetics and mechanism were analyzed in terms of the progression of the peptide nonbonded contact formation, helicity, end-to-end distance and energetics. It was found that (1) localized turns characterized by i, i + 3 hydrogen bonds form in the early stages of the coil-to-helix transition, (2) the peptide first collapses and then becomes somewhat more extended in the final stage of helix formation, and (3) 310-helix formation does not appear to be a necessary step in the transition from coil to helix. These conclusions are in agreement with the results of more computationally intensive direct molecular dynamics simulations. Proteins 1999;36:249-261.

Novel methods of sampling phase space in the simulation of biological systems.

New advances in molecular dynamics and Monte Carlo simulations have led to impressive increases in the speed of sampling phase space for complex biological systems. These methods have been combined with new fast algorithms for computing long range electrostatic interactions for new polarizable force fields. In addition, new methods for sampling low energy molecular conformations allow the rapid determination of thermodynamically dominant regions on the potential-energy surface. Accurate measures of the rate of phase-space sampling should allow both the optimization and the comparison of methods for a particular problem of interest.

Molecular dynamics simulation of NO recombination to myoglobin mutants.

Molecular dynamics simulations on two coupled electronic surfaces are employed to investigate the geminate recombination of nitric oxide to mutants of sperm whale myoglobin. A model for the ground and the excited states is constructed based on experimental data. The crossing between the surfaces is treated using the Landau-Zener formula. The reaction probability and the recombination curves are calculated directly by histogramming the results of an ensemble of trajectories. The experimental trend is reproduced in which the picosecond recombination rate of different mutants increases in the order Phe29 > Leu29 > Val29 > Ala29. Furthermore, in accord with the experiment on significantly longer time scales an opposite trend is obtained, in which the recombination rate for Ala29 is larger than for Phe29. These results are explained by constrained diffusion of the ligand in the heme pocket. The average and the transient volume of the heme pocket is modified by the 29 mutants.

Theoretical probes of conformational fluctuations in S-peptide and RNase A/3'-UMP enzyme product complex.

The dynamic properties of the RNase A/3'-UMP enzyme/product complex and the S-peptide of RNase A have been investigated by molecular dynamics simulations using suitable generalization of ideas introduced to probe the energy landscape in structural glasses. We introduce two measures, namely, the kinetic energy fluctuation metric and the force metric, both of which are used to calculate the time needed for sampling the conformation space of the molecules. The calculation of the fluctuation metric requires a single trajectory whereas the force metric is computed using two independent trajectories. The vacuum MD simulations show that for both systems the time required for kinetic energy equipartitioning is surprisingly long even at high temperatures. We show that the force metric is a powerful means of probing the nature and relative importance of conformational substates which determine the dynamics at low temperatures. In particular the time dependence of the non-bonded force metric is used to demonstrate that at low temperatures the system is predominantly localized in a single cluster of conformational substates. The force metric is used to show that relaxation of long range (in sequence space) interactions must be mediated by a sequence of local dihedral angle transitions. We also argue that the time needed for compact structure formation is intimately related to the time needed for the relaxation of the dihedral angle degrees of freedom. The time for non-bonded interactions, which drive protein molecules to fold under appropriate conditions, to relax becomes extremely long as the temperature is lowered suggesting that the formation of maximally compact structure in proteins must be a very slow process.

Exploring the energy landscape in proteins.

We present two methods to probe the energy landscape and motions of proteins in the context of molecular dynamics simulations of the helix-forming S-peptide of RNase A and the RNase A-3'-UMP enzyme-product complex. The first method uses the generalized ergodic measure to compute the rate of conformational space sampling. Using the dynamics of nonbonded forces as a means of probing the time scale for ergodicity to be obtained, we argue that even in a relatively short time (< 10 psec) several different conformational substrates are sampled. At longer times, barriers on the order of a few kcal/mol (1 cal = 4.184 J) are involved in the large-scale motion of proteins. We also present an approximate method for evaluating the distribution of barrier heights g(EB) using the instantaneous normal-mode spectra of a protein. For the S-peptide, we show that g(EB) is adequately represented by a Poisson distribution. By comparing with previous work on other systems, we suggest that the statistical characteristics of the energy landscape may be a "universal" feature of all proteins.

The interpretation of site-directed mutagenesis experiments by linear free energy relations.

Fersht and co-workers have applied a linear free energy relation (Brønsted equation) to analyze site-directed mutagenesis experiments involving the enzyme tyrosyl-tRNA synthetase and have suggested that the Brønsted exponent is linearly correlated with the value of the reaction coordinate at the transition state. We point out that when the mutants differ solely through the formation or deletion of a hydrogen bond away from the reaction center, a linear free energy relation is expected only in limiting cases for which the Brønsted relation exponent is 0, 1 or infinity. The results may be correlated with a conformational coordinate but not with the development of the reaction coordinate per se.