Thermal folding and mechanical unfolding pathways of protein secondary structures

Unfolding mechanism and free energy landscape of single, stable, alpha helices at low pull speeds

Single alpha helices (SAHs) stable in isolated form are often found in motor proteins where they bridge functional domains. Understanding the mechanical response of SAHs is thus critical to understand their function. The quasi-static force-extension relation of a small number of SAHs is known from single-molecule experiments. Unknown, or still controversial, are the molecular scale details behind those observations. We show that the deformation mechanism of SAHs pulled from the termini at pull speeds approaching the quasi-static limit differs from that of typical helices found in proteins, which are stable only when interacting with other protein domains. Using molecular dynamics simulations with atomistic resolution at low pull speeds previously inaccessible to simulation, we show that SAHs start unfolding from the termini at all pull speeds we investigated. Unfolding proceeds residue-by-residue and hydrogen bond breaking is not the main event determining the barrier to unfolding. We use the molecular simulation data to test the cooperative sticky chain model. This model yields excellent fits of the force-extension curves and quantifies the distance, xE = 0.13 nm, to the transition state, the natural frequency of bond vibration, ν0 = 0.82 ns-1, and the height, V0 = 2.9 kcal mol-1, of the free energy barrier associated with the deformation of single residues. Our results demonstrate that the sticky chain model could advantageously be used to analyze experimental force-extension curves of SAHs and other biopolymers.

The relative stability of trpzip1 and its mutants determined by computation and experiment

The single-point mutations of tprzip1 are indicated at left, and their relative energetics are compared at right.

Mechanical unfolding of alpha- and beta-helical protein motifs

Alpha-helices and beta-sheets are the two most common secondary structure motifs in proteins. Beta-helical structures merge features of the two motifs, containing two or three beta-sheet faces connected by loops or turns in a single protein. Beta-helical structures form the basis of proteins with diverse mechanical functions such as bacterial adhesins, phage cell-puncture devices, antifreeze proteins, and extracellular matrices. Alpha-helices are commonly found in cellular and extracellular matrix components, whereas beta-helices such as curli fibrils are more common as bacterial and biofilm matrix components. It is currently not known whether it may be advantageous to use one helical motif over the other for different structural and mechanical functions. To better understand the mechanical implications of using different helix motifs in networks, here we use Steered Molecular Dynamics (SMD) simulations to mechanically unfold multiple alpha- and beta-helical proteins at constant velocity at the single molecule scale. We focus on the energy dissipated during unfolding as a means of comparison between proteins and work normalized by protein characteristics (initial and final length, # H-bonds, # residues, etc.). We find that although alpha-helices such as keratin and beta-helices CsgA and CsgB can require similar amounts of work to unfold, the normalized work per hydrogen bond, initial end to end length, and number of residues is greater for beta-helices at the same pulling rate. To explain this, we analyze the orientation of the backbone alpha carbons and backbone hydrogen bonds during unfolding. We find that the larger width and shorter height of beta-helices results in smaller angles between the protein backbone and the pulling direction during unfolding. As subsequent strands are separated from the beta-helix core, the angle between the backbone and the pulling direction diminishes. This marks a transition where beta-sheet hydrogen bonds become loaded predominantly in a collective shearing mode, which requires a larger rupture force. This finding underlines the importance of geometry in optimizing resistance to mechanical unfolding in proteins. The helix radius is identified here as an important parameter that governs how much sacrificial energy dissipation capacity can be stored in protein networks, where beta-helices offer unique properties.

Binding of interferon reduces the force of unfolding for interferon receptor 1

Differential signaling of the type I interferon receptor (IFNAR) has been correlated with the ability of its subunit, IFNAR1, to differentially recognize a large spectrum of different ligands, which involves intricate conformational re-arrangements of multiple interacting domains. To shed light onto the structural determinants governing ligand recognition, we compared the force-induced unfolding of the IFNAR1 ectodomain when bound to interferon and when free, using the atomic force microscope and steered molecular dynamics simulations. Unexpectedly, we find that IFNAR1 is easier to mechanically unfold when bound to interferon than when free. Analysis of the structures indicated that the origin of the reduction in unfolding forces is a conformational change in IFNAR1 induced by ligand binding.

As Simple As Possible, but Not Simpler: Exploring the Fidelity of Coarse-Grained Protein Models for Simulated Force Spectroscopy

Mechanical unfolding of a single domain of loop-truncated superoxide dismutase protein has been simulated via force spectroscopy techniques with both all-atom (AA) models and several coarse-grained models having different levels of resolution: A Gō model containing all heavy atoms in the protein (HA-Gō), the associative memory, water mediated, structure and energy model (AWSEM) which has 3 interaction sites per amino acid, and a Gō model containing only one interaction site per amino acid at the Cα position (Cα-Gō). To systematically compare results across models, the scales of time, energy, and force had to be suitably renormalized in each model. Surprisingly, the HA-Gō model gives the softest protein, exhibiting much smaller force peaks than all other models after the above renormalization. Clustering to render a structural taxonomy as the protein unfolds showed that the AA, HA-Gō, and Cα-Gō models exhibit a single pathway for early unfolding, which eventually bifurcates repeatedly to multiple branches only after the protein is about half-unfolded. The AWSEM model shows a single dominant unfolding pathway over the whole range of unfolding, in contrast to all other models. TM alignment, clustering analysis, and native contact maps show that the AWSEM pathway has however the most structural similarity to the AA model at high nativeness, but the least structural similarity to the AA model at low nativeness. In comparison to the AA model, the sequence of native contact breakage is best predicted by the HA-Gō model. All models consistently predict a similar unfolding mechanism for early force-induced unfolding events, but diverge in their predictions for late stage unfolding events when the protein is more significantly disordered.

Steered molecular dynamics analysis of the role of cofilin in increasing the flexibility of actin filaments

Cofilin is one of the most essential regulatory proteins and participates in the process of disassembling actin filaments. Cofilin induces conformational changes to actin filaments, and both the bending and torsional rigidity of the filament. In this study, we investigate the effects of cofilin on the mechanical properties of actin filaments using computational methods. Three models defined by their number of bound cofilins are constructed using coarse-grained MARTINI force field, and they are then extended with steered molecular dynamics simulation. After obtaining the stress-strain curves of the models, we calculate their Young's moduli and other mechanical properties that have not yet been determined for actin filaments. We analyze the cause of the different behaviors of the three models based on their atomistic geometrical differences. Finally, it is demonstrated that cofilin binding causes changes in the distances, angles, and stabilities of the residues in actin filaments.

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.

Simulation of the opening and closing of Hsp70 chaperones by coarse-grained molecular dynamics

Heat-shock proteins 70 (Hsp70s) are key molecular chaperones which assist in the folding and refolding/disaggregation of proteins. Hsp70s, which consist of a nucleotide-binding domain (NBD, consisting of NBD-I and NBD-II subdomains) and a substrate-binding domain [SBD, further split into the β-sheet (SBD-β) and α-helical (SBD-α) subdomains], occur in two major conformations having (a) a closed SBD, in which the SBD and NBD domains do not interact, (b) an open SBD, in which SBD-α interacts with NBD-I and SBD-β interacts with the top parts of NBD-I and NBD-II. In the SBD-closed conformation, SBD is bound to a substrate protein, with release occurring after transition to the open conformation. While the transition from the closed to the open conformation is triggered efficiently by binding of adenosine triphosphate (ATP) to the NBD, it also occurs, although less frequently, in the absence of ATP. The reverse transition occurs after ATP hydrolysis. Here, we report canonical and multiplexed replica exchange simulations of the conformational dynamics of Hsp70s using a coarse-grained molecular dynamics approach with the UNRES force field. The simulations were run in the following three modes: (i) with the two halves of the NBD unrestrained relative to each other, (ii) with the two halves of the NBD restrained in an "open" geometry as in the SBD-closed form of DnaK (2KHO), and (iii) the two halves of NBD restrained in a "closed" geometry as in known experimental structures of ATP-bound NBD forms of Hsp70. Open conformations, in which the SBD interacted strongly with the NBD, formed spontaneously during all simulations; the number of transitions was largest in simulations carried out with the "closed" NBD domain, and smallest in those carried out with the "open" NBD domain; this observation is in agreement with the experimentally-observed influence of ATP-binding on the transition of Hsp70's from the SBD-closed to the SBD-open form. Two kinds of open conformations were observed: one in which SBD-α interacts with NBD-I and SBD-β interacts with the top parts of NBD-I and NBD-II (as observed in the structures of nucleotide exchange factors), and another one in which this interaction pattern is swapped. A third type of motion, in which SBD-α binds to NBD without dissociating from SBD-β was also observed. It was found that the first stage of interdomain communication (approach of SBD-β, to NBD) is coupled with the rotation of the long axes of NBD-I and NBD-II towards each other. To the best of our knowledge, this is the first successful simulation of the full transition of an Hsp70 from the SBD-closed to the SBD-open conformation.

Protein high-force pulling simulations yield low-force results

All-atom explicit-solvent molecular dynamics simulations are used to pull with extremely large constant force (750–3000 pN) on three small proteins. The introduction of a nondimensional timescale permits direct comparison of unfolding across all forces. A crossover force of approximately 1100 pN divides unfolding dynamics into two regimes. At higher forces, residues sequentially unfold from the pulling end while maintaining the remainder of the protein force-free. Measurements of hydrodynamic viscous stresses are made easy by the high speeds of unfolding. Using an exact low-Reynolds-number scaling, these measurements can be extrapolated to provide, for the first time, an estimate of the hydrodynamic force on low-force unfolding. Below 1100 pN, but surprisingly still at extremely large applied force, intermediate states and cooperative unfoldings as seen at much lower forces are observed. The force-insensitive persistence of these structures indicates that decomposition into unfolded fragments requires a large fluctuation. This finding suggests how proteins are constructed to resist transient high force. The progression of helix and sheet unfolding is also found to be insensitive to force. The force-insensitivity of key aspects of unfolding opens the possibility that numerical simulations can be accelerated by high applied force while still maintaining critical features of unfolding.

Determination of effective potentials for the stretching of C(α) ⋯ C(α) virtual bonds in polypeptide chains for coarse-grained simulations of proteins from ab initio energy surfaces of N-methylacetamide and N-acetylpyrrolidine

The potentials of mean force (PMF's) for the deformation of the C(α) ⋯ C(α) virtual bonds in polypeptide chains were determined from the diabatic energy surfaces of N-methylacetamide (modeling regular peptide groups) and N-acetylpyrrolidine (modeling the peptide groups preceding proline), calculated at the Møller-Plesset (MP2) ab initio level of theory with the 6-31G(d,p) basis set. The energy surfaces were expressed in the C(α) ⋯ C(α) virtual-bond length (d) and the H-N-C(α) ⋯ C' improper dihedral angle (α) that describes the pyramidicity of the amide nitrogen, or in the C(α)-C'(O)-N-C(α) dihedral angle (ω) and the angle α. For each grid point, the potential energy was minimized with respect to all remaining degrees of freedom. The PMF's obtained from the (d, α) energy surfaces produced realistic free-energy barriers to the trans-cis transition (10 kcal/mol and 13 kcal/mol for the regular and proline peptide groups, respectively, compared to 12.6 - 13.9 kcal/mol and 17.3 - 19.6 kcal/mol determined experimentally for glycylglycine and N-acylprolines, respectively), while those obtained from the (ω, α) energy maps produced either low-quality PMF curves when direct Boltzmann summation was implemented to compute the PMF's or too-flat curves with too-low free-energy barriers to the trans-cis transition if harmonic extrapolation was used to estimate the contributions to the partition function. An analytical bimodal logarithmic-Gaussian expression was fitted to the PMF's, and the potentials were implemented in the UNRES force field. Test Langevin-dynamics simulations were carried out for the Gly-Gly and Gly-Pro dipeptides, which showed a 10(6)-fold increase of the simulated rate of the trans-cis isomerization with respect to that measured experimentally; effectively the same result was obtained with the analytical Kramers theory of reaction rate applied to the UNRES representation of the peptide groups. Application of Kramers' theory to compute the rate constants from the all-atom ab initio energy surfaces of the model compounds studied resulted in isomerization rates close to the experimental values, which demonstrates that the increase of the isomerization rate in UNRES simulations results solely from averaging out the secondary degrees of freedom.

Implementation of molecular dynamics and its extensions with the coarse-grained UNRES force field on massively parallel systems; towards millisecond-scale simulations of protein structure, dynamics, and thermodynamics

We report the implementation of our united-residue UNRES force field for simulations of protein structure and dynamics with massively parallel architectures. In addition to coarse-grained parallelism already implemented in our previous work, in which each conformation was treated by a different task, we introduce a fine-grained level in which energy and gradient evaluation are split between several tasks. The Message Passing Interface (MPI) libraries have been utilized to construct the parallel code. The parallel performance of the code has been tested on a professional Beowulf cluster (Xeon Quad Core), a Cray XT3 supercomputer, and two IBM BlueGene/P supercomputers with canonical and replica-exchange molecular dynamics. With IBM BlueGene/P, about 50 % efficiency and 120-fold speed-up of the fine-grained part was achieved for a single trajectory of a 767-residue protein with use of 256 processors/trajectory. Because of averaging over the fast degrees of freedom, UNRES provides an effective 1000-fold speed-up compared to the experimental time scale and, therefore, enables us to effectively carry out millisecond-scale simulations of proteins with 500 and more amino-acid residues in days of wall-clock time.

The folding of knotted proteins: insights from lattice simulations

We carry out systematic Monte Carlo simulations of Gō lattice proteins to investigate and compare the folding processes of two model proteins whose native structures differ from each other due to the presence of a trefoil knot located near the terminus of one of the protein chains. We show that the folding time of the knotted fold is larger than that of the unknotted protein and that this difference in folding time is particularly striking in the temperature region below the optimal folding temperature. Both proteins display similar folding transition temperatures, which is indicative of similar thermal stabilities. By using the folding probability reaction coordinate as an estimator of folding progression we have found out that the formation of the knot is mainly a late folding event in our shallow knot system.

Experimental and computational characterization of biological liquid crystals: a review of single-molecule bioassays

Quantitative understanding of the mechanical behavior of biological liquid crystals such as proteins is essential for gaining insight into their biological functions, since some proteins perform notable mechanical functions. Recently, single-molecule experiments have allowed not only the quantitative characterization of the mechanical behavior of proteins such as protein unfolding mechanics, but also the exploration of the free energy landscape for protein folding. In this work, we have reviewed the current state-of-art in single-molecule bioassays that enable quantitative studies on protein unfolding mechanics and/or various molecular interactions. Specifically, single-molecule pulling experiments based on atomic force microscopy (AFM) have been overviewed. In addition, the computational simulations on single-molecule pulling experiments have been reviewed. We have also reviewed the AFM cantilever-based bioassay that provides insight into various molecular interactions. Our review highlights the AFM-based single-molecule bioassay for quantitative characterization of biological liquid crystals such as proteins.

Large Protein Dynamics Described by Hierarchical-Component Mode Synthesis

Protein dynamics has played a pivotal role in understanding the biological function of protein. For investigation of such dynamics, normal-mode analysis (NMA) has been broadly employed with atomistic model and/or coarse-grained models such as elastic network model (ENM). For large protein complexes, NMA with even ENM encounters the expensive computational process such as diagonalization of Hessian (stiffness) matrix. Here, we suggest the hierarchical-component mode synthesis (hCMS), which allows the fast computation of low-frequency normal modes related to conformational change. Specifically, a large protein structure is regarded as a combination of several structural units, for which the eigen-value problem is utilized for obtaining the frequencies and their normal modes for each structural unit, and consequently, such frequencies and normal modes are assembled with geometrical constraint for interface between structural units in order to find the low-frequency normal modes of a large protein complex. It is shown that hCMS is able to provide the normal modes with accuracy, quantitatively comparable to those of original NMA. This implies that hCMS may enable the computationally efficient analysis of large protein dynamics.

Effects of molecular crowding on stretching of polymers in poor solvent

We consider a linear polymer chain in a disordered environment modeled by percolation clusters on a square lattice. The disordered environment is meant to roughly represent molecular crowding as seen in cells. The model may be viewed as the simplest representation of biopolymers in a cell. We show the existence of intermediate states during stretching arising as a consequence of molecular crowding. In the constant distance ensemble the force-extension curves exhibit oscillations. We observe the emergence of two or more peaks in the probability distribution curves signaling the coexistence of different states and indicating that the transition is discontinuous unlike what is observed in the absence of molecular crowding.

Application of Multiplexed Replica Exchange Molecular Dynamics to the UNRES Force Field: Tests with alpha and alpha+beta Proteins

The replica exchange (RE) method is increasingly used to improve sampling in molecular dynamics (MD) simulations of biomolecular systems. Recently, we implemented the united-residue UNRES force field for mesoscopic MD. Initial results from UNRES MD simulations show that we are able to simulate folding events that take place in a microsecond or even a millisecond time scale. To speed up the search further, we applied the multiplexing replica exchange molecular dynamics (MREMD) method. The multiplexed variant (MREMD) of the RE method, developed by Rhee and Pande, differs from the original RE method in that several trajectories are run at a given temperature. Each set of trajectories run at a different temperature constitutes a layer. Exchanges are attempted not only within a single layer but also between layers. The code has been parallelized and scales up to 4000 processors. We present a comparison of canonical MD, REMD, and MREMD simulations of protein folding with the UNRES force-field. We demonstrate that the multiplexed procedure increases the power of replica exchange MD considerably and convergence of the thermodynamic quantities is achieved much faster.

Hierarchies, multiple energy barriers, and robustness govern the fracture mechanics of alpha-helical and beta-sheet protein domains

The fundamental fracture mechanisms of biological protein materials remain largely unknown, in part, because of a lack of understanding of how individual protein building blocks respond to mechanical load. For instance, it remains controversial whether the free energy landscape of the unfolding behavior of proteins consists of multiple, discrete transition states or the location of the transition state changes continuously with the pulling velocity. This lack in understanding has thus far prevented us from developing predictive strength models of protein materials. Here, we report direct atomistic simulation that over four orders of magnitude in time scales of the unfolding behavior of alpha-helical (AH) and beta-sheet (BS) domains, the key building blocks of hair, hoof, and wool as well as spider silk, amyloids, and titin. We find that two discrete transition states corresponding to two fracture mechanisms exist. Whereas the unfolding mechanism at fast pulling rates is sequential rupture of individual hydrogen bonds (HBs), unfolding at slow pulling rates proceeds by simultaneous rupture of several HBs. We derive the hierarchical Bell model, a theory that explicitly considers the hierarchical architecture of proteins, providing a rigorous structure-property relationship. We exemplify our model in a study of AHs, and show that 3-4 parallel HBs per turn are favorable in light of the protein's mechanical and thermodynamical stability, in agreement with experimental findings that AHs feature 3.6 HBs per turn. Our results provide evidence that the molecular structure of AHs maximizes its robustness at minimal use of building materials.

Mechanical unfolding revisited through a simple but realistic model

Single-molecule experiments and their application to probe the mechanical resistance and related properties of proteins provide a new dimension in our knowledge of these important and complex biological molecules. Single-molecule techniques may not have yet overridden solution experiments as a method of choice to characterize biophysical and biological properties of proteins, but have stimulated a debate and contributed considerably to bridge theory and experiment. Here we demonstrate this latter contribution by illustrating the reach of some theoretical findings using a solvable but nontrivial molecular model whose properties are analogous to those of the corresponding experimental systems. In particular, we show the relationship between the thermodynamic and the mechanical properties of a protein. The simulations presented here also illustrate how forced and spontaneous unfolding occur through different pathways and that folding and unfolding rates at equilibrium cannot in general be obtained from forced unfolding experiments or simulations. We also study the relationship between the energy surface and the mechanical resistance of a protein and show how a simple analysis of the native state can predict much of the mechanical properties of a protein.

Molecular dynamics with the United-residue force field: ab initio folding simulations of multichain proteins

The implementation of molecular dynamics with the united-residue (UNRES) force field is extended to treat multichain proteins. Constant temperature was maintained in the simulations with Berendsen or Langevin thermostats. The method was tested on three alpha-helical proteins (1G6U and GCN4-p1, each with two chains, and 1C94, with four chains). Simulations were carried out for both the isolated single chains and the multichain complexes. The proteins were folded by starting from the extended conformation with random initial velocities and with the chains parallel to each other. No symmetry constraints or structure information were included for the single chains or the multichain complexes. In the case of single-chain simulations, a high percentage of the trajectories (100% for 1G6U, 90% for GCN4-p1, and 80% for 1C94) converged to nativelike structures (assumed as the experimental structure of a monomer in the multichain complex), showing that, for the proteins studied in this work with the UNRES force field, the interactions between chains are not critical for stabilization of the individual chains. In the case of multichain simulations, the native structures of the 1G6U and GCN4-p1 complexes, but not that of 1C94, are predicted successfully. The association of the subunits does not follow a unique mechanism; the monomers were observed to fold both before and simultaneously with their association.

Protein-Folding Dynamics: Overview of Molecular Simulation Techniques

Molecular dynamics (MD) is an invaluable tool with which to study protein folding in silico. Although just a few years ago the dynamic behavior of a protein molecule could be simulated only in the neighborhood of the experimental conformation (or protein unfolding could be simulated at high temperature), the advent of distributed computing, new techniques such as replica-exchange MD, new approaches (based on, e.g., the stochastic difference equation), and physics-based reduced models of proteins now make it possible to study protein-folding pathways from completely unfolded structures. In this review, we present algorithms for MD and their extensions and applications to protein-folding studies, using all-atom models with explicit and implicit solvent as well as reduced models of polypeptide chains.

Modification and optimization of the united-residue (UNRES) potential energy function for canonical simulations. I. Temperature dependence of the effective energy function and tests of the optimization method with single training proteins

We report the modification and parametrization of the united-residue (UNRES) force field for energy-based protein structure prediction and protein folding simulations. We tested the approach on three training proteins separately: 1E0L (beta), 1GAB (alpha), and 1E0G (alpha + beta). Heretofore, the UNRES force field had been designed and parametrized to locate native-like structures of proteins as global minima of their effective potential energy surfaces, which largely neglected the conformational entropy because decoys composed of only lowest-energy conformations were used to optimize the force field. Recently, we developed a mesoscopic dynamics procedure for UNRES and applied it with success to simulate protein folding pathways. However, the force field turned out to be largely biased toward -helical structures in canonical simulations because the conformational entropy had been neglected in the parametrization. We applied the hierarchical optimization method, developed in our earlier work, to optimize the force field; in this method, the conformational space of a training protein is divided into levels, each corresponding to a certain degree of native-likeness. The levels are ordered according to increasing native-likeness; level 0 corresponds to structures with no native-like elements, and the highest level corresponds to the fully native-like structures. The aim of optimization is to achieve the order of the free energies of levels, decreasing as their native-likeness increases. The procedure is iterative, and decoys of the training protein(s) generated with the energy function parameters of the preceding iteration are used to optimize the force field in a current iteration. We applied the multiplexing replica-exchange molecular dynamics (MREMD) method, recently implemented in UNRES, to generate decoys; with this modification, conformational entropy is taken into account. Moreover, we optimized the free-energy gaps between levels at temperatures corresponding to a predominance of folded or unfolded structures, as well as to structures at the putative folding-transition temperature, changing the sign of the gaps at the transition temperature. This enabled us to obtain force fields characterized by a single peak in the heat capacity at the transition temperature. Furthermore, we introduced temperature dependence to the UNRES force field; this is consistent with the fact that it is a free-energy and not a potential energy function. beta

Molecular simulation of the reversible mechanical unfolding of proteins

In this work we have combined a Wang-Landau sampling scheme [F. Wang and D. Landau, Phys. Rev. Lett. 86, 2050 (2001)] with an expanded ensemble formalism to yield a simple and powerful method for computing potentials of mean force. The new method is implemented to investigate the mechanical deformation of proteins. Comparisons are made with analytical results for simple model systems such as harmonic springs and Rouse chains. The method is then illustrated on a model 15-residue alanine molecule in an implicit solvent. Results for mechanical unfolding of this oligopeptide are compared to those of steered molecular dynamics calculations.

Mechanical properties of the domains of titin in a Go-like model

Comparison of properties of three domains of titin, I1, I27, and I28, in a simple geometry-based model shows that despite a high structural homology between their native states different domains show similar but distinguishable mechanical properties. Folding properties of the separate domains are predicted to be diversified which reflects sensitivity of the kinetics to the details of native structures. The Go-like model corresponding to the experimentally resolved native structure of the I1 domain is found to provide the biggest thermodynamic and mechanical stability compared to the other domains studied here. We analyze elastic, thermodynamic, and kinetic properties of several structures corresponding to the I28 domain as obtained through homology-based modeling. We discuss the ability of the models of the I28 domain to reproduce experimental results qualitatively. A strengthening of contacts that involve hydrophobic amino acids does not affect theoretical comparisons of the domains. Tandem linkages of up to five identical or different domains unravel in a serial fashion at low temperatures. We study the nature of the intermediate state that arises in the early stages of the serial unraveling and find it to qualitatively agree with the results of Marszalek et al.

Kinetic Studies of Folding of the B-domain of Staphylococcal Protein A with Molecular Dynamics and a United-residue (UNRES) Model of Polypeptide Chains

Langevin dynamics is used with our physics-based united-residue (UNRES) force field to study the folding pathways of the B-domain of staphylococcal protein A (1BDD (alpha; 46 residues)). With 400 trajectories of protein A started from the extended state (to gather meaningful statistics), and simulated for more than 35 ns each, 380 of them folded to the native structure. The simulations were carried out at the optimal folding temperature of protein A with this force field. To the best of our knowledge, this is the first simulation study of protein-folding kinetics with a physics-based force field in which reliable statistics can be gathered. In all the simulations, the C-terminal alpha-helix forms first. The ensemble of the native basin has an average RMSD value of 4 A from the native structure. There is a stable intermediate along the folding pathway, in which the N-terminal alpha-helix is unfolded; this intermediate appears on the way to the native structure in less than one-fourth of the folding pathways, while the remaining ones proceed directly to the native state. Non-native structures persist until the end of the simulations, but the native-like structures dominate. To express the kinetics of protein A folding quantitatively, two observables were used: (i) the average alpha-helix content (averaged over all trajectories within a given time window); and (ii) the fraction of conformations (averaged over all trajectories within a given time window) with Calpha RMSD values from the native structure less than 5 A (fraction of completely folded structures). The alpha-helix content grows quickly with time, and its variation fits well to a single-exponential term, suggesting fast two-state kinetics. On the other hand, the fraction of folded structures changes more slowly with time and fits to a sum of two exponentials, in agreement with the appearance of the intermediate, found when analyzing the folding pathways. This observation demonstrates that different qualitative and quantitative conclusions about folding kinetics can be drawn depending on which observable is monitored.

Calculation of the entropy and free energy from monte carlo simulations of a peptide stretched by an external force

Hypothetical scanning Monte Carlo (HSMC) is a method for calculating the absolute entropy, S, and free energy, F, from a trajectory generated by any simulation technique. HSMC was applied initially to fluids (argon and water) and later to peptides and self-avoiding walks on a lattice. In this paper we make a step further and apply it to a model of decaglycine (at T = 300 K) in vacuum with constant bond lengths where external stretching forces are exerted at the end points; the changes in S and F are calculated as the forces are increased. The molecule is placed initially in a helical structure, which is changed to an extended structure after a short simulation time due to the exerted forces. This study has relevance to problems in polymers (e.g., rubber elasticity) and to the analysis of experiments where individual molecules are stretched by atomic force microscopy (AFM), for example. The results for S and F are accurate and are significantly better than those obtained by the quasi-harmonic approximation and the local states method. However, the molecule is quite stiff due to the strong bond angle potentials and the extensions are small even for relatively large forces. Correspondingly, as the force is increased the decrease in the entropy is relatively small while the potential energy is enhanced significantly. Still, differences, TDeltaS, for different forces are obtained with very good accuracy of approximately 0.2 kcal/mol.

Mechanical resistance of proteins explained using simple molecular models

Recent experiments have demonstrated that proteins unfold when two atoms are mechanically pulled apart, and that this process is different to when heated or when a chemical denaturant is added to the solution. Experiments have also shown that the response of proteins to external forces is very diverse, some of them being "hard," and others "soft." Mechanical resistance originates from the presence of barriers on the energy landscape; together, experiment and simulation have demonstrated that unfolding occurs through alternative pathways when different pairs of atoms undergo mechanical extension. Here we use simulation to probe the mechanical resistance of six structurally diverse proteins when pulled in different directions. For this, we use two very different models: a detailed, transferable one, and a coarse-grained, structure-based one. The coarse-grained model gives results that are surprisingly similar to the detailed one and qualitatively agree with experiment; i.e., the mechanical resistance of different proteins or of a single protein pulled in different directions can be predicted by simulation. The results demonstrate the importance of pulling direction relative to the local topology in determining mechanical stability, and rationalize the effect of the location of importation/degradation tags on the rates of mitochondrial import or protein degradation in vivo.

Dissecting the mechanical unfolding of ubiquitin

The unfolding behavior of ubiquitin under the influence of a stretching force recently was investigated experimentally by single-molecule constant-force methods. Many observed unfolding traces had a simple two-state character, whereas others showed clear evidence of intermediate states. Here, we use Monte Carlo simulations to investigate the force-induced unfolding of ubiquitin at the atomic level. In agreement with experimental data, we find that the unfolding process can occur either in a single step or through intermediate states. In addition to this randomness, we find that many quantities, such as the frequency of occurrence of intermediates, show a clear systematic dependence on the strength of the applied force. Despite this diversity, one common feature can be identified in the simulated unfolding events, which is the order in which the secondary-structure elements break. This order is the same in two- and three-state events and at the different forces studied. The observed order remains to be verified experimentally but appears physically reasonable.

Molecular dynamics with the united-residue model of polypeptide chains. I. Lagrange equations of motion and tests of numerical stability in the microcanonical mode

The Lagrange formalism was implemented to derive the equations of motion for the physics-based united-residue (UNRES) force field developed in our laboratory. The C(alpha)...C(alpha) and C(alpha)...SC (SC denoting a side-chain center) virtual-bond vectors were chosen as variables. The velocity Verlet algorithm was adopted to integrate the equations of motion. Tests on the unblocked Ala(10) polypeptide showed that the algorithm is stable in short periods of time up to the time step of 1.467 fs; however, even with the shorter time step of 0.489 fs, some drift of the total energy occurs because of momentary jumps of the acceleration. These jumps are caused by numerical instability of the forces arising from the U(rot) component of UNRES that describes the energetics of side-chain-rotameric states. Test runs on the Gly(10) sequence (in which U(rot) is not present) and on the Ala(10) sequence with U(rot) replaced by a simple numerically stable harmonic potential confirmed this observation; oscillations of the total energy were observed only up to the time step of 7.335 fs, and some drift in the total energy or instability of the trajectories started to appear in long-time (2 ns and longer) trajectories only for the time step of 9.78 fs. These results demonstrate that the present U(rot) components (which are statistical potentials derived from the Protein Data Bank) must be replaced with more numerically stable functions; this work is under way in our laboratory. For the purpose of our present work, a nonsymplectic variable-time-step algorithm was introduced to reduce the energy drift for regular polypeptide sequences. The algorithm scales down the time step at a given point of a trajectory if the maximum change of acceleration exceeds a selected cutoff value. With this algorithm, the total energy is reasonably conserved up to a time step of 2.445 fs, as tested on the unblocked Ala(10) polypeptide. We also tried a symplectic multiple-time-step reversible RESPA algorithm and achieved satisfactory energy conservation for time steps up to 7.335 fs. However, at present, it appears that the reversible RESPA algorithm is several times more expensive than the variable-time-step algorithm because of the necessity to perform additional matrix multiplications. We also observed that, because Ala(10) folds and unfolds within picoseconds in the microcanonical mode, this suggests that the effective (event-based) time unit in UNRES dynamics is much larger than that of all-atom dynamics because of averaging over the fast-moving degrees of freedom in deriving the UNRES potential.

Molecular dynamics with the united-residue model of polypeptide chains. II. Langevin and Berendsen-bath dynamics and tests on model alpha-helical systems

The implementation of molecular dynamics (MD) with our physics-based protein united-residue (UNRES) force field, described in the accompanying paper, was extended to Langevin dynamics. The equations of motion are integrated by using a simplified stochastic velocity Verlet algorithm. To compare the results to those with all-atom simulations with implicit solvent in which no explicit stochastic and friction forces are present, we alternatively introduced the Berendsen thermostat. Test simulations on the Ala(10) polypeptide demonstrated that the average kinetic energy is stable with about a 5 fs time step. To determine the correspondence between the UNRES time step and the time step of all-atom molecular dynamics, all-atom simulations with the AMBER 99 force field and explicit solvent and also with implicit solvent taken into account within the framework of the generalized Born/surface area (GBSA) model were carried out on the unblocked Ala(10) polypeptide. We found that the UNRES time scale is 4 times longer than that of all-atom MD simulations because the degrees of freedom corresponding to the fastest motions in UNRES are averaged out. When the reduction of the computational cost for evaluation of the UNRES energy function is also taken into account, UNRES (with hydration included implicitly in the side chain-side chain interaction potential) offers about at least a 4000-fold speed up of computations relative to all-atom simulations with explicit solvent and at least a 65-fold speed up relative to all-atom simulations with implicit solvent. To carry out an initial full-blown test of the UNRES/MD approach, we ran Berendsen-bath and Langevin dynamics simulations of the 46-residue B-domain of staphylococcal protein A. We were able to determine the folding temperature at which all trajectories converged to nativelike structures with both approaches. For comparison, we carried out ab initio folding simulations of this protein at the AMBER 99/GBSA level. The average CPU time for folding protein A by UNRES molecular dynamics was 30 min with a single Alpha processor, compared to about 152 h for all-atom simulations with implicit solvent. It can be concluded that the UNRES/MD approach will enable us to carry out microsecond and, possibly, millisecond simulations of protein folding and, consequently, of the folding process of proteins in real time.

Theoretical studies of the kinetics of mechanical unfolding of cross-linked polymer chains and their implications for single-molecule pulling experiments

We have used kinetic Monte Carlo simulations to study the kinetics of unfolding of cross-linked polymer chains under mechanical loading. As the ends of a chain are pulled apart, the force transmitted by each cross-link increases until it ruptures. The stochastic cross-link rupture process is assumed to be governed by first order kinetics with a rate that depends exponentially on the transmitted force. We have performed random searches to identify optimal cross-link configurations whose unfolding requires a large applied force (measure of strength) and/or large dissipated energy (measure of toughness). We found that such optimal chains always involve cross-links arranged to form parallel strands. The location of those optimal strands generally depends on the loading rate. Optimal chains with a small number of cross-links were found to be almost as strong and tough as optimal chains with a large number of cross-links. Furthermore, optimality of chains with a small number of cross-links can be easily destroyed by adding cross-links at random. The present findings are relevant for the interpretation of single molecule force probe spectroscopy studies of the mechanical unfolding of "load-bearing" proteins, whose native topology often involves parallel strand arrangements similar to the optimal configurations identified in the study.

Ab initio simulations of protein-folding pathways by molecular dynamics with the united-residue model of polypeptide chains

We report the application of Langevin dynamics to the physics-based united-residue (UNRES) force field developed in our laboratory. Ten trajectories were run on seven proteins [PDB ID codes 1BDD (alpha; 46 residues), 1GAB (alpha; 47 residues), 1LQ7 (alpha; 67 residues), 1CLB (alpha; 75 residues), 1E0L (beta; 28 residues), and 1E0G (alpha+beta; 48 residues), and 1IGD (alpha+beta; 61 residues)] with the UNRES force field parameterized by using our recently developed method for obtaining a hierarchical structure of the energy landscape. All alpha-helical proteins and 1E0G folded to the native-like structures, whereas 1IGD and 1E0L yielded mostly nonnative alpha-helical folds although the native-like structures are lowest in energy for these two proteins, which can be attributed to neglecting the entropy factor in the current parameterization of UNRES. Average folding times for successful folding simulations were of the order of nanoseconds, whereas even the ultrafast-folding proteins fold only in microseconds, which implies that the UNRES time scale is approximately three orders of magnitude larger than the experimental time scale because the fast motions of the secondary degrees of freedom are averaged out. Folding with Langevin dynamics required 2-10 h of CPU time on average with a single AMD Athlon MP 2800+ processor depending on the size of the protein. With the advantage of parallel processing, this process leads to the possibility to explore thousands of folding pathways and to predict not only the native structure but also the folding scenario of a protein together with its quantitative kinetic and thermodynamic characteristics.

Stepwise unfolding of collapsed polymers

Motivated by recent experimental data on DNA stretching in presence of polyvalent counterions, we study the force-induced unfolding of a homopolymer on and off lattice. In the fixed force ensemble the globule unravels via a series of steps due to surface effects which play an important role for finite-size chains. This holds both for flexible and stiff polymers. We discuss in a qualitative way how this result may impact on the interpretation of DNA stretching experiments showing peaks in the characteristic curves, by extracting from the raw data the corresponding elongation- versus-force characteristic curves. Furthermore, approximate analytical and numerical calculations, valid in a quasi-equilibrium fixed stretch ensemble, and if the initial low-temperature state is ordered in a spool, show that the average force versus elongation displays peaks related to the geometry of the initial configuration. We finally argue how the proposed mechanisms identified for the arising of peaks may couple in the experiments, and comment on the role of dynamic effects.

Stretching of homopolymers and contact order

Mechanical stretching of self-interacting homopolymers is studied through molecular dynamics simulations in which the polymers are pulled with constant speed. At temperatures below the compactification temperature, the force-extension curves show a plateau that corresponds to the situation in which the polymer adopts "ball-string" configurations. The dependence of rupture distances on contact order and the effects of temperature are similar to those found in the case of model proteins. The dependence of behavior on the pulling speed is logarithmic. In the entropic limit, above the compactification temperature, the rupturing of contacts shows a monotonic decrease of extension with the contact order. The attainment of this limit depends on the system size and the pulling speed.

Thermal effects in stretching of Go-like models of titin and secondary structures

The effect of temperature on mechanical unfolding of proteins is studied using a Go-like model with a realistic contact map and Lennard-Jones contact interactions. The behavior of the I27 domain of titin and its serial repeats is contrasted to that of simple secondary structures. In all cases, thermal fluctuations accelerate the unraveling process, decreasing the unfolding force nearly linearly at low temperatures. However, differences in bonding geometry lead to different sensitivity to temperature and different changes in the unfolding pattern. Due to its special native-state geometry, titin is much more thermally and elastically stable than the secondary structures. At low temperatures, serial repeats of titin show a parallel unfolding of all domains to an intermediate state, followed by serial unfolding of the domains. At high temperatures, all domains unfold simultaneously, and the unfolding distance decreases monotonically with the contact order, that is, the sequence distance between the amino acids that form the native contact.

Stretching of proteins in the entropic limit

Mechanical stretching of six proteins is studied through molecular dynamics simulations. The model is Go-like, with Lennard-Jones interactions at native contacts. Low-temperature unfolding scenarios are remarkably complex and sensitive to small structural changes. Thermal fluctuations reduce the peak forces and the number of metastable states during unfolding. The unfolding pathways also simplify as temperature rises. In the entropic limit, all proteins show a monotonic decrease of the extension where bonds rupture with their separation along the backbone (contact order).

Stretching of a polymer below the theta point

The unfolding of a polymer below the theta point when pulled by an external force is studied both in d=2 on the lattice and in d=3 off the lattice. At T=0 and for finite length chains, it is found that the globule unfolds via multiple steps, corresponding to transitions between different minima, in both cases. In d=3 one of these intermediates is a regular helix. In the infinite length limit, these steps have a qualitative effect only in d=2. The phase diagram in d=2 is determined via the transfer matrix. To rationalize these results, energy-entropy and renormalization group arguments are given.

Protein folding: bringing theory and experiment closer together

The ability to perform enzymatic function by harnessing random molecular motion into self-organized protein structures is one of the most fascinating results of evolution. A close interplay between theory and experiment is driving the progress in understanding the principles that determine the behaviour of proteins. New techniques that significantly increase the amount of information obtainable from experimental data have been recently proposed; it is now becoming possible to describe at atomic resolution the events that take place during the folding process. Successful predictions of these events are being reported at an increasing rate and general principles are being outlined.

Universality classes in folding times of proteins

Molecular dynamics simulations in simplified models allow one to study the scaling properties of folding times for many proteins together under a controlled setting. We consider three variants of the Go models with different contact potentials and demonstrate scaling described by power laws and no correlation with the relative contact order parameter. We demonstrate existence of at least three kinetic universality classes that are correlated with the types of structure: the alpha-, alpha-beta-, and beta- proteins have the scaling exponents of approximately 1.7, 2.5, and 3.2, respectively. The three classes merge into one when the contact range is truncated at a reasonable value. We elucidate the role of the potential associated with the chirality of a protein.