An atomically detailed study of the folding pathways of protein A with the stochastic difference equation

The B domain of protein A retains residual structures in 6 M guanidinium chloride as revealed by hydrogen/deuterium-exchange NMR spectroscopy

The characterization of residual structures persistent in unfolded proteins is an important issue in studies of protein folding, because the residual structures present, if any, may form a folding initiation site and guide the subsequent folding reactions. Here, we studied the residual structures of the isolated B domain (BDPA) of staphylococcal protein A in 6 M guanidinium chloride. BDPA is a small three-helix-bundle protein, and until recently its folding/unfolding reaction has been treated as a simple two-state process between the native and the fully unfolded states. We employed a dimethylsulfoxide (DMSO)-quenched hydrogen/deuterium (H/D)-exchange 2D NMR techniques with the use of spin desalting columns, which allowed us to investigate the H/D-exchange behavior of individually identified peptide amide (NH) protons. We obtained H/D-exchange protection factors of the 21 NH protons that form an α-helical hydrogen bond in the native structure, and the majority of these NH protons were significantly protected with a protection factor of 2.0-5.2 in 6 M guanidinium chloride, strongly suggesting that these weakly protected NH protons form much stronger hydrogen bonds under native folding conditions. The results can be used to deduce the structure of an early folding intermediate, when such an intermediate is shown by other methods. Among three native helical regions, the third helix in the C-terminal side was highly protected and stabilized by side-chain salt bridges, probably acting as the folding initiation site of BDPA. The present results are discussed in relation to previous experimental and computational findings on the folding mechanisms of BDPA.

Expression and Characterization of Intein-Cyclized Trimer of Staphylococcus aureus Protein A Domain Z

Staphylococcus aureus protein A (SpA) is an IgG Fc-binding virulence factor that is widely used in antibody purification and as a scaffold to develop affinity molecules. A cyclized SpA Z domain could offer exopeptidase resistance, reduced chromatographic ligand leaching after single-site endopeptidase cleavage, and enhanced IgG binding properties by preorganization, potentially reducing conformational entropy loss upon binding. In this work, a Z domain trimer (Z3) was cyclized using protein intein splicing. Interactions of cyclic and linear Z3 with human IgG1 were characterized by differential scanning fluorimetry (DSF), surface plasmon resonance (SPR), and isothermal titration calorimetry (ITC). DSF showed a 5 ℃ increase in IgG1 melting temperature when bound by each Z3 variant. SPR showed the dissociation constants of linear and cyclized Z3 with IgG1 to be 2.9 nM and 3.3 nM, respectively. ITC gave association enthalpies for linear and cyclic Z3 with IgG1 of -33.0 kcal/mol and -32.7 kcal/mol, and -T∆S of association 21.2 kcal/mol and 21.6 kcal/mol, respectively. The compact cyclic Z3 protein contains 2 functional binding sites and exhibits carboxypeptidase Y-resistance. The results suggest cyclization as a potential approach toward more stable SpA-based affinity ligands, and this analysis may advance our understanding of protein engineering for ligand and drug development.

Integrating Rigidity Analysis into the Exploration of Protein Conformational Pathways Using RRT* and MC

To understand how proteins function on a cellular level, it is of paramount importance to understand their structures and dynamics, including the conformational changes they undergo to carry out their function. For the aforementioned reasons, the study of large conformational changes in proteins has been an interest to researchers for years. However, since some proteins experience rapid and transient conformational changes, it is hard to experimentally capture the intermediate structures. Additionally, computational brute force methods are computationally intractable, which makes it impossible to find these pathways which require a search in a high-dimensional, complex space. In our previous work, we implemented a hybrid algorithm that combines Monte-Carlo (MC) sampling and RRT*, a version of the Rapidly Exploring Random Trees (RRT) robotics-based method, to make the conformational exploration more accurate and efficient, and produce smooth conformational pathways. In this work, we integrated the rigidity analysis of proteins into our algorithm to guide the search to explore flexible regions. We demonstrate that rigidity analysis dramatically reduces the run time and accelerates convergence.

Investigation of alkaline effects on Protein A affinity ligands and resins using high resolution mass spectrometry

In this study, the enhanced alkaline stability of Protein A ligands and resins designed by protein engineering approaches is demonstrated. High throughput PreDictor™ plates were used to evaluate and compare the human Immunoglobulin G (IgG) static binding capacities (SBC) of MabSelect SuRe™ and MabSelect™ PrismA affinity chromatography (AC) resins after continuous incubation in 0.1-2.0 M NaOH for 1-72 h. The alkaline effect on the Protein A affinity ligand was studied by high resolution mass spectrometry (MS). The IgG binding capacity of the investigated AC resins show expected declining trends with increasing NaOH concentrations and incubation times. The decrease is larger for MabSelect SuRe than for MabSelect PrismA and occur at lower NaOH concentrations. MabSelect SuRe display high remaining binding capacity even after 72 h continuous incubation in 0.1 M NaOH, while higher concentrations induce an accentuated decline with incubation time. The MabSelect PrismA resin shows almost no effect on the binding capacity even after 72 h incubation in 0.5 M NaOH. Decline in capacity is first observed after 48 h incubation in 1.0 M NaOH, thus displaying the extreme alkaline stability of the PrismA affinity ligand. The MS analysis of the ligands, including a Protein A single B-domain, SuRe-domain and PrismA-domain clearly illustrate the increasing alkaline stability (B-domain < SuRe < PrismA) as the ligand has been refined using a protein engineering approach. Deamidation and ligand degradation could be monitored in relation to NaOH incubation conditions. Enzymatic digestion of MabSelect SuRe and MabSelect PrismA resins after alkaline incubation and LC-MS/MS peptide mapping facilitates identification and quantification of specific deamidation sites on the affinity ligand.

Dramatic Shape Changes Occur as Cytochrome c Folds

Extensive experimental studies on the folding of cytochrome c (Cyt c) make this small protein an ideal target for atomic detailed simulations for the purposes of quantitatively characterizing the structural transitions and the associated time scales for folding to the native state from an ensemble of unfolded states. We use previously generated atomically detailed folding trajectories by the stochastic difference equation in length to calculate the time-dependent changes in the small-angle X-ray scattering (SAXS) profiles. Excellent agreement is obtained between experiments and simulations for the time-dependent SAXS spectra, allowing us to identify the structures of the folding intermediates, which shows that Cyt c reaches the native state by a sequential folding mechanism. Using the ensembles of structures along the folding pathways, we show that compaction and the sphericity of Cyt c change dramatically from the prolate ellipsoid shape in the unfolded state to the spherical native state. Our data, which are in unprecedented quantitative agreement with all aspects of time-resolved SAXS experiments, show that hydrophobic collapse and amide group protection coincide on the 100 microseconds time scale, which is in accordance with ultrafast hydrogen/deuterium exchange studies. Based on these results, we propose that compaction of polypeptide chains, accompanied by dramatic shape changes, is a universal characteristic of globular proteins, regardless of the underlying folding mechanism.

Transition path theory from biased simulations

Transition Path Theory (TPT) provides a rigorous framework to investigate the dynamics of rare thermally activated transitions. In this theory, a central role is played by the forward committor function q⁺(x), which provides the ideal reaction coordinate. Furthermore, the reactive dynamics and kinetics are fully characterized in terms of two time-independent scalar and vector distributions. In this work, we develop a scheme which enables all these ingredients of TPT to be efficiently computed using the short non-equilibrium trajectories generated by means of a specific combination of enhanced path sampling techniques. In particular, first we further extend the recently introduced self-consistent path sampling algorithm in order to compute the committor q⁺(x). Next, we show how this result can be exploited in order to define efficient algorithms which enable us to directly sample the transition path ensemble.

The effects of implicit modeling of nonpolar solvation on protein folding simulations

Implicit solvent models, in which the polar and nonpolar solvation free-energies of solute molecules are treated separately, have been widely adopted for molecular dynamics simulation of protein folding. While the development of the implicit models is mainly focused on the methodological improvement and key parameter optimization for polar solvation, nonpolar solvation has been either ignored or described by a simplistic surface area (SA) model. In this work, we performed the folding simulations of multiple β-hairpin and α-helical proteins with varied surface tension coefficients embedded in the SA model to clearly demonstrate the effects of nonpolar solvation treated by a popular SA model on protein folding. The results indicate that the change in the surface tension coefficient does not alter the ability of implicit solvent simulations to reproduce a protein native structure but indeed controls the components of the equilibrium conformational ensemble and modifies the energetic characterization of the folding transition pathway. The suitably set surface tension coefficient can yield explicit solvent simulations and/or experimentally suggested folding mechanism of protein. In addition, the implicit treatment of both polar and nonpolar components of solvation free-energy contributes to the overestimation of the secondary structure in implicit solvent simulations.

All-atom calculation of protein free-energy profiles

The Bias Functional (BF) approach is a variational method which enables one to efficiently generate ensembles of reactive trajectories for complex biomolecular transitions, using ordinary computer clusters. For example, this scheme was applied to simulate in atomistic detail the folding of proteins consisting of several hundreds of amino acids and with experimental folding time of several minutes. A drawback of the BF approach is that it produces trajectories which do not satisfy microscopic reversibility. Consequently, this method cannot be used to directly compute equilibrium observables, such as free energy landscapes or equilibrium constants. In this work, we develop a statistical analysis which permits us to compute the potential of mean-force (PMF) along an arbitrary collective coordinate, by exploiting the information contained in the reactive trajectories calculated with the BF approach. We assess the accuracy and computational efficiency of this scheme by comparing its results with the PMF obtained for a small protein by means of plain molecular dynamics.

Augmenting the anisotropic network model with torsional potentials improves PATH performance, enabling detailed comparison with experimental rate data

PATH algorithms for identifying conformational transition states provide computational parameters-time to the transition state, conformational free energy differences, and transition state activation energies-for comparison to experimental data and can be carried out sufficiently rapidly to use in the "high throughput" mode. These advantages are especially useful for interpreting results from combinatorial mutagenesis experiments. This report updates the previously published algorithm with enhancements that improve correlations between PATH convergence parameters derived from virtual variant structures generated by RosettaBackrub and previously published kinetic data for a complete, four-way combinatorial mutagenesis of a conformational switch in Tryptophanyl-tRNA synthetase.

Minimum action transition paths connecting minima on an energy surface

Dynamics is essential to the biological functions of many bio-molecules, yet our knowledge of dynamics remains fragmented. Experimental techniques for studying bio-molecules either provide high resolution information on static conformations of the molecule or provide low-resolution, ensemble information that does not shed light on single molecule dynamics. In parallel, bio-molecular dynamics occur at time scale that are not yet attainable through detailed simulation methods. These limitations are especially noticeable when studying transition paths. To address this issue, we report in this paper two methods that derive meaningful trajectories for proteins between two of their conformations. The first method, MinActionPath, uses approximations of the potential energy surface for the molecule to derive an analytical solution of the equations of motion related to the concept of minimum action path. The second method, RelaxPath, follows the same principle of minimum action path but implements a more sophisticated potential, including a mixed elastic potential and a collision term to alleviate steric clashes. Using this new potential, the equations of motion cannot be solved analytically. We have introduced a relaxation method for solving those equations. We describe both the theories behind the two methods and their implementations, focusing on the specific techniques we have used that make those implementations amenable to study large molecular systems. We have illustrated the performance of RelaxPath on simple 2D systems. We have also compared MinActionPath and RelaxPath to other methods for generating transition paths on a well suited test set of large proteins, for which the end points of the trajectories as well as an intermediate conformation between those end points are known. We have shown that RelaxPath outperforms those other methods, including MinActionPath, in its ability to generate trajectories that get close to the known intermediates. We have also shown that the structures along the RelaxPath trajectories remain protein-like. Open source versions of the two programs MinActionPath and RelaxPath are available by request.

The Renormalization Group and Its Applications to Generating Coarse-Grained Models of Large Biological Molecular Systems

Understanding the dynamics of biomolecules is the key to understanding their biological activities. Computational methods ranging from all-atom molecular dynamics simulations to coarse-grained normal-mode analyses based on simplified elastic networks provide a general framework to studying these dynamics. Despite recent successes in studying very large systems with up to a 100,000,000 atoms, those methods are currently limited to studying small- to medium-sized molecular systems due to computational limitations. One solution to circumvent these limitations is to reduce the size of the system under study. In this paper, we argue that coarse-graining, the standard approach to such size reduction, must define a hierarchy of models of decreasing sizes that are consistent with each other, i.e., that each model contains the information of the dynamics of its predecessor. We propose a new method, Decimate, for generating such a hierarchy within the context of elastic networks for normal-mode analysis. This method is based on the concept of the renormalization group developed in statistical physics. We highlight the details of its implementation, with a special focus on its scalability to large systems of up to millions of atoms. We illustrate its application on two large systems, the capsid of a virus and the ribosome translation complex. We show that highly decimated representations of those systems, containing down to 1% of their original number of atoms, still capture qualitatively and quantitatively their dynamics. Decimate is available as an OpenSource resource.

Perspective: Computer simulations of long time dynamics

Atomically detailed computer simulations of complex molecular events attracted the imagination of many researchers in the field as providing comprehensive information on chemical, biological, and physical processes. However, one of the greatest limitations of these simulations is of time scales. The physical time scales accessible to straightforward simulations are too short to address many interesting and important molecular events. In the last decade significant advances were made in different directions (theory, software, and hardware) that significantly expand the capabilities and accuracies of these techniques. This perspective describes and critically examines some of these advances.

Folding Mechanism of Proteins Im7 and Im9: Insight from All-Atom Simulations in Implicit and Explicit Solvent

Im7 and Im9 are evolutionary related proteins with almost identical native structures. In spite of their structural similarity, experiments show that Im7 folds through a long-lived on-pathway intermediate, while Im9 folds according to two-state kinetics. In this work, we use a recently developed enhanced path sampling method to generate many folding trajectories for these proteins, using realistic atomistic force fields, in both implicit and explicit solvent. Overall, our results are in good agreement with the experimental ϕ values and with the result of ϕ-value-restrained molecular dynamics (MD) simulations. However, our implicit solvent simulations fail to predict a qualitative difference in the folding pathways of Im7 and Im9. In contrast, our simulations in explicit solvent correctly reproduce the fact that only protein Im7 folds through a on-pathway intermediate. By analyzing our atomistic trajectories, we provide a physical picture which explains the observed difference in the folding kinetics of these chains.

A modified PATH algorithm rapidly generates transition states comparable to those found by other well established algorithms

PATH rapidly computes a path and a transition state between crystal structures by minimizing the Onsager-Machlup action. It requires input parameters whose range of values can generate different transition-state structures that cannot be uniquely compared with those generated by other methods. We outline modifications to estimate these input parameters to circumvent these difficulties and validate the PATH transition states by showing consistency between transition-states derived by different algorithms for unrelated protein systems. Although functional protein conformational change trajectories are to a degree stochastic, they nonetheless pass through a well-defined transition state whose detailed structural properties can rapidly be identified using PATH.

Fast helix formation in the B domain of protein A revealed by site-specific infrared probes

Comparison of experimental and computational protein folding studies can be difficult because of differences in structural resolution. Isotope-edited infrared spectroscopy offers a direct measure of structural changes involved in protein folding at the single-residue level. Here we demonstrate the increased resolution of site-specific infrared probes to the peptide backbone in the B domain of staphylococcal protein A (BdpA). (13)C═(18)O-labeled methionine was incorporated into each of the helices using recombinant protein expression. Laser-induced temperature jumps coupled with infrared spectroscopy were used to probe changes in the peptide backbone on the submillisecond time scale. The relaxation kinetics of the buried helices, solvated helices, and labeled positions were measured independently by probing the corresponding bands assigned in the amide I region. Using these wavelength-dependent measurements, we observe a fast nanosecond phase and slower microsecond phase at each position. We find at least partial formation of helices 1-3 in the fast intermediate state that precedes the transition state. These measurements provide direct, time-resolved experimental evidence of the early formation of partial helical structure in helices 1 and 3, supporting folding models proposed by computer simulations.

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.

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.

From coarse-grained to atomic-level characterization of protein dynamics: transition state for the folding of B domain of protein A

Atomic-level molecular dynamics simulations are widely used for the characterization of the structural dynamics of proteins; however, they are limited to shorter time scales than the duration of most of the relevant biological processes. Properly designed coarse-grained models that trade atomic resolution for efficient sampling allow access to much longer time-scales. In-depth understanding of the structural dynamics, however, must involve atomic details. In this study, we tested a method for the rapid reconstruction of all-atom models from α carbon atom positions in the application to convert a coarse-grained folding trajectory of a well described model system: the B domain of protein A. The results show that the method and the spatial resolution of the resulting coarse-grained models enable computationally inexpensive reconstruction of realistic all-atom models. Additionally, by means of structural clustering, we determined the most persistent ensembles of the key folding step, the transition state. Importantly, the analysis of the overall structural topologies suggests a dominant folding pathway. This, together with the all-atom characterization of the obtained ensembles, in the form of contact maps, matches the experimental results well.

Hidden protein folding pathways in free-energy landscapes uncovered by network analysis

A network analysis is used to uncover hidden folding pathways in free-energy landscapes usually defined in terms of such arbitrary order parameters as root-mean-square deviation from the native structure, radius of gyration, etc. The analysis has been applied to molecular dynamics (MD) trajectories of the B-domain of staphylococcal protein A, generated with the coarse-grained united-residue (UNRES) force field in a broad range of temperatures (270K ≤ T ≤ 325K). Thousands of folding pathways have been identified at each temperature. Out of these many folding pathways, several most probable ones were selected for investigation of the conformational transitions during protein folding. Unlike other conformational space network (CSN) methods, a node in the CSN variant implemented in this work is defined according to the nativelikeness class of the structure, which defines the similarity of segments of the compared structures in terms of secondary-structure, contact-pattern, and local geometry, as well as the overall geometric similarity of the conformation under consideration to that of the reference (experimental) structure. Our previous findings, regarding the folding model and conformations found at the folding-transition temperature for protein A (Maisuradze et al., J. Am. Chem. Soc. 132, 9444, 2010), were confirmed by the conformational space network analysis. In the methodology and in the analysis of the results, the shortest path identified by using the shortest-path algorithm corresponds to the most probable folding pathway in the conformational space network.

Dominant folding pathways of a WW domain

We investigate the folding mechanism of the WW domain Fip35 using a realistic atomistic force field by applying the Dominant Reaction Pathways approach. We find evidence for the existence of two folding pathways, which differ by the order of formation of the two hairpins. This result is consistent with the analysis of the experimental data on the folding kinetics of WW domains and with the results obtained from large-scale molecular dynamics simulations of this system. Free-energy calculations performed in two coarse-grained models support the robustness of our results and suggest that the qualitative structure of the dominant paths are mostly shaped by the native interactions. Computing a folding trajectory in atomistic detail only required about one hour on 48 Central Processing Units. The gain in computational efficiency opens the door to a systematic investigation of the folding pathways of a large number of globular proteins.

Fluctuations in the ensemble of reaction pathways

The dominant reaction pathway is a rigorous framework to microscopically compute the most probable trajectories, in nonequilibrium transitions. In the low-temperature regime, such dominant pathways encode the information about the reaction mechanism and can be used to estimate nonequilibrium averages of arbitrary observables. On the other hand, at sufficiently high temperatures, the stochastic fluctuations around the dominant paths become important and have to be taken into account. In this work, we develop a technique to systematically include the effects of such stochastic fluctuations, to order k(B)T. This method is used to compute the probability for a transition to take place through a specific reaction channel and to evaluate the reaction rate.

Simulation of chaperonin effect on protein folding: a shift from nucleation-condensation to framework mechanism

The iterative annealing mechanism (IAM) of chaperonin-assisted protein folding is explored in a framework of a well-established coarse-grained protein modeling tool, which enables the study of protein dynamics in a time-scale well beyond classical all-atom molecular mechanics. The chaperonin mechanism of action is simulated for two paradigm systems of protein folding, B domain of protein A (BdpA) and B1 domain of protein G (GB1), and compared to chaperonin-free simulations presented here for BdpA and recently published for GB1. The prediction of the BdpA transition state ensemble (TSE) is in perfect agreement with experimental findings. It is shown that periodic distortion of the polypeptide chains by hydrophobic chaperonin interactions can promote rapid folding and leads to a decrease in folding temperature. It is also demonstrated how chaperonin action prevents kinetically trapped conformations and modulates the observed folding mechanisms from nucleation-condensation to a more framework-like.

Molecular dynamics at low time resolution

The internal dynamics of macromolecular systems is characterized by widely separated time scales, ranging from fraction of picoseconds to nanoseconds. In ordinary molecular dynamics simulations, the elementary time step Δt used to integrate the equation of motion needs to be chosen much smaller of the shortest time scale in order not to cut-off physical effects. We show that in systems obeying the overdamped Langevin equation, it is possible to systematically correct for such discretization errors. This is done by analytically averaging out the fast molecular dynamics which occurs at time scales smaller than Δt, using a renormalization group based technique. Such a procedure gives raise to a time-dependent calculable correction to the diffusion coefficient. The resulting effective Langevin equation describes by construction the same long-time dynamics, but has a lower time resolution power, hence it can be integrated using larger time steps Δt. We illustrate and validate this method by studying the diffusion of a point-particle in a one-dimensional toy model and the denaturation of a protein.

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.

Design of an optimized scaffold for affibody molecules

Affibody molecules are non-immunoglobulin-derived affinity proteins based on a three-helical bundle protein domain. Here, we describe the design process of an optimized Affibody molecule scaffold with improved properties and a surface distinctly different from that of the parental scaffold. The improvement was achieved by applying an iterative process of amino acid substitutions in the context of the human epidermal growth factor receptor 2 (HER2)-specific Affibody molecule Z(HER2:342). Replacements in the N-terminal region, loop 1, helix 2 and helix 3 were guided by extensive structural modeling using the available structures of the parent Z domain and Affibody molecules. The effect of several single substitutions was analyzed followed by combination of up to 11 different substitutions. The two amino acid substitutions N23T and S33K accounted for the most dramatic improvements, including increased thermal stability with elevated melting temperatures of up to +12 degrees C. The optimized scaffold contains 11 amino acid substitutions in the nonbinding surface and is characterized by improved thermal and chemical stability, as well as increased hydrophilicity, and enables generation of identical Affibody molecules both by chemical peptide synthesis and by recombinant bacterial expression. A HER2-specific Affibody tracer, [MMA-DOTA-Cys61]-Z(HER2:2891)-Cys (ABY-025), was produced by conjugating MMA-DOTA (maleimide-monoamide-1,4,7,10-tetraazacyclododecane-1,4,7,10-tetraacetic acid) to the peptide produced either chemically or in Escherichia coli. ABY-025 showed high affinity and specificity for HER2 (equilibrium dissociation constant, K(D), of 76 pM) and detected HER2 in tissue sections of SKOV-3 xenograft and human breast tumors. The HER2-binding capacity was fully retained after three cycles of heating to 90 degrees C followed by cooling to room temperature. Furthermore, the binding surfaces of five Affibody molecules targeting other proteins (tumor necrosis factor alpha, insulin, Taq polymerase, epidermal growth factor receptor or platelet-derived growth factor receptor beta) were grafted onto the optimized scaffold, resulting in molecules with improved thermal stability and a more hydrophilic nonbinding surface.

How conformational transition depends on hydrophobicity of elastin-like polypeptides

The three-dimensional structures of elastin-like polypeptides Val1-Pro2-Gly3-Xaa4-Gly5 were investigated by using the multicanonical Monte Carlo (MC) simulation procedure. By substituting different amino acids in the fourth position of the sequence, the thermodynamical variables are calculated in vacuo and in solvent to determine the hydrophobicity dependence of the conformational transition temperatures of the peptides. Resultant hydrophobicity scale is in good agreement with many hydrophobicity scales.

Relation between free energy landscapes of proteins and dynamics

By examining the molecular dynamics (MD) of protein folding trajectories, generated with the coarse-grained UNRES force field, for the B-domain of staphylococcal protein A and the triple β-strand WW domain from the Formin binding protein 28 (FBP), by principal component analysis (PCA), it is demonstrated how different free energy landscapes (FELs) and folding pathways of trajectories can be, even though they appear to be very similar by visual inspection of the time-dependence of the root-mean-square deviation (rmsd). Approaches to determine the minimal dimensionality of FELs for a correct description of protein folding dynamics are discussed. The correlation between the amplitude of the fluctuations of proteins and the dimensionality of the FELs is shown. The advantage of internal coordinate PCA over Cartesian PCA for small proteins is also illustrated.

Nearly symmetrical proteins: folding pathways and transition states

The folding pathways of the B domain of protein A have been the subject of many experimental and computational studies. Based on a statistical mechanical model, it has been suggested that the native state symmetry leads to multiple pathways, highly dependent on temperature and denaturant concentration. Experiments, however, have not confirmed this scenario. By considering four nearly symmetrical proteins, one of them being the above molecule, here we show that, if contact energies are properly taken into account, a different picture emerges from kinetic simulations of the above-mentioned model. This is characterized by a dominant folding pathway, which is consistent with the most recent experimental results. Given the simplicity of the model, we also report on a direct sampling of the transition state.

Free-energy-driven folding and thermodynamics of the 67-residue protein GS-alpha3W--a large-scale Monte Carlo study

Utilizing the computational power of a few thousand processors on a BlueGene/P, we have explored the folding mechanism of the 67-residue protein GS-alpha(3)W. Results from our large-scale simulation indicate a diffusion-collision mechanism for folding. However, the lower-than-expected frequency of native-like configurations at physiological temperatures indicates shortcomings of our energy function. Our results suggest that computational studies of large proteins call for redevelopment and reparametrization of force fields that in turn require extensive simulations only possible with the newly available supercomputers with computing powers reaching the petaflop range.

Path-integral method for predicting relative binding affinities of protein-ligand complexes

We present a novel approach for computing biomolecular interaction binding affinities based on a simple path integral solution of the Fokker-Planck equation. Computing the free energy of protein-ligand interactions can expedite structure-based drug design. Traditionally, the problem is seen through the lens of statistical thermodynamics. The computations can become, however, prohibitively long for the change in the free energy upon binding to be determined accurately. In this work, we present a different approach based on a stochastic kinetic formalism. Inspired by Feynman's path integral formulation, we extend the theory to classical interacting systems. The ligand is modeled as a Brownian particle subjected to the effective nonbonding interaction potential of the receptor. This allows the calculation of the relative binding affinities of interacting biomolecules in water to be computed as a function of the ligand's diffusivity and the curvature of the potential surface in the vicinity of the binding minimum. The calculation is thus exceedingly rapid. In test cases, the correlation coefficient between actual and computed free energies is >0.93 for accurate data sets.

Dominant reaction pathways in high-dimensional systems

This paper is devoted to the development of a theoretical and computational framework denominated dominant reaction pathways (DRPs) to efficiently sample the statistically significant thermally activated reaction pathways, in multidimensional systems. The DRP approach is consistently derived from the Langevin equation through a systematic expansion in the thermal energy, k(B)T. Its main advantage with respect to existing simulation techniques is that it provides a natural and rigorous framework to perform the path sampling using constant displacement steps, rather than constant time steps. In our previous work, we have shown how to obtain the set of most probable reaction pathways, i.e., the lowest order in the k(B)T expansion. In this work, we show how to compute the corrections to the leading order due to stochastic fluctuations around the most probable trajectories. We also discuss how to obtain predictions for the evolution of arbitrary observables and how to generate conformations, which are representative of the transition state ensemble. We illustrate how our method works in practice by studying the diffusion of a point particle in a two-dimensional funneled external potential.

Quantifying the structural requirements of the folding transition state of protein A and other systems

The B-domain of protein A is a small three-helix bundle that has been the subject of considerable experimental and theoretical investigation. Nevertheless, a unified view of the structure of the transition-state ensemble (TSE) is still lacking. To characterize the TSE of this surprisingly challenging protein, we apply a combination of psi analysis (which probes the role of specific side-chain to side-chain contacts) and kinetic H/D amide isotope effects (which measures hydrogen-bond content), building upon previous studies using mutational phi analysis (which probes the energetic influence of side-chain substitutions). The second helix is folded in the TSE, while helix formation appears just at the carboxy and amino termini of the first and third helices, respectively. The experimental data suggest a homogenous yet plastic TS with a native-like topology. This study generalizes our earlier conclusion, based on two larger alpha/beta proteins, that the TSEs of most small proteins achieve approximately 70% of their native state's relative contact order. This high percentage limits the degree of possible TS heterogeneity and requires a reevaluation of the structural content of the TSE of other proteins, especially when they are characterized as small or polarized.

From helix-coil transitions to protein folding

An evolution of procedures to simulate protein structure and folding pathways is described. From an initial focus on the helix-coil transition and on hydrogen-bonding and hydrophobic interactions, our original attempts to determine protein structure and folding pathways were based on an experimental approach. Experiments on the oxidative folding of reduced bovine pancreatic ribonuclease A (RNase A) led to a mechanism by which the molecule folded to the native structure by a minimum of four different pathways. The experiments with RNase A were followed by development of a molecular mechanics approach, first, making use of global optimization procedures and then with molecular dynamics (MD), evolving from an all-atom to a united-residue model. This hierarchical MD approach facilitated probing of the folding trajectory to longer time scales than with all-atom MD, and hence led to the determination of complete folding trajectories, thus far for a protein containing as many as 75 amino acid residues. With increasing refinement of the computational procedures, the computed results are coming closer to experimental observations, providing an understanding as to how physics directs the folding process.

Computer simulations of the refolding of sperm whale apomyoglobin from high-temperature denaturated state

The refolding mechanism of apomyoglobin (apoMb) subsequent to high-temperature unfolding has been examined using computer simulations with atomic level detail. The folding of this protein has been extensively studied experimentally, providing a large database of folding parameters which can be probed using simulations. In the present study, 4-folding trajectories of apoMb were computed starting from coiled structures. A crystal structure of sperm whale myoglobin taken from the Protein Data Bank was used to construct the final native conformation by removal of the heme group followed by energy optimization. The initial unfolded conformations were obtained from high-temperature molecular dynamics simulations. Room-temperature refolding trajectories at neutral pH were obtained using the stochastic difference equation in length algorithm. The folding trajectories were compared with experimental results and two previous molecular dynamics studies at low pH. In contrast to the previous simulations, an extended intermediate with large helical content was not observed. In the present study, a structural collapse occurs without formation of helices or native contacts. Once the protein structure is more compact (radius of gyration<18 A) secondary and tertiary structures appear. These results suggest that apoMb follows a different folding pathway after high-temperature denaturation.

The complex folding pathways of protein A suggest a multiple-funnelled energy landscape

Folding proteins into their native states requires the formation of both secondary and tertiary structures. Many questions remain, however, as to whether these form into a precise order, and various pictures have been proposed that place the emphasis on the first or the second level of structure in describing folding. One of the favorite test models for studying this question is the B domain of protein A, which has been characterized by numerous experiments and simulations. Using the activation-relaxation technique coupled with a generic energy model (optimized potential for efficient peptide structure prediction), we generate more than 50 folding trajectories for this 60-residue protein. While the folding pathways to the native state are fully consistent with the funnel-like description of the free energy landscape, we find a wide range of mechanisms in which secondary and tertiary structures form in various orders. Our nonbiased simulations also reveal the presence of a significant number of non-native beta and alpha conformations both on and off pathway, including the visit, for a non-negligible fraction of trajectories, of fully ordered structures resembling the native state of nonhomologous proteins.

Backbone and side-chain ordering in a small protein

We investigate the relation between backbone and side-chain ordering in a small protein. For this purpose, we have performed multicanonical simulations of the villin headpiece subdomain HP-36, an often used toy model in protein studies. Concepts of circular statistics are introduced to analyze side-chain fluctuations. In contrast to earlier studies on homopolypeptides [Wei et al., J. Phys. Chem. B 111, 4244 (2007)], we do not find collective effects leading to a separate transition. Rather, side-chain ordering is spread over a wide temperature range. Our results indicate a thermal hierarchy of ordering events, with side-chain ordering appearing at temperatures below the helix-coil transition but above the folding transition. We conjecture that this thermal hierarchy reflects an underlying temporal order, and that side-chain ordering facilitates the search for the correct backbone topology.

An analytical study of the interplay between geometrical and energetic effects in protein folding

Analytical studies have several advantages for an understanding of the mechanisms of protein folding such as the interplay between geometrical and energetic effects. In this paper, we introduce a Gaussian filament with a C(alpha) structure-based (Go) potential as a new theoretical scheme based on a Hamiltonian approach. This model takes into account geometrical information in a realistic fashion without the need of phenomenological descriptions. In order to make this model more appropriate for comparison with protein folding simulations and experiments, we introduce a many-body interaction into the potential term to enhance cooperativity. We apply our new analytical model to a beta-hairpin-type peptide and compare our results with a molecular dynamics simulation of a structure-based model.

Simulating temperature jumps for protein folding

We present a new computational methodology aimed to calculate the thermodynamics and kinetics of peptide folding. We focus in particular on temperature jump experiments of folding rates and show how a combination of replica exchange molecular dynamics (REMD) followed by multiplexed molecular dynamics starting from structures taken from the REMD runs can be used to extract properties in line with experiments. A model system, alanine20, is studied in this article as a proof of principle and used to describe the methodology.

Universality and diversity of folding mechanics for three-helix bundle proteins

In this study we evaluate, at full atomic detail, the folding processes of two small helical proteins, the B domain of protein A and the Villin headpiece. Folding kinetics are studied by performing a large number of ab initio Monte Carlo folding simulations using a single transferable all-atom potential. Using these trajectories, we examine the relaxation behavior, secondary structure formation, and transition-state ensembles (TSEs) of the two proteins and compare our results with experimental data and previous computational studies. To obtain a detailed structural information on the folding dynamics viewed as an ensemble process, we perform a clustering analysis procedure based on graph theory. Moreover, rigorous p(fold) analysis is used to obtain representative samples of the TSEs and a good quantitative agreement between experimental and simulated Phi values is obtained for protein A. Phi values for Villin also are obtained and left as predictions to be tested by future experiments. Our analysis shows that the two-helix hairpin is a common partially stable structural motif that gets formed before entering the TSE in the studied proteins. These results together with our earlier study of Engrailed Homeodomain and recent experimental studies provide a comprehensive, atomic-level picture of folding mechanics of three-helix bundle proteins.

Protein folding: then and now

Over the past three decades the protein folding field has undergone monumental changes. Originally a purely academic question, how a protein folds has now become vital in understanding diseases and our abilities to rationally manipulate cellular life by engineering protein folding pathways. We review and contrast past and recent developments in the protein folding field. Specifically, we discuss the progress in our understanding of protein folding thermodynamics and kinetics, the properties of evasive intermediates, and unfolded states. We also discuss how some abnormalities in protein folding lead to protein aggregation and human diseases.

Optimized folding simulations of protein A

We describe optimized parallel tempering simulations of the 46-residue B-fragment of protein A. Native-like configurations with a root-mean-square deviation of approximately 3 A to the experimentally determined structure (Protein Data Bank identifier 1BDD) are found. However, at biologically relevant temperatures such conformations appear with only approximately 10 % frequency in our simulations. Possible shortcomings in our energy function are discussed.

Revisiting the carboxylic acid dimers in aqueous solution: interplay of hydrogen bonding, hydrophobic interactions, and entropy

Carboxylic acid dimers are useful model systems for understanding the interplay of hydrogen bonding, hydrophobic effects, and entropy in self-association and assembly. Through extensive sampling with a classical force field and careful free energy analysis, it is demonstrated that both hydrogen bonding and hydrophobic interactions are indeed important for dimerization of carboxylic acids (except formic acid). The dimers are only weakly ordered, and the degree of ordering increases with stronger hydrophobic interactions between longer alkyl chains. Comparison of calculated and experimental dimerization constants reveals a systematic tendency for excessive self-aggregation in current classical force fields. Qualitative and quantitative information on the thermodynamics of hydrogen bonding and hydrophobic interactions derived from these simulations is in excellent agreement with existing results from experiment and theory. These results provide a verification from first principles of previous estimations based on two statistical mechanical hydrophobic theories. We also revisit and clarify the fundamental statistical thermodynamics formalism for calculating absolute binding constants, external entropy, and solvation entropy changes upon association from detailed free energy simulations. This analysis is believed to be useful for a wide range of applications including computational studies of protein-ligand and protein-protein binding.

MinActionPath: maximum likelihood trajectory for large-scale structural transitions in a coarse-grained locally harmonic energy landscape

The non-linear problem of simulating the structural transition between two known forms of a macromolecule still remains a challenge in structural biology. The problem is usually addressed in an approximate way using 'morphing' techniques, which are linear interpolations of either the Cartesian or the internal coordinates between the initial and end states, followed by energy minimization. Here we describe a web tool that implements a new method to calculate the most probable trajectory that is exact for harmonic potentials; as an illustration of the method, the classical Calpha-based Elastic Network Model (ENM) is used both for the initial and the final states but other variants of the ENM are also possible. The Langevin equation under this potential is solved analytically using the Onsager and Machlup action minimization formalism on each side of the transition, thus replacing the original non-linear problem by a pair of linear differential equations joined by a non-linear boundary matching condition. The crossover between the two multidimensional energy curves around each state is found numerically using an iterative approach, producing the most probable trajectory and fully characterizing the transition state and its energy. Jobs calculating such trajectories can be submitted on-line at: http://lorentz.dynstr.pasteur.fr/joel/index.php.