Einstein Model of a Graph to Characterize Protein Folded/Unfolded States.

The folded structures of proteins can be accurately predicted by deep learning algorithms from their amino-acid sequences. By contrast, in spite of decades of research studies, the prediction of folding pathways and the unfolded and misfolded states of proteins, which are intimately related to diseases, remains challenging. A two-state (folded/unfolded) description of protein folding dynamics hides the complexity of the unfolded and misfolded microstates. Here, we focus on the development of simplified order parameters to decipher the complexity of disordered protein structures. First, we show that any connected, undirected, and simple graph can be associated with a linear chain of atoms in thermal equilibrium. This analogy provides an interpretation of the usual topological descriptors of a graph, namely the Kirchhoff index and Randic resistance, in terms of effective force constants of a linear chain. We derive an exact relation between the Kirchhoff index and the average shortest path length for a linear graph and define the free energies of a graph using an Einstein model. Second, we represent the three-dimensional protein structures by connected, undirected, and simple graphs. As a proof of concept, we compute the topological descriptors and the graph free energies for an all-atom molecular dynamics trajectory of folding/unfolding events of the proteins Trp-cage and HP-36 and for the ensemble of experimental NMR models of Trp-cage. The present work shows that the local, nonlocal, and global force constants and free energies of a graph are promising tools to quantify unfolded/disordered protein states and folding/unfolding dynamics. In particular, they allow the detection of transient misfolded rigid states.

Lysosomal enzyme tripeptidyl peptidase 1 destabilizes fibrillar Aß by multiple endoproteolytic cleavages within the ß-sheet domain.

Accumulation of amyloid-beta (Aß), which is associated with Alzheimer's disease, can be caused by excess production or insufficient clearance. Because of its ß-sheet structure, fibrillar Aß is resistant to proteolysis, which would contribute to slow degradation of Aß plaques in vivo. Fibrillar Aß can be internalized by microglia, which are the scavenger cells of the brain, but the fibrils are degraded only slowly in microglial lysosomes. Cathepsin B is a lysosomal protease that has been shown to proteolyze fibrillar Aß. Tripeptidyl peptidase 1 (TPP1), a lysosomal serine protease, possesses endopeptidase activity and has been shown to cleave peptides between hydrophobic residues. Herein, we demonstrate that TPP1 is able to proteolyze fibrillar Aß efficiently. Mass spectrometry analysis of peptides released from fibrillar Aß digested with TPP1 reveals several endoproteolytic cleavages including some within ß-sheet regions that are important for fibril formation. Using molecular dynamics simulations, we demonstrate that these cleavages destabilize fibrillar ß-sheet structure. The demonstration that TPP1 can degrade fibrillar forms of Aß provides insight into the turnover of fibrillar Aß and may lead to new therapeutic methods to increase degradation of Aß plaques.

Sequence-, structure-, and dynamics-based comparisons of structurally homologous CheY-like proteins.

We recently introduced a physically based approach to sequence comparison, the property factor method (PFM). In the present work, we apply the PFM approach to the study of a challenging set of sequences-the bacterial chemotaxis protein CheY, the N-terminal receiver domain of the nitrogen regulation protein NT-NtrC, and the sporulation response regulator Spo0F. These are all response regulators involved in signal transduction. Despite functional similarity and structural homology, they exhibit low sequence identity. PFM sequence comparison demonstrates a statistically significant qualitative difference between the sequence of CheY and those of the other two proteins that is not found using conventional alignment methods. This difference is shown to be consonant with structural characteristics, using distance matrix comparisons. We also demonstrate that residues participating strongly in native contacts during unfolding are distributed differently in CheY than in the other two proteins. The PFM result is also in accord with dynamic simulation results of several types. Molecular dynamics simulations of all three proteins were carried out at several temperatures, and it is shown that the dynamics of CheY are predicted to differ from those of NT-NtrC and Spo0F. The predicted dynamic properties of the three proteins are in good agreement with experimentally determined B factors and with fluctuations predicted by the Gaussian network model. We pinpoint the differences between the PFM and traditional sequence comparisons and discuss the informatic basis for the ability of the PFM approach to detect physical differences between these sequences that are not apparent from traditional alignment-based comparison.

Preventing fibril formation of a protein by selective mutation.

The origins of formation of an intermediate state involved in amyloid formation and ways to prevent it are illustrated with the example of the Formin binding protein 28 (FBP28) WW domain, which folds with biphasic kinetics. Molecular dynamics of protein folding trajectories are used to examine local and global motions and the time dependence of formation of contacts between C(α)s and C(ß)s of selected pairs of residues. Focus is placed on the WT FBP28 WW domain and its six mutants (L26D, L26E, L26W, E27Y, T29D, and T29Y), which have structures that are determined by high-resolution NMR spectroscopy. The origins of formation of an intermediate state are elucidated, viz. as formation of hairpin 1 by a hydrophobic collapse mechanism causing significant delay of formation of both hairpins, especially hairpin 2, which facilitates the emergence of an intermediate state. It seems that three-state folding is a major folding scenario for all six mutants and WT. Additionally, two-state and downhill folding scenarios were identified in ∼ 15% of the folding trajectories for L26D and L26W, in which both hairpins are formed by the Matheson-Scheraga mechanism much faster than in three-state folding. These results indicate that formation of hairpins connecting two antiparallel ß-strands determines overall folding. The correlations between the local and global motions identified for all folding trajectories lead to the identification of the residues making the main contributions in the formation of the intermediate state. The presented findings may provide an understanding of protein folding intermediates in general and lead to a procedure for their prevention.

Accounting for a mirror-image conformation as a subtle effect in protein folding.

By using local (free-energy profiles along the amino acid sequence and (13)C(α) chemical shifts) and global (principal component) analyses to examine the molecular dynamics of protein-folding trajectories, generated with the coarse-grained united-residue force field, for the B domain of staphylococcal protein A, we are able to (i) provide the main reason for formation of the mirror-image conformation of this protein, namely, a slow formation of the second loop and part of the third helix (Asp29-Asn35), caused by the presence of multiple local conformational states in this portion of the protein; (ii) show that formation of the mirror-image topology is a subtle effect resulting from local interactions; (iii) provide a mechanism for how protein A overcomes the barrier between the metastable mirror-image state and the native state; and (iv) offer a plausible reason to explain why protein A does not remain in the metastable mirror-image state even though the mirror-image and native conformations are at least energetically compatible.

Folding kinetics of WW domains with the united residue force field for bridging microscopic motions and experimental measurements.

To demonstrate the utility of the coarse-grained united-residue (UNRES) force field to compare experimental and computed kinetic data for folding proteins, we have performed long-time millisecond-timescale canonical Langevin molecular dynamics simulations of the triple ß-strand from the Formin binding protein 28 WW domain and six nonnatural variants, using UNRES. The results have been compared with available experimental data in both a qualitative and a quantitative manner. Complexities of the folding pathways, which cannot be determined experimentally, were revealed. The folding mechanisms obtained from the simulated folding kinetics are in agreement with experimental results, with a few discrepancies for which we have accounted. The origins of single- and double-exponential kinetics and their correlations with two- and three-state folding scenarios are shown to be related to the relative barrier heights between the various states. The rate constants obtained from time profiles of the fractions of the native, intermediate, and unfolded structures, and the kinetic equations fitted to them, correlate with the experimental values; however, they are about three orders of magnitude larger than the experimental ones for most of the systems. These differences are in agreement with the timescale extension derived by scaling down the friction of water and averaging out the fast degrees of freedom when passing from all-atom to a coarse-grained representation. Our results indicate that the UNRES force field can provide accurate predictions of folding kinetics of these WW domains, often used as models for the study of the mechanisms of proein folding.

Anomalous diffusion and dynamical correlation between the side chains and the main chain of proteins in their native state.

Structural fluctuations of a protein are essential for a protein to function and fold. By using molecular dynamics (MD) simulations of the model α/ß protein VA3 in its native state, the coupling between the main-chain (MC) motions [represented by coarse-grained dihedral angles (CGDAs) γ(n) based on four successive C(α) atoms (n - 1, n, n + 1, n + 2) along the amino acid sequence] and its side-chain (SC) motions [represented by CGDAs δ(n) formed by the virtual bond joining two consecutive C(α) atoms (n, n + 1) and the bonds joining these C(α) atoms to their respective C(ß) atoms] was analyzed. The motions of SCs (δ(n)) and MC (γ(n)) over time occur on similar free-energy profiles and were found to be subdiffusive. The fluctuations of the SCs (δ(n)) and those of the MC (γ(n)) are generally poorly correlated on a ps time-scale with a correlation increasing with time to reach a maximum value at about 10 ns. This maximum value is close to the correlation between the δ(n)(t) and γ(n)(t) time-series extracted from the entire duration of the MD runs (400 ns) and varies significantly along the amino acid sequence. High correlations between the SC and MC motions [δ(t) and γ(t) time-series] were found only in flexible regions of the protein for a few residues which contribute the most to the slowest collective modes of the molecule. These results are a possible indication of the role of the flexible regions of proteins for the biological function and folding.

Kinks, loops, and protein folding, with protein A as an example.

The dynamics and energetics of formation of loops in the 46-residue N-terminal fragment of the B-domain of staphylococcal protein A has been studied. Numerical simulations have been performed using coarse-grained molecular dynamics with the united-residue (UNRES) force field. The results have been analyzed in terms of a kink (heteroclinic standing wave solution) of a generalized discrete nonlinear Schrödinger (DNLS) equation. In the case of proteins, the DNLS equation arises from a C(α)-trace-based energy function. Three individual kink profiles were identified in the experimental three-α-helix structure of protein A, in the range of the Glu16-Asn29, Leu20-Asn29, and Gln33-Asn44 residues, respectively; these correspond to two loops in the native structure. UNRES simulations were started from the full right-handed α-helix to obtain a clear picture of kink formation, which would otherwise be blurred by helix formation. All three kinks emerged during coarse-grained simulations. It was found that the formation of each is accompanied by a local free energy increase; this is expressed as the change of UNRES energy which has the physical sense of the potential of mean force of a polypeptide chain. The increase is about 7 kcal/mol. This value can thus be considered as the free energy barrier to kink formation in full α-helical segments of polypeptide chains. During the simulations, the kinks emerge, disappear, propagate, and annihilate each other many times. It was found that the formation of a kink is initiated by an abrupt change in the orientation of a pair of consecutive side chains in the loop region. This resembles the formation of a Bloch wall along a spin chain, where the C(α) backbone corresponds to the chain, and the amino acid side chains are interpreted as the spin variables. This observation suggests that nearest-neighbor side chain-side chain interactions are responsible for initiation of loop formation. It was also found that the individual kinks are reflected as clear peaks in the principal modes of the analyzed trajectory of protein A, the shapes of which resemble the directional derivatives of the kinks along the chain. These observations suggest that the kinks of the DNLS equation determine the functionally important motions of proteins.

Nonexponential decay of internal rotational correlation functions of native proteins and self-similar structural fluctuations.

Structural fluctuations of a protein are essential for the function of native proteins and for protein folding. To understand how the main chain in the native state of a protein fluctuates on different time scales, we examined the rotational correlation functions (RCFs), C(t), of the backbone N-H bonds and of the dihedral angles γ built on four consecutive C(α) atoms. Using molecular dynamics simulations of a model α/ß protein (VA3) in its native state, we demonstrate that these RCFs decay as stretched exponentials, ln[C(t)] ≈ D(α)t(α) with a constant D(α) and an exponent α (0 < α < 0.35) varying with the free-energy profiles (FEPs) along the amino acid sequence. The probability distributions of the fluctuations of the main chain computed at short time scale (1 ps) were identical to those computed at large time scale (1 ns) if the time is rescaled by a factor depending on α < 1. This self-similar property and the nonexponential decays (α ≠ 1) of the RCFs are described by a rotational diffusion equation with a time-dependent diffusion coefficient D(t) = αD(α)t(α-1). The present findings agree with observations of subdiffusion (α < 1) of fluorescent probes within a protein molecule. The subdiffusion of (15)N-H bonds did not affect the value of the order parameter S(2) extracted from the NMR relaxation data by assuming normal diffusion (α = 1) of (15)N-H bonds on a nanosecond time scale. However, we found that the RCF does not converge to S(2) on the nanosecond time scale for residues with multiple-minima FEPs.

Mechanistic Kinetic Model Reveals How Amyloidogenic Hydrophobic Patches Facilitate the Amyloid-ß Fibril Elongation.

Abnormal aggregation of amyloid ß (Aß) peptides into fibrils plays a critical role in the development of Alzheimer's disease. A two-stage "dock-lock" model has been proposed for the Aß fibril elongation process. However, the mechanisms of the Aß monomer-fibril binding process have not been elucidated with the necessary molecular-level precision, so it remains unclear how the lock phase dynamics leads to the overall in-register binding of the Aß monomer onto the fibril. To gain mechanistic insights into this critical step during the fibril elongation process, we used molecular dynamics (MD) simulations with a physics-based coarse-grained UNited-RESidue (UNRES) force field and sampled extensively the dynamics of the lock phase process, in which a fibril-bound Aß(9-40) peptide rearranged to establish the native docking conformation. Analysis of the MD trajectories with Markov state models was used to quantify the kinetics of the rearrangement process and the most probable pathways leading to the overall native docking conformation of the incoming peptide. These revealed a key intermediate state in which an intra-monomer hairpin is formed between the central core amyloidogenic patch 18VFFA21 and the C-terminal hydrophobic patch 34LMVG37. This hairpin structure is highly favored as a transition state during the lock phase of the fibril elongation. We propose a molecular mechanism for facilitation of the Aß fibril elongation by amyloidogenic hydrophobic patches.

Probing Protein Aggregation Using the Coarse-Grained UNRES Force Field.

Protein aggregation is the cause of many, often lethal, diseases, including the Alzheimer's, Parkinson's, and Huntington's diseases, and familial amyloidosis. Theoretical investigation of the mechanism of this process, including the structures of the oligomeric intermediates which are the most toxic, is difficult because of long time scale of aggregation. Coarse-grained models, which enable us to extend the simulation time scale by three or more orders of magnitude, are, therefore, of great advantage in such studies. In this chapter, we describe the application of the physics-based UNited RESidue (UNRES) force field developed in our laboratory to study protein aggregation, in both free simulations and simulations of aggregation propagation from an existing template (seed), and illustrate it with the examples of Aß-peptide aggregation and Aß-peptide-assisted aggregation of the peptides derived from the repeat domains of tau (TauRD).

Exploring Structural Flexibility and Stability of α-Synuclein by the Landau-Ginzburg-Wilson Approach.

α-Synuclein (αS) is the principal protein component of the Lewy body and Lewy neurite deposits that are found in the brains of the victims of one of the most prevalent neurodegenerative disorders, Parkinson's disease. αS can be qualified as a chameleon protein because of the large number of different conformations that it is able to adopt: it is disordered under physiological conditions in solution, in equilibrium with a minor α-helical tetrameric form in the cytoplasm, and is α-helical when bound to a cell membrane. Also, in vitro, αS forms polymorphic amyloid fibrils with unique arrangements of cross-ß-sheet motifs. Therefore, it is of interest to elucidate the origins of the structural flexibility of αS and what makes αS stable in different conformations. We address these questions here by analyzing the experimental structures of the micelle-bound, tetrameric, and fibrillar αS in terms of a kink (heteroclinic standing wave solution) of a generalized discrete nonlinear Schrödinger equation. It is illustrated that without molecular dynamics simulations the kinks are capable of identifying the key residues causing structural flexibility of αS. Also, the stability of the experimental structures of αS is investigated by simulating heating/cooling trajectories using the Glauber algorithm. The findings are consistent with experiments.

Wild-Type α-Synuclein and Variants Occur in Different Disordered Dimers and Pre-Fibrillar Conformations in Early Stage of Aggregation.

α-Synuclein is a 140 amino-acid intrinsically disordered protein mainly found in the brain. Toxic α-synuclein aggregates are the molecular hallmarks of Parkinson's disease. In vitro studies showed that α-synuclein aggregates in oligomeric structures of several 10th of monomers and into cylindrical structures (fibrils), comprising hundred to thousands of proteins, with polymorphic cross-ß-sheet conformations. Oligomeric species, formed at the early stage of aggregation remain, however, poorly understood and are hypothezised to be the most toxic aggregates. Here, we studied the formation of wild-type (WT) and mutant (A30P, A53T, and E46K) dimers of α-synuclein using coarse-grained molecular dynamics. We identified two principal segments of the sequence with a higher propensity to aggregate in the early stage of dimerization: residues 36-55 and residues 66-95. The transient α-helices (residues 53-65 and 73-82) of α-synuclein monomers are destabilized by A53T and E46K mutations, which favors the formation of fibril native contacts in the N-terminal region, whereas the helix 53-65 prevents the propagation of fibril native contacts along the sequence for the WT in the early stages of dimerization. The present results indicate that dimers do not adopt the Greek key motif of the monomer fold in fibrils but form a majority of disordered aggregates and a minority (9-15%) of pre-fibrillar dimers both with intra-molecular and intermolecular ß-sheets. The percentage of residues in parallel ß-sheets is by increasing order monomer < disordered dimers < pre-fibrillar dimers. Native fibril contacts between the two monomers are present in the NAC domain for WT, A30P, and A53T and in the N-domain for A53T and E46K. Structural properties of pre-fibrillar dimers agree with rupture-force atomic force microscopy and single-molecule Förster resonance energy transfer available data. This suggests that the pre-fibrillar dimers might correspond to the smallest type B toxic oligomers. The probability density of the dimer gyration radius is multi-peaks with an average radius that is 10 Å larger than the one of the monomers for all proteins. The present results indicate that even the elementary α-synuclein aggregation step, the dimerization, is a complicated phenomenon that does not only involve the NAC region.

How main-chains of proteins explore the free-energy landscape in native states.

Understanding how a single native protein diffuses on its free-energy landscape is essential to understand protein kinetics and function. The dynamics of a protein is complex, with multiple relaxation times reflecting a hierarchical free-energy landscape. Using all-atom molecular dynamics simulations of an alpha/beta protein (crambin) and a beta-sheet polypeptide (BS2) in their "native" states, we demonstrate that the mean-square displacement of dihedral angles, defined by 4 successive C(alpha) atoms, increases as a power law of time, t(alpha), with an exponent alpha between 0.08 and 0.39 (alpha = 1 corresponds to Brownian diffusion), at 300 K. Residues with low exponents are located mainly in well-defined secondary elements and adopt 1 conformational substate. Residues with high exponents are found in loops/turns and chain ends and exist in multiple conformational substates, i.e., they move on multiple-minima free-energy profiles.

How Useful can the Voigt Profile be in Protein Folding Processes?

The analytical expression for the Voigt profile, along with its simplified forms for the Gaussian and Lorentzian dominance, is presented. The applicability of the Voigt profile in the description of anomalous diffusion phenomena, ubiquitous in different fields of science including protein folding, is discussed. It is shown that the Voigt profile is a good descriptor of the processes occurring in protein folding and in the native state. The usefulness of the Voigt profile in deriving important information of the diffusive motions in proteins from a quasielastic incoherent neutron scattering experiments is illustrated.

Missense Mutations Modify the Conformational Ensemble of the α-Synuclein Monomer Which Exhibits a Two-Phase Characteristic.

α-Synuclein is an intrinsically disordered protein occurring in different conformations and prone to aggregate in ß-sheet structures, which are the hallmark of the Parkinson disease. Missense mutations are associated with familial forms of this neuropathy. How these single amino-acid substitutions modify the conformations of wild-type α-synuclein is unclear. Here, using coarse-grained molecular dynamics simulations, we sampled the conformational space of the wild type and mutants (A30P, A53P, and E46K) of α-synuclein monomers for an effective time scale of 29.7 ms. To characterize the structures, we developed an algorithm, CUTABI (CUrvature and Torsion based of Alpha-helix and Beta-sheet Identification), to identify residues in the α-helix and ß-sheet from Cα -coordinates. CUTABI was built from the results of the analysis of 14,652 selected protein structures using the Dictionary of Secondary Structure of Proteins (DSSP) algorithm. DSSP results are reproduced with 93% of success for 10 times lower computational cost. A two-dimensional probability density map of α-synuclein as a function of the number of residues in the α-helix and ß-sheet is computed for wild-type and mutated proteins from molecular dynamics trajectories. The density of conformational states reveals a two-phase characteristic with a homogeneous phase (state B, ß-sheets) and a heterogeneous phase (state HB, mixture of α-helices and ß-sheets). The B state represents 40% of the conformations for the wild-type, A30P, and E46K and only 25% for A53T. The density of conformational states of the B state for A53T and A30P mutants differs from the wild-type one. In addition, the mutant A53T has a larger propensity to form helices than the others. These findings indicate that the equilibrium between the different conformations of the α-synuclein monomer is modified by the missense mutations in a subtle way. The α-helix and ß-sheet contents are promising order parameters for intrinsically disordered proteins, whereas other structural properties such as average gyration radius, R g , or probability distribution of R g cannot discriminate significantly the conformational ensembles of the wild type and mutants. When separated in states B and HB, the distributions of R g are more significantly different, indicating that global structural parameters alone are insufficient to characterize the conformational ensembles of the α-synuclein monomer.

Investigation of Phosphorylation-Induced Folding of an Intrinsically Disordered Protein by Coarse-Grained Molecular Dynamics.

Apart from being the most common mechanism of regulating protein function and transmitting signals throughout the cell, phosphorylation has an ability to induce disorder-to-order transition in an intrinsically disordered protein. In particular, it was shown that folding of the intrinsically disordered protein, eIF4E-binding protein isoform 2 (4E-BP2), can be induced by multisite phosphorylation. Here, the principles that govern the folding of phosphorylated 4E-BP2 (pT37pT46 4E-BP218-62) are investigated by analyzing canonical and replica exchange molecular dynamics trajectories, generated with the coarse-grained united-residue force field, in terms of local and global motions and the time dependence of formation of contacts between Cαs of selected pairs of residues. The key residues involved in the folding of the pT37pT46 4E-BP218-62 are elucidated by this analysis. The correlations between local and global motions are identified. Moreover, for a better understanding of the physics of the formation of the folded state, the experimental structure of the pT37pT46 4E-BP218-62 is analyzed in terms of a kink (heteroclinic standing wave solution) of a generalized discrete nonlinear Schrödinger equation. It is shown that without molecular dynamics simulations the kinks are able to identify not only the phosphorylated sites of protein, the key players in folding, but also the reasons for the weak stability of the pT37pT46 4E-BP218-62.

Evidence, from simulations, of a single state with residual native structure at the thermal denaturation midpoint of a small globular protein.

The folding of the B-domain of staphylococcal protein A has been studied by coarse-grained canonical and multiplexed replica-exchange molecular dynamics simulations with the UNRES force field in a broad range of temperatures (270 K < or = T < or = 350 K). In canonical simulations, the folding was found to occur either directly to the native state or through kinetic traps, mainly the topological mirror image of the native three-helix bundle. The latter folding scenario was observed more frequently at low temperatures. With increase of temperature, the frequency of the transitions between the folded and misfolded/unfolded states increased and the folded state became more diffuse with conformations exhibiting increased root-mean-square deviations from the experimental structure (from about 4 A at T = 300 K to 8.7 A at T = 325 K). An analysis of the equilibrium conformational ensemble determined from multiplexed replica exchange simulations at the folding-transition temperature (T(f) = 325 K) showed that the conformational ensemble at this temperature is a collection of conformations with residual secondary structures, which possess native or near-native clusters of nonpolar residues in place, and not a 50-50% mixture of fully folded and fully unfolded conformations. These findings contradict the quasi-chemical picture of two- or multistate protein folding, which assumes an equilibrium between the folded, unfolded, and intermediate states, with equilibrium shifting with temperature but with the native conformations remaining essentially unchanged. Our results also suggest that long-range hydrophobic contacts are the essential factor to keep the structure of a protein thermally stable.

Determination of side-chain-rotamer and side-chain and backbone virtual-bond-stretching potentials of mean force from AM1 energy surfaces of terminally-blocked amino-acid residues, for coarse-grained simulations of protein structure and folding. II. Results, comparison with statistical potentials, and implementation in the UNRES force field.

Using the harmonic-approximation approach of the accompanying article and AM1 energy surfaces of terminally blocked amino-acid residues, we determined physics-based side-chain rotamer potentials and the side-chain virtual-bond-deformation potentials of 19 natural amino-acid residues with side chains. The potentials were approximated by analytical formulas and implemented in the UNRES mesoscopic dynamics program. For comparison, the corresponding statistical potentials were determined from 19,682 high-resolution protein structures. The low free-energy region of both the AM1-derived and the statistical potentials is determined by the valence geometry and the L-chirality, and its size increases with side-chain flexibility and decreases with increasing virtual-bond-angle theta. The differences between the free energies of rotamers are greater for the AM1-derived potentials compared with the statistical potentials and, for alanine and other residues with small side chains, a region corresponding to the C(ax)(7) conformation has remarkably low free-energy for the AM1-derived potentials, as opposed to the statistical potentials. These differences probably result from the interactions between neighboring residues and indicate the need for introduction of cooperative terms accounting for the coupling between side-chain rotamer and backbone interactions. Both AM1-derived and statistical virtual-bond-deformation potentials are multimodal for flexible side chains and are topologically similar; however, the regions of minima of the statistical potentials are much narrower, which probably results from imposing restraints in structure determination. The force field with the new potentials was preliminarily optimized using the FBP WW domain (1E0L) and the engrailed homeodomain (1ENH) as training proteins and assessed to be reasonably transferable.

Investigation of protein folding by coarse-grained molecular dynamics with the UNRES force field.

Coarse-grained molecular dynamics simulations offer a dramatic extension of the time-scale of simulations compared to all-atom approaches. In this article, we describe the use of the physics-based united-residue (UNRES) force field, developed in our laboratory, in protein-structure simulations. We demonstrate that this force field offers about a 4000-times extension of the simulation time scale; this feature arises both from averaging out the fast-moving degrees of freedom and reduction of the cost of energy and force calculations compared to all-atom approaches with explicit solvent. With massively parallel computers, microsecond folding simulation times of proteins containing about 1000 residues can be obtained in days. A straightforward application of canonical UNRES/MD simulations, demonstrated with the example of the N-terminal part of the B-domain of staphylococcal protein A (PDB code: 1BDD, a three-alpha-helix bundle), discerns the folding mechanism and determines kinetic parameters by parallel simulations of several hundred or more trajectories. Use of generalized-ensemble techniques, of which the multiplexed replica exchange method proved to be the most effective, enables us to compute thermodynamics of folding and carry out fully physics-based prediction of protein structure, in which the predicted structure is determined as a mean over the most populated ensemble below the folding-transition temperature. By using principal component analysis of the UNRES folding trajectories of the formin-binding protein WW domain (PDB code: 1E0L; a three-stranded antiparallel beta-sheet) and 1BDD, we identified representative structures along the folding pathways and demonstrated that only a few (low-indexed) principal components can capture the main structural features of a protein-folding trajectory; the potentials of mean force calculated along these essential modes exhibit multiple minima, as opposed to those along the remaining modes that are unimodal. In addition, a comparison between the structures that are representative of the minima in the free-energy profile along the essential collective coordinates of protein folding (computed by principal component analysis) and the free-energy profile projected along the virtual-bond dihedral angles gamma of the backbone revealed the key residues involved in the transitions between the different basins of the folding free-energy profile, in agreement with existing experimental data for 1E0L .