Peptide chain dynamics in light and heavy water: zooming in on internal friction

Frictional effects due to the chain itself, rather than the solvent, may have a significant effect on protein dynamics. Experimentally, such "internal friction" has been investigated by studying folding or binding kinetics at varying solvent viscosity; however, the molecular origin of these effects is hard to pinpoint. We consider the kinetics of disordered glycine-serine and α-helix forming alanine peptides and a coarse-grained protein folding model in explicit-solvent molecular dynamics simulations. By varying the solvent mass over more than two orders of magnitude, we alter only the solvent viscosity and not the folding free energy. Folding dynamics at the near-vanishing solvent viscosities accessible by this approach suggests that solvent and internal friction effects are intrinsically entangled. This finding is rationalized by calculation of the polymer end-to-end distance dynamics from a Rouse model that includes internal friction. An analysis of the friction profile along different reaction coordinates, extracted from the simulation data, demonstrates that internal as well as solvent friction varies substantially along the folding pathways and furthermore suggests a connection between friction and the formation of hydrogen bonds upon folding.     

Solvent viscosity dependence of the folding rate of a small protein: distributed computing study

By using distributed computing techniques and a supercluster of more than 20,000 processors we simulated folding of a 20-residue Trp Cage miniprotein in atomistic detail with implicit GB/SA solvent at a variety of solvent viscosities (gamma). This allowed us to analyze the dependence of folding rates on viscosity. In particular, we focused on the low-viscosity regime (values below the viscosity of water). In accordance with Kramers' theory, we observe approximately linear dependence of the folding rate on 1/gamma for values from 1-10(-1)x that of water viscosity. However, for the regime between 10(-4)-10(-1)x that of water viscosity we observe power-law dependence of the form k approximately gamma(-1/5). These results suggest that estimating folding rates from molecular simulations run at low viscosity under the assumption of linear dependence of rate on inverse viscosity may lead to erroneous results.     

Numerical simulation of the effect of solvent viscosity on the motions of a beta-peptide heptamer

This report examines the effect of a decrease in solvent viscosity on the simulated folding behaviour of a beta-peptide heptamer in methanol. Simulations of the molecular dynamics of the heptamer H-beta3-HVal-beta3-HAla-beta3-HLeu-(S,S)-beta3-HAla(alphaMe)-beta3-HVal-beta3-HAla-beta3-HLeu-OH in methanol, with an explicit representation of the methanol molecules, were performed for 80 ns at various solvent viscosities. The simulations indicate that at a solvent viscosity of one third of that of methanol, only the dynamic aspects of the folding process are altered, and that the rate of folding is increased. At a viscosity of one tenth of that of methanol, insufficient statistics are obtained within the 80 ns period. We suggest that 80 ns is an insufficient time to reach conformational equilibrium at very low viscosity because the dependence of the folding rate of a beta-peptide on solvent viscosity has two regimes; a result that was observed in another computational study for alpha-peptides.     

The influence of thermostats and manostats on reverse nonequilibrium molecular dynamics calculations of fluid viscosities

Reverse nonequilibrium molecular dynamics to calculate the shear viscosity of Lennard-Jones liquids was extended to simulations at constant number of particles, constant volume, and constant pressure using a Berendsen thermostat and a Berendsen manostat. Using additional systems such as water and hexane, we also report on the performance of shear viscosity calculations of systems with electrostatic and nontrivial intramolecular interactions when a manostat is applied. We compare the shear viscosities of simulations using no coupling, only temperature coupling, and temperature and pressure coupling and characterize discrepancies, where observed. From this, we deduce guidelines for when and how manostats can be usefully applied in reverse nonequilibrium simulations.     

Pulling direction as a reaction coordinate for the mechanical unfolding of single molecules

The folding and unfolding kinetics of single molecules, such as proteins or nucleic acids, can be explored by mechanical pulling experiments. Determining intrinsic kinetic information, at zero stretching force, usually requires an extrapolation by fitting a theoretical model. Here, we apply a recent theoretical approach describing molecular rupture in the presence of force to unfolding kinetic data obtained from coarse-grained simulations of ubiquitin. Unfolding rates calculated from simulations over a broad range of stretching forces, for different pulling directions, reveal a remarkable "turnover" from a force-independent process at low force to a force-dependent process at high force, akin to the "roll-over" in unfolding rates sometimes seen in studies using chemical denaturant. While such a turnover in rates is unexpected in one dimension, we demonstrate that it can occur for dynamics in just two dimensions. We relate the turnover to the quality of the pulling direction as a reaction coordinate for the intrinsic folding mechanism. A novel pulling direction, designed to be the most relevant to the intrinsic folding pathway, results in the smallest turnover. Our results are in accord with protein engineering experiments and simulations which indicate that the unfolding mechanism at high force can differ from the intrinsic mechanism. The apparent similarity between extrapolated and intrinsic rates in experiments, unexpected for different unfolding barriers, can be explained if the turnover occurs at low forces.     

Viscosity scaling for the glassy phase of protein folding

Although commendable progress has been made in the understanding of the physics of protein folding, a key unresolved issue is whether Kramers' diffusion model of chemical reactions is generally applicable to activated barrier crossing events during folding. To examine the solvent viscosity effect on the folding transition of native-like trapped intermediates, laser flash photolysis has been used to measure the microsecond folding kinetics of a natively folded state of CO-liganded ferrocytochrome c (M-state) in the 1-250 cP range of glycerol viscosity at pH 7.0, 20 degrees C. The single rate coefficient for the folding of the M-state to the native state of the protein (i.e., the M --> N folding process) decreases initially when the solvent viscosity is low (<10 cP), but saturates at higher viscosity, indicating that Kramers model is not general enough for scaling the viscosity dependence of post-transition folding involving glassy dynamics. Analysis based on the Grote-Hynes idea of time dependent friction in conjunction with defect diffusion dynamics can account for the observed non-Kramers scaling.     

Turn-directed folding dynamics of β-hairpin-forming de novo decapeptide Chignolin

Realistic mechanistic pictures of β-hairpin formation, offering valuable insights into some of the key early events in protein folding, are accessible through short designed polypeptides as they allow atomic-level scrutiny through simulations. Here, we present a detailed picture of the dynamics and mechanism of β-hairpin formation of Chignolin, a de novo decapeptide, using extensive, unbiased molecular dynamics simulations. The results provide clear evidence for turn-directed broken-zipper folding and reveal details of turn nucleation and cooperative progression of turn growth, hydrogen-bond formations, and eventual packing of the hydrophobic core. Further, we show that, rather than driving folding through hydrophobic collapse, cross-strand side-chain packing could in fact be rate-limiting as packing frustrations can delay formation of the native hydrophobic core prior to or during folding and even cause relatively long-living misfolded or partially folded states that may nucleate aggregative events in more complex situations. The results support the increasing evidence for turn-centric folding mechanisms for β-hairpin formation suggested recently for GB1 and Peptide 1 based on experiments and simulations but also point to the need for similar examinations of polypeptides with larger numbers of cross-strand hydrophobic residues.     

Solvent effect on the folding dynamics and structure of E6-associated protein characterized from ab initio protein folding simulations

Solvent effect on protein conformation and folding mechanism of E6-associated protein (E6ap) peptide are investigated using a recently developed charge update scheme termed as adaptive hydrogen bond-specific charge (AHBC). On the basis of the close agreement between the calculated helix contents from AHBC simulations and experimental results, we observed based on the presented simulations that the two ends of the peptide may simultaneously take part in the formation of the helical structure at the early stage of folding and finally merge to form a helix with lowest backbone RMSD of about 0.9 A? in 40% 2,2,2-trifluoroethanol solution. However, in pure water, the folding may start at the center of the peptide sequence instead of at the two opposite ends. The analysis of the free energy landscape indicates that the solvent may determine the folding clusters of E6ap, which subsequently leads to the different final folded structure. The current study demonstrates new insight to the role of solvent in the determination of protein structure and folding dynamics.     

Diffusion models of protein folding

In theory and in the analysis of experiments, protein folding is often described as diffusion along a single coordinate. We explore here the application of a one-dimensional diffusion model to interpret simulations of protein folding, where the parameters of a model that "best" describes the simulation trajectories are determined using a Bayesian analysis. We discuss the requirements for such a model to be a good approximation to the global dynamics, and several methods for testing its accuracy. For example, one test considers the effect of an added bias potential on the fitted free energies and diffusion coefficients. Such a bias may also be used to extend our approach to determining parameters for the model to systems that would not normally explore the full coordinate range on accessible time scales. Alternatively, the propagators predicted from the model at different "lag" times may be compared with observations from simulation. We then present some applications of the model to protein folding, including Kramers-like turnover in folding rates of coarse-grained models, the effect of non-native interactions on folding, and the effect of the chosen coordinate on the observed position-dependence of the diffusion coefficients. Lastly, we consider how our results are useful for the interpretation of experiments, and how this type of Bayesian analysis may eventually be applied directly to analyse experimental data.     

The influence of macromolecular crowding on HIV-1 protease internal dynamics

High macromolecular concentrations, or crowded conditions, have been shown to affect a wide variety of molecular processes, including diffusion, association and dissociation, and protein folding and stability. Here, we model the effect of macromolecular crowding on the internal dynamics of a protein, HIV-1 protease, using Brownian dynamics simulations. HIV-1 protease possesses a pair of flaps which are postulated to open in the early stages of its catalytic mechanism. Compared to low concentrations, close-packed concentrations of repulsive crowding agents are found to significantly reduce the fraction of time that the protease flaps are open. Macromolecular crowding is likely to have a major effect on in vivo enzyme activity, and may play an important regulatory role in the viral life cycle.     

Variationally determined free energy profiles for structural models of proteins: characteristic temperatures for folding and trapping

Characterizing the phase diagram for proteins is important both for laboratory studies and for the development of structure prediction algorithms. Using a variational scheme, we calculated the generic features of the protein thermostability over a large range of temperatures for a set of more than 50 different proteins using a model based on native structure alone. Focusing on a specific system, protein G, we further examined, using a more realistic model that includes the nonnative interaction, the thermostability of both the native state and a collection of trap structures. By surveying the native structures for many proteins and by paying closer attention to the various trap structures of protein G, we obtained an overall understanding of the folding dynamics far from the conditions usually focused on; namely, those near the folding temperature alone. Two characteristic temperatures (shown to scale with folding temperature in general) signal drastic changes in the folding mechanism. The variational calculations suggest that most proteins would, indeed, fold in a barrierless manner below a critical temperature analogous to a spinodal in crystallization. For fixed interaction strengths, this temperature, however, seems to be generally very low, approximately 50% of the equilibrium folding temperature. Likewise, native proteins, in general, would unfold in a completely barrierless way at a temperature 25% above folding temperature according to these variational calculations. We also studied the distribution of free energy profiles for escape from a set of trap structures generated by simulations.     

Rotational and spin viscosities of water: Application to nanofluidics

In this paper we evaluate the rotational viscosity and the two spin viscosities for liquid water using equilibrium molecular dynamics. Water is modeled via the flexible SPC/Fw model where the Coulomb interactions are calculated via the Wolf method which enables the long simulation times required. We find that the rotational viscosity is independent of the temperature in the range from 284 to 319 K. The two spin viscosities, on the other hand, decrease with increasing temperature and are found to be two orders of magnitude larger than that estimated by Bonthuis et al. [Phys. Rev. Lett. 103, 144503 (2009)] We apply the results from molecular dynamics simulations to the extended Navier-Stokes equations that include the coupling between intrinsic angular momentum and linear momentum. For a flow driven by an external field the coupling will reduce the flow rate significantly for nanoscale geometries. The coupling also enables conversion of rotational electrical energy into fluid linear momentum and we find that in order to obtain measurable flow rates the electrical field strength must be in the order of 0.1?MV?m(-1) and rotate with a frequency of more than 100 MHz.     

Numerical Simulation of the Effect of Solvent Viscosity on the Motions of a β‐Peptide Heptamer

This report examines the effect of a decrease in solvent viscosity on the simulated folding behaviour of a β‐peptide heptamer in methanol. Simulations of the molecular dynamics of the heptamer H‐β³‐HVal‐β³‐HAla‐β³‐HLeu‐(S,S)‐β³‐HAla(αMe)‐β³‐HVal‐β³‐HAla‐β³‐HLeu‐OH in methanol, with an explicit representation of the methanol molecules, were performed for 80 ns at various solvent viscosities. The simulations indicate that at a solvent viscosity of one third of that of methanol, only the dynamic aspects of the folding process are altered, and that the rate of folding is increased. At a viscosity of one tenth of that of methanol, insufficient statistics are obtained within the 80 ns period. We suggest that 80 ns is an insufficient time to reach conformational equilibrium at very low viscosity because the dependence of the folding rate of a β‐peptide on solvent viscosity has two regimes; a result that was observed in another computational study for α‐peptides.     

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.     

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.     

Folding dynamics of Trp-cage in the presence of chemical interference and macromolecular crowding. I

Proteins fold and function in the crowded environment of the cell's interior. In the recent years it has been well established that the so-called "macromolecular crowding" effect enhances the folding stability of proteins by destabilizing their unfolded states for selected proteins. On the other hand, chemical and thermal denaturation is often used in experiments as a tool to destabilize a protein by populating the unfolded states when probing its folding landscape and thermodynamic properties. However, little is known about the complicated effects of these synergistic perturbations acting on the kinetic properties of proteins, particularly when large structural fluctuations, such as protein folding, have been involved. In this study, we have first investigated the folding mechanism of Trp-cage dependent on urea concentration by coarse-grained molecular simulations where the impact of urea is implemented into an energy function of the side chain and/or backbone interactions derived from the all-atomistic molecular dynamics simulations with urea through a Boltzmann inversion method. In urea solution, the folding rates of a model miniprotein Trp-cage decrease and the folded state slightly swells due to a lack of contact formation between side chains at the terminal regions. In addition, the equilibrium m-values of Trp-cage from the computer simulations are in agreement with experimental measurements. We have further investigated the combined effects of urea denaturation and macromolecular crowding on Trp-cage's folding mechanism where crowding agents are modeled as hard-spheres. The enhancement of folding rates of Trp-cage is most pronounced by macromolecular crowding effect when the extended conformations of Trp-cast dominate at high urea concentration. Our study makes quantitatively testable predictions on protein folding dynamics in a complex environment involving both chemical denaturation and macromolecular crowding effects.     

Coupling between lysozyme and glycerol dynamics: microscopic insights from molecular-dynamics simulations

We explore possible molecular mechanisms behind the coupling of protein and solvent dynamics using atomistic molecular-dynamics simulations. For this purpose, we analyze the model protein lysozyme in glycerol, a well-known protein-preserving agent. We find that the dynamics of the hydrogen bond network between the solvent molecules in the first shell and the surface residues of the protein controls the structural relaxation (dynamics) of the whole protein. Specifically, we find a power-law relationship between the relaxation time of the aforementioned hydrogen bond network and the structural relaxation time of the protein obtained from the incoherent intermediate scattering function. We demonstrate that the relationship between the dynamics of the hydrogen bonds and the dynamics of the protein appears also in the dynamic transition temperature of the protein. A study of the dynamics of glycerol as a function of the distance from the surface of the protein indicates that the viscosity seen by the protein is not the one of the bulk solvent. The presence of the protein suppresses the dynamics of the surrounding solvent. This implies that the protein sees an effective viscosity higher than the one of the bulk solvent. We also found significant differences in the dynamics of surface and core residues of the protein. The former is found to follow the dynamics of the solvent more closely than the latter. These results allowed us to propose a molecular mechanism for the coupling of the solvent-protein dynamics.     

Improved theoretical description of protein folding kinetics from rotations in the phase space of relevant order parameters

A method is introduced to construct a better approximation for the reaction coordinate for protein folding from known order parameters. The folding of a two-state off-lattice alpha helical Go-type protein is studied using molecular dynamics simulations. Folding times are computed directly from simulation, as well as theoretically using an equation derived by considering Brownian-type dynamics for the putative reaction coordinate. Theoretical estimates of the folding time using the number of native contacts (Qn) as a reaction coordinate were seen to differ quite significantly from the true folding time of the protein. By considering the properties of the bimodal free energy surface of this protein as a function of Qn and another relevant coordinate for folding Q (the total number of contacts), we show that by introducing a rotation in the phase space of the order parameters Q and Qn, we can construct a new reaction coordinate q that leads to a fivefold improvement in the estimate of the folding rate. This new coordinate q, resulting from the rotation, lies along the line connecting the unfolded and folded ensemble minima of the free energy map plotted as a function of the original order parameters Q and Qn. Possible reasons for the remaining discrepancy between the folding time computed theoretically and from folding simulations are discussed.     

Application of the accelerated molecular dynamics simulations to the folding of a small protein

In this paper, we further explore the applicability of the accelerated molecular dynamics simulation method using a bias potential. The method is applied to both simple model systems and real multidimensional systems. The method is also compared to replica exchange simulations in folding a small protein, Trp cage, using an all atom potential for the protein and an implicit model for the solvent. We show that the bias potential method allows quick searches of folding pathways. We also show that the choice of the bias potential has significant influence on the efficiency of the bias potential method.     

Contact pair dynamics during folding of two small proteins: chicken villin head piece and the Alzheimer protein beta-amyloid

The folding of an extended protein to its unique native state requires establishment of specific, predetermined, often distant, contacts between amino acid residue pairs. The dynamics of contact pair formation between various hydrophobic residues during folding of two different small proteins, the chicken villin head piece (HP-36) and the Alzheimer protein beta-amyloid (betaA-40), are investigated by Brownian dynamics (BD) simulations. These two proteins represent two very different classes-HP-36 being globular while betaA-40 is nonglobular, stringlike. Hydropathy scale and nonlocal helix propensity of amino acids are used to model the complex interaction potential among the various amino acid residues. The minimalistic model we use here employs a connected backbone chain of atoms of equal size while an amino acid is attached to each backbone atom as an additional atom of differing sizes and interaction parameters, determined by the characteristics of each amino acid. Even for such simple models, we find that the low-energy structures obtained by BD simulations of both the model proteins mimic the native state of the real protein rather well, with a best root-mean-square deviation of 4.5 A for HP-36. For betaA-40 (where a single well-defined structure is not available), the simulated structures resemble the reported ensemble rather well, with the well-known beta-bend correctly reproduced. We introduce and calculate a contact pair distance time correlation function, C(P) (ij)(t), to quantify the dynamical evolution of the pair contact formation between the amino acid residue pairs i and j. The contact pair time correlation function exhibits multistage dynamics, including a two stage fast collapse, followed by a slow (microsecond long) late stage dynamics for several specific pairs. The slow late stage dynamics is in accordance with the findings of Sali et al. Analysis of the individual trajectories shows that the slow decay is due to the attempt of the protein to form energetically more favorable pair contacts to replace the less favorable ones. This late stage contact formation is a highly cooperative process, involving participation of several pairs and thus entropically unfavorable and expected to face a large free energy barrier. This is because any new pair contact formation among hydrophobic pairs will require breaking of several contacts, before the favorable ones can be formed. This aspect of protein folding dynamics is similar to relaxation in glassy liquids, where also alpha relaxation requires highly cooperative process of hopping. The present analysis suggests that waiting time for the necessary pair contact formation may obey the Poissonian distribution. We also study the dynamics of Forster energy transfer during folding between two tagged amino acid pairs. This dynamics can be studied by fluorescence resonance energy transfer (FRET). It is found that suitably placed donor-acceptor pairs can capture the slow dynamics during folding. The dynamics probed by FRET is predicted to be nonexponential.