The free energy of expansion and contraction: treatment of arbitrary systems using the Jarzynski equality

Thermodynamic integration, free energy perturbation, and slow change techniques have long been utilised in the calculation of free energy differences between two states of a system that has undergone some transformation. With the introduction of the Jarzynski equality and the Crooks relation, new approaches are possible. This paper investigates an important phenomenon - systems undergoing a change in volume/density - and derives both the Jarzynski equality and Crooks relation of such systems using a statistical mechanical approach. These results apply to systems with arbitrary particle interactions and densities. The application of this approach to the expansion/compression of particles confined within a vessel with a piston and within a periodic system is considered.     

Ordering of limits in the Jarzynski equality

We consider the sampling problems encountered in computing free-energy differences using Jarzynski's nonequilibrium work relation [Phys. Rev. Lett. 56, 2690 (1997)]. This relation expresses the free-energy change of a system, on which finite-time work is done, as an average over all possible trajectories of the system. This average can then be expressed as a cumulant expansion of the work. We study the scaling of these cumulants with an appropriately defined measure of phase-space accessibility epsilon and particle number N for two simple changes in state. We find that the cumulant expansion is slowly convergent for predominantly entropic processes, those where epsilon is considerably altered during the course of the process. An accurate determination of the free-energy change requires some knowledge of the behavior of the tails of the work distribution associated with the process. Jarzynski's irreversible work relation is only valid with the correct ordering of the infinite limits of N and epsilon, clarifying the regime of its applicability.     

Free energy surfaces from nonequilibrium processes without work measurement

Recent developments in statistical mechanics have allowed the estimation of equilibrium free energies from the statistics of work measurements during processes that drive the system out of equilibrium. Here a different class of processes is considered, wherein the system is prepared and released from a nonequilibrium state, and no external work is involved during its observation. For such "clamp-and-release" processes, a simple strategy for the estimation of equilibrium free energies is offered. The method is illustrated with numerical simulations and analyzed in the context of tethered single-molecule experiments.     

Estimation of protein folding probability from equilibrium simulations

The assumption that similar structures have similar folding probabilities (p(fold)) leads naturally to a procedure to evaluate p(fold) for every snapshot saved along an equilibrium folding-unfolding trajectory of a structured peptide or protein. The procedure utilizes a structurally homogeneous clustering and does not require any additional simulation. It can be used to detect multiple folding pathways as shown for a three-stranded antiparallel beta-sheet peptide investigated by implicit solvent molecular dynamics simulations.     

Basin hopping simulations for all-atom protein folding

We investigate different protocols of the basin hopping technique for de novo protein folding. Using the protein free-energy force field PFF01 we report the reproducible all-atom folding of the 20-amino-acid tryptophan-cage protein [Protein Data Bank (PDB) code: 112y] and of the recently discovered 26-amino-acid potassium channel blocker (PDB code: 1wqc), which exhibits an unusual fold. We find that simulations with increasing cycle length and random starting temperatures perform best in comparison with other parametrizations. The basin hopping technique emerges as a simple but very efficient and robust workhorse for all-atom protein folding.     

Free energy of nascent-chain folding in the translocon

During their synthesis, many water-soluble proteins and nearly all membrane proteins transit through a protein-conducting channel in the membrane, the Sec translocon, from where they are inserted into the lipid bilayer. Increasing evidence indicates that folding of the nascent protein begins already within the ribosomal exit tunnel in a sequence- and environment-dependent fashion. To examine the effects of the translocon on the nascent-chain folding, we have calculated the potential of mean force for α-helix formation of a 10-alanine oligopeptide as a function of its position within the translocon channel. We find that the predominant conformational states, α-helical and extended, reflect those found for the peptide in water. However, the translocon, via its surface properties and its variable diameter, shifts the equilibrium in favor of the α-helical state. Thus, we suggest that the translocon facilitates not only the insertion of membrane proteins into the bilayer but also their folding.     

Accelerated molecular dynamics simulations of protein folding

Folding of four fast‐folding proteins, including chignolin, Trp‐cage, villin headpiece and WW domain, was simulated via accelerated molecular dynamics (aMD). In comparison with hundred‐of‐microsecond timescale conventional molecular dynamics (cMD) simulations performed on the Anton supercomputer, aMD captured complete folding of the four proteins in significantly shorter simulation time. The folded protein conformations were found within 0.2–2.1 Å of the native NMR or X‐ray crystal structures. Free energy profiles calculated through improved reweighting of the aMD simulations using cumulant expansion to the second‐order are in good agreement with those obtained from cMD simulations. This allows us to identify distinct conformational states (e.g., unfolded and intermediate) other than the native structure and the protein folding energy barriers. Detailed analysis of protein secondary structures and local key residue interactions provided important insights into the protein folding pathways. Furthermore, the selections of force fields and aMD simulation parameters are discussed in detail. Our work shows usefulness and accuracy of aMD in studying protein folding, providing basic references in using aMD in future protein‐folding studies. © 2015 Wiley Periodicals, Inc.     

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.     

Estimation of protein folding rate from Monte Carlo simulations and entropy capacity

The problem of protein self-organization is one of the most important problems of molecular biology nowadays. Despite the recent success in the understanding of general principles of protein folding, details of this process are yet to be elucidated. Moreover, the prediction of protein folding rates has its own practical value due to the fact that aggregation directly depends on the rate of protein folding. The time of folding has been calculated for 67 proteins with known experimental data at the point of thermodynamic equilibrium between unfolded and native states using a Monte Carlo model where each residue is considered to be either folded as in the native state or completely disordered. The times of folding for 67 proteins which reach the native state within the limit of 10(8) Monte Carlo steps are in a good correlation with the experimentally measured folding rate at the mid-transition point (the correlation coefficient is -0.82). Theoretical consideration of a capillarity model for the process of protein folding demonstrates that the difference in the folding rate for proteins sharing more spherical and less spherical folds is the result of differences in the conformational entropy due to a larger surface of the boundary between folded and unfolded phases in the transition state for proteins with more spherical fold. The capillarity model allows us to predict the folding rate at the same level of correlation as by Monte Carlo simulations. The calculated model entropy capacity (conformational entropy per residue divided by the average contact energy per residue) for 67 proteins correlates by about 78% with the experimentally measured folding rate at the mid-transition point.     

Free energy calculations using flexible-constrained, hard-constrained and non-constrained molecular dynamics simulations.

A comparison of different treatments of bond-stretching interactions in molecular dynamics simulation is presented. Relative free energies from simulations using rigid bonds maintained with the SHAKE algorithm, using partially rigid bonds maintained with a recently introduced flexible constraints algorithm, and using fully flexible bonds are compared in a multi-configurational thermodynamic integration calculation of changing liquid water into liquid methanol. The formula for the free energy change due to a changing flexible constraint in a flexible constraint simulation is derived. To allow for a more direct comparison between these three methods, three different pairs of models for water and methanol were used: a flexible model (simulated without constraints and with flexible constraints), a rigid model (simulated with standard hard constraints), and an alternative flexible model (simulated with flexible constraints and standard hard constraints) in which the ideal or constrained bond lengths correspond to the average bond lengths obtained from a short simulation of the unconstrained flexible model. The particular treatment of the bonds induces differences of up to 2 % in the liquid densities, whereas (excess) free energy differences of up to 5.7 (4.3) kJ mol(-1) are observed. These values are smaller than the differences observed between the three different pairs of methanol/water models: up to 5 % in density and up to 8.5 kJ mol(-1) in (excess) free energy.     

Free energy of a trans-membrane pore calculated from atomistic molecular dynamics simulations

Atomistic molecular dynamics simulations of a lipid bilayer were performed to calculate the free energy of a trans-membrane pore as a function of its radius. The free energy was calculated as a function of a reaction coordinate using a potential of mean constraint force. The pore radius was then calculated from the reaction coordinate using Monte Carlo particle insertions. The main characteristics of the free energy that comes out of the simulations are a quadratic shape for a radius less than about 0.3 nm, a linear shape for larger radii than this, and a rather abrupt change without local minima or maxima between the two regions. In the outer region, a line tension can be calculated, which is consistent with the experimentally measured values. Further, this line tension can be rationalized and understood in terms of the energetic cost for deforming a part of the lipid bilayer into a hydrophilic pore. The region with small radii can be described and understood in terms of statistical mechanics of density fluctuations. In the region of crossover between a quadratic and linear free energy there was some hysteresis associated with filling and evacuation of the pore with water. The metastable prepore state hypothesized to interpret the experiments was not observed in this region.     

Topology-based models and NMR structures in protein folding simulations

Topology-based interaction potentials are simplified models that use the native contacts in the folded structure of a protein to define an energetically unfrustrated folding funnel. They have been widely used to analyze the folding transition and pathways of different proteins through computer simulations. Obviously, they need a reliable, experimentally determined folded structure to define the model interactions. In structures elucidated through NMR spectroscopy, a complex treatment of the raw experimental data usually provides a series of models, a set of different conformations compatible with the available experimental data. Here, we use an efficient coarse-grained simulation technique to independently consider the contact maps from every different NMR model in a protein whose structure has been resolved by the use of NMR spectroscopy. For lambda-Cro repressor, a homodimeric protein, we have analyzed its folding characteristics with a topology-based model. We have focused on the competition between the folding of the individual chains and their binding to form the final quaternary structure. From 20 different NMR models, we find a predominant three-state folding behavior, in agreement with experimental data on the folding pathway for this protein. Individual NMR models, however, show distinct characteristics, which are analyzed both at the level of the interplay between tertiary/quaternary structure formation and also regarding the thermal stability of the tertiary structure of every individual chain.     

Free Energy Calculations with Non-Equilibrium Methods: Applications of the Jarzynski Relationship

We present an introduction to the Jarzynski relationship that makes a strong connection, for a thermodynamic transformation, between the distribution of non-equilibrium work values and the corresponding equilibrium free energy differences. The relationship is discussed in the context of sampling issues, high level parallel computing and convergence criteria. We discuss three different applications by our group: mechanical unfolding of peptides, mixed quantum/classical free energy calculations in enzymes, and ligand escape pathways     

A novel method for analyzing energy relaxation in condensed phases using nonequilibrium molecular dynamics simulations: application to the energy relaxation of intermolecular motions in liquid water

We present a novel method to investigate energy relaxation processes in condensed phases using nonequilibrium molecular dynamics simulations. This method can reveal details of the time evolution of energy relaxation like two-color third-order IR spectroscopy. Nonetheless, the computational cost of this method is significantly lower than that of third-order response functions. We apply this method to the energy relaxation of intermolecular motions in liquid water. We show that the intermolecular energy relaxation in water is characterized by four energy transfer processes. The structural changes of the liquid associated with the energy relaxation are also analyzed by the nonequilibrium molecular dynamics technique.     

Free energy perturbation simulations of cation binding to valinomycin

Experimental values of the free energies of cation binding to the cyclic depsipeptide molecule, valinomycin, obtained from Pedersen-type salt extraction measurements, provide data against which it is possible to test the adequacy of the procedures and force fields of the molecular dynamics algorithms, MOLARIS and GROMOS. These data are then used to assess appropriate values for the partial charges of the ester carbonyl oxygen and carbon. Valinomycin was chosen because it has only one kind of ion-binding ligand and because the cation is sufficiently enfolded by the molecule in the ion-complexes that the overall size and shape of the complex is virtually the same regardless of the species of cation bound. For such an isosteric complex, the experimentally measured selectivities are sufficiently similar in a wide variety of solvent environments that thedifferences in free energies measured between the different ion-valinomycin complexes by two-phase salt extraction experiments into dichloromethane can be taken as equivalent to the differences in free energies in vacuo. Thesedifferences were therefore compared with those computed for ion-valinomycin complexationin vacuo by Free Energy Perturbation/Molecular Dynamics (FEP/MD) simulations using the MOLARIS and GROMOS programs. Starting with a set of Lennard-Jones 6–12 parameters for the monovalent cations assessed for aqueous solution we explored the effect of varying the partial charges of the ester carbonyl ligands on binding free energy differences (i.e. the selectivity) among Na, K, Rb, and Cs. The computed selectivity was found to depend strongly on the value of partial charge, following a typical Eisenman Selectivity Pattern in which the correct selectivity sequence and magnitude occurred only over a very narrow range of partial charge (around 0.33 and 0.6 for the standard carbonyls of MOLARIS and GROMOS, respectively). Using MOLARIS we explored the effect of varying the size of the ester carbonyl ligands by comparing the standard carbonyl of MOLARIS with the somewhat smaller carbonyls of GROMOS and found an equally satisfactory ability to reproduce the experimental data with a partial charge value of 0.41. These results validate the use of both the MOLARIS and GROMOS force fields as starting points for quantitative calculations of ion-binding in more complex molecules (e.g., ion-binding sites and channels in proteins).This paper is dedicated to the memory of the late Dr C. J. Pedersen.     

All-atom structure prediction and folding simulations of a stable protein

We present results from all-atom, fully unrestrained ab initio folding simulations for a stable protein with nontrivial secondary structure elements and a hydrophobic core. The construct, "trpcage", is a 20-residue sequence optimized by the Andersen group at University of Washington and is currently the smallest protein that displays two-state folding properties. Compared over the well-defined regions of the experimental structure, our prediction has a remarkably low 0.97 A Calpha root-mean-square-deviation (rmsd) and 1.4 A for all heavy atoms. The simulated structure family displays additional features that are suggested by experimental data, yet are not evident in the family of NMR-derived structures.     

Multiple molecule effects on the cooperativity of protein folding transitions in simulations

Though molecular simulation of proteins has made notable contributions to the study of protein folding and kinetics, disagreement between simulation and experiment still exists. One of the criticisms levied against simulation is its failure to reproduce cooperative protein folding transitions. This weakness has been attributed to many factors such as a lack of polarizability and adequate capturing of solvent effects. This work, however, investigates how increasing the number of proteins simulated simultaneously can affect the cooperativity of folding transitions--a topic that has received little attention previously. Two proteins are studied in this work: phage T4 lysozyme (Protein Data Bank (PDB) ID: 7LZM) and phage 434 repressor (PDB ID: 1R69). The results show that increasing the number of proteins molecules simulated simultaneously leads to an increase in the macroscopic cooperativity for transitions that are inherently cooperative on the molecular level but has little effect on the cooperativity of other transitions. Taken as a whole, the results identify one area of consideration to improving simulations of protein folding.