Article Information

• PDF (190 kb)
• Full Text with Large Figures
• Cited by in Scopus (3)

PubMed

• Articles by Bryan A. Patel
• Articles by Peter J. Rossky

Related Articles

• Stability of Triple-Helical Poly(dT)-Poly(dA)-Poly(dT) DNA w...

  Stability of Triple-Helical Poly(dT)-Poly(dA)-Poly(dT) DNA with Counterions Biophysical Journal, Volume 75, Issue 1, 1 July 1998, Pages 70-91 Voichita M. Dadarlat and V.K. Saxena Abstract Structural conformation of triple-helical poly(dT)-poly(dA)-poly(dT) has been a very controversial issue recently. Earlier investigations, based on fiber diffraction data and molecular modeling, indicated an A-form conformation with sugar pucker. On the other hand, Raman, solution infrared spectral, and NMR studies show a B-form structure with sugars. In accordance with these experimental results, a theoretical model with B-form, sugars was proposed in 1993. In the present work we investigate the dynamics and stability of the two conformations within the effective local field approach applied to the normal mode calculations for the system. The presence of counterions was explicitly taken into account. Stable equilibrium positions for the counterions were calculated by analyzing the normal mode dynamics and free energy of the system. The breathing modes of the triple helix are shifted to higher frequencies over those of the double helix by 4–16cm. The characteristic marker band for the B conformation at 835cm is split up into two marker bands at 830 and 835cm. A detailed comparison of the normal modes and the free energies indicates that the B-form structure, with sugar pucker, is more stable than the A-form structure. The normal modes and the corresponding dipole moments are found to be in close agreement with recent spectroscopic findings. Abstract | Full Text | PDF (307 kb)

• A Quantitative Study of λ-Phage SWITCH and Its Components

  A Quantitative Study of λ-Phage SWITCH and Its Components Biophysical Journal, Volume 92, Issue 8, 15 April 2007, Pages 2685-2693 Chunbo Lou, Xiaojing Yang, Xili Liu, Bin He and Qi Ouyang Abstract We propose what we believe is a new model to quantitatively describe the -phage SWITCH system. The model incorporates facilitated transfer mechanism of transcription factor, which can be simplified into a two-step reaction. We first sequentially obtain two indispensable parameters by fitting our model to experimental data of two simple systems, and then apply them to study the natural -SWITCH system. By incorporating the facilitated transfer mechanism, we find that in RecA host , the wild-type lysogenic state is in a monostable regime rather than in a bistable regime. Furthermore, the model explains the weak role of Cro protein and probably sheds light on the evolution of -Cro protein, which is known to be structurally distinct from the other Cros in lambdoid family members. Abstract | Full Text | PDF (651 kb)

• Reaction Coordinates for the Flipping of Genetic Switches

  Reaction Coordinates for the Flipping of Genetic Switches Biophysical Journal, Volume 94, Issue 9, 1 May 2008, Pages 3413-3423 Marco J. Morelli, Sorin Tănase-Nicola, Rosalind J. Allen and Pieter Rein ten Wolde Abstract We present a detailed analysis, based on the forward flux sampling simulation method, of the switching dynamics and stability of two models of genetic toggle switches, consisting of two mutually repressing genes encoding transcription factors (TFs); in one model (the exclusive switch), the two transcription factors mutually exclude each other's binding, while in the other model (general switch), the two TFs can bind simultaneously to the shared operator region. We assess the role of two pairs of reactions that influence the stability of these switches: TF-TF homodimerization and TF-DNA association/dissociation. In both cases, the switch flipping rate increases with the rate of TF dimerization, while it decreases with the rate of TF-operator binding. We factorize the flipping rate into the product of the probability (*) of finding the system at the dividing surface (separatrix) between the two stable states, and a kinetic prefactor . In the case of the exclusive switch, the rate of TF-operator binding affects both (*) and , while the rate of TF dimerization affects only . The general switch displays a higher flipping rate than the exclusive switch, and both TF-operator binding and TF dimerization affect , , and (*). To elucidate this, we analyze the transition state ensemble. For the exclusive switch, the transition state ensemble is strongly affected by the rate of TF-operator binding, but unaffected by varying the rate of TF-TF binding. Thus, varying the rate of TF-operator binding can drastically change the pathway of switching, while changing the rate of dimerization changes the switching rate without altering the mechanism. The switching pathways of the general switch are highly robust to changes in the rate constants of both TF-operator and TF-TF binding, even though these rate constants do affect the flipping rate; this feature is unique for nonequilibrium systems. Abstract | Full Text | PDF (582 kb)

• …more

Copyright © 2007 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 93, Issue 12, 4116-4127, 15 December 2007

doi:10.1529/biophysj.107.108530

Biophysical Theory and Modeling

A Water-Explicit Lattice Model of Heat-, Cold-, and Pressure-Induced Protein Unfolding

Bryan A. Patel^*, ,  , Pablo G. Debenedetti^*, Frank H. Stillinger^† and Peter J. Rossky^‡

^* Department of Chemical Engineering, Princeton University, Princeton, New Jersey
^† Department of Chemistry, Princeton University, Princeton, New Jersey
^‡ Department of Chemistry and Biochemistry, Institute for Theoretical Chemistry, University of Texas at Austin, Austin, Texas

Address reprint requests to Bryan A. Patel, Tel.: 609-258-5413.

Globular proteins remain in their functional native state in a limited range of temperatures and pressures, unfolding into denatured states at both high and low temperature and high pressure 1,2. Fig. 1 shows the characteristic parabolic shape of a protein phase diagram, for the specific case of Staphylococcal nuclease 3. The stability of the native state is the result of several contributions, including van der Waals, electrostatic, and hydrogen-bonding interactions, as well as configurational entropy 4. Walter Kauzmann first suggested that the hydrophobic effect plays a central role in protein stability on the basis of several observations 5. First, proteins tend to sequester the majority of their hydrophobic residues in the core of the molecule, avoiding exposure to water 6. Second, natural water-soluble proteins are unable to fold into their native states in nonpolar solvents 7. Finally, the free energy of transfer of a hydrophobic solute from its pure phase into water shows the same qualitative temperature dependence as the protein unfolding free energy. The liquid-hydrocarbon model of the protein core suggests that the temperature-dependent behavior of the enthalpy and entropy of unfolding is due largely to hydrophobic hydration phenomena 8.

Display large version of this figure

Figure 1

Phase diagram for Staphylococcal nuclease from Fourier transform infra-red (FTIR) spectroscopy, small angle x-ray scattering (SAXS) and differential scanning calorimetry (DSC) experiments. Adapted with permission from Ravindra and Winter 3.

The relevance of this simple picture to high-pressure protein stability is less certain. Kauzmann pointed out inconsistencies between the volume effects associated with protein unfolding and hydrocarbon transfer, thereby questioning whether the liquid-hydrocarbon picture can account for pressure denaturation. The volume change of hydrocarbon transfer is negative at low pressure and generally positive at high pressure, while the volume change of protein unfolding is positive at low pressure and negative at high pressure 9. Calculations from an information theory model of hydrophobic interactions have reconciled these pictures by, in effect, inverting the liquid-hydrocarbon model for pressure effects; pressure denaturation is viewed as the penetration of the protein cavities by water molecules, more akin to the transfer of a water molecule into a pure hydrophobic phase 10,11.

Alternatively, recent theoretical work on the hydration of hydrophobic solutes of different sizes is also of relevance to hydrophobic homopolymer collapse 12,13,14,15. The hydration of large hydrophobic solutes is accompanied by density depletion at their interface with water, which can then resemble the interface between vapor and liquid water 12. Correspondingly, water density fluctuations at the interface would give rise to bubble formation, and bubbles of a critical size would span intermolecular distances and induce the collapse of large hydrophobic polymers 13. From this viewpoint, the stability of this collapsed state is reduced as the solvent moves away from vapor-liquid coexistence at low temperature and high pressure, giving rise to cold and pressure denaturation 12,14. This framework is consistent with the approach adopted here, in that the complexity of real proteins is reduced to an exclusive focus on the thermodynamics of hydrophobicity for the interpretation of the phase diagram of protein folding. A recent molecular dynamics simulation of a hydrophobic polymer in water has manifested the existence of a dewetted polymer-water interface in extended polymer conformations 15. The polymer also exhibited signatures of both warm and cold denaturation, with compact states remaining marginally stable only within an intermediate range of temperature, while extended configurations predominate at high and low temperatures.

Atomically detailed simulations of proteins have provided substantial insight into the mechanisms of protein denaturation and the formation of native-state structure. The choices of force field and whether to model the solvent implicitly or explicitly dictate the types of physical situations that these simulations can tackle. Explicit-solvent molecular dynamics simulations are commonly used to study protein dynamics and hydration for small protein fragments of <40 residues. These studies have provided insight into protein folding pathways and intermediates for several different peptides and proteins, such as the villin headpiece 16,17, a β-hairpin 18, and BBA5 19. Other studies have examined the role of water in hydrophobic collapse of the protein core and protein-water hydrogen bonding for larger proteins like the BphC enzyme 20 and ribonuclease A 21. High-pressure, explicit-solvent simulations utilizing water insertions into the protein interior have investigated the effect of water penetration into the hydrophobic core on pressure denaturation 22.

Modeling the solvent implicitly drastically reduces simulation complexity, allowing larger proteins to be simulated for longer times. These simulations are used to explore the folding dynamics of proteins 23,24,25 and the temperature dependence of protein folding 26,27. Implicit solvent approximations are often used in protein structure prediction, where faster computation is necessary and folding dynamics are less important 28,29.

In contrast to the wealth of existing simulations to study short-time dynamics, few detailed simulations have been done to study protein stability over a wide range of temperatures and pressures, a limitation resulting from the long times required to simulate large-scale conformational changes. One example of such a study is a replica exchange molecular dynamics simulation of a small peptide in TIP3P water that demonstrated both pressure and thermal denaturation, but also showed that pressure did not disrupt the studied α-helical structures 30,31.

The reduced description of interactions used in lattice models of proteins offer the opportunity to examine protein thermodynamics without incurring the large computational expense of atomically detailed simulations. The HP model 32 is a lattice heteropolymer model with hydrophobic (H) and polar (P) monomers. It implicitly incorporates the hydrophobic effect through an attraction between hydrophobic monomers on the protein chain. This approach was used to study the designability and uniqueness of protein native states and to explore how sequence affects the folded structure 32,33. While a broad thermal denaturation transition does occur in the HP model, it does not manifest cold unfolding. Inherently, the model’s ground state is also the protein’s native state, which precludes the existence of a denatured state which is more stable than the native state at lower temperatures. An extension of the HP model 34,35,36,37 to include water explicitly (the HPW model) uses the Muller-Lee-Graziano water model 38,39 to describe the solvent degrees of freedom. The Muller-Lee-Graziano model uses a bimodal description of water, distinguishing between bulk and solvation shell waters in their energy and entropy. From exact enumeration of the configurations of short heteropolymers, Caldarelli and co-workers showed that both a warm and a cold denaturation transition are manifest in the HPW model 34.

Pressure effects have been less commonly treated in theoretical models. A mean-field model for heteropolymer collapse by Dill 40,41 was recently extended by Cheung and co-workers 42 to include pressure-induced weakening of the hydrophobic association of protein monomers. The model exhibits cold-, pressure-, and heat-unfolding but, because it invokes the mean-field approximation, it can only examine the effect of average sequence hydrophobicity on protein thermodynamics. Another model dealing with pressure effects was developed by Marques and co-workers 43. Placing an all-hydrophobic homopolymer and explicit water on a compressible two-dimensional square lattice, they observed cold, thermal, and pressure denaturation. However, their Hamiltonian coupled the compactness of the protein to the bulk water structure, linking a local property of the protein to an average bulk property of water. This feature of the Hamiltonian favors compact protein conformations at conditions where bulk water is highly hydrogen-bonded, in effect forcing the correct outcome a priori, rather than obtaining folding as a result of microscopic interactions. Furthermore, the model predicts that cold-unfolding occurs only at high pressure, a fact not supported by experiments 3,44.

Here we revise the model of Marques et al., and present a different lattice model of an all-hydrophobic protein in water. We distinguish between bulk and interfacial hydrogen bonding, incorporating ideas from Frank’s iceberg model 45. From this simple treatment of protein-water interactions, we can investigate how hydrophobic solvation affects the stability of proteins. A physically meaningful treatment of hydrophobic interactions allows us to reproduce many of the experimentally observed features of protein conformational stability in the pressure-temperature plane, including cold denaturation at ambient pressure.

The outline of this article is as follows. We begin with a description of the microscopic model and discuss the flat histogram simulation techniques used to compute the model’s density of states. We then examine the shape of the resulting protein phase diagrams and discuss the driving forces behind each of the denaturation transitions. Finally, we present the main conclusions, as well as possible extensions to the model, some of which we are currently investigating. The Appendices contain derivations of biased trial moves used in the simulations and analytical approximations used to interpret the results as well as illustrative results for the volumetric properties of small model hydrophobic solutes.

Each site on a two-dimensional square lattice is occupied either by a water molecule or a protein monomer. The water molecules have four hydrogen-bonding arms, each associated with a neighboring lattice site. The orientation of a bonding arm on water molecule i associated with adjacent site j is described by the variable σ_ij. Each bonding arm can adopt q orientations, thus σ_ij can take on values between 1 and q. The orientations of bonding arms on the same water molecule are uncorrelated with each other. A hydrogen bond forms between two neighboring water molecules 〈i, j〉 when their bonding arms are properly oriented, satisfying the condition σ_ij=σ_ji in the original version of the model 46, or a modified condition (to be described below) in this work. To account for the lower local density associated with the formation of hydrogen bonds, we treat the lattice as compressible, and the total volume expands uniformly by an amount Δυ upon formation of a hydrogen bond. The Hamiltonian for the hydrogen-bonding interaction in the original water model 46,47 is

(1)

The protein is modeled as a self-avoiding walk with covalently attached monomers occupying nearest-neighbor sites on the lattice. The protein is a homopolymer and each monomer is nonpolar. The protein has no self-interaction aside from excluded volume effects. The nonpolar protein’s only interaction with water is through its indirect effect on water-water hydrogen-bonding.

We redefine the hydrogen-bonding interaction of the pure water model by differentiating water molecules into two classes: bulk and interfacial. Frank and Evans identified the tendency of water molecules to form ordered “icelike” cages around nonpolar solutes 45. These interfacial water molecules avoid orientations in which their hydrogen-bonding arms point toward the hydrophobe and remain in low entropy configurations in the first solvation shell. However, in these restricted configurations the interfacial waters form stronger hydrogen bonds than observed in bulk 49. While a bimodal description of hydrogen bonding was also used in the HPW model, as already noted, that water model does not reproduce water’s unusual thermodynamics nor does it treat pressure effects 34,35,36,37.

To adapt these principles to our model, we must first recast the orientational criteria for hydrogen bonding to generalize the model. Hydrogen bonding is said to exist when two nearest-neighbor water molecules have their bonding arms oriented such that they are within some tolerance λ, or |σ_ij−σ_ji|≤λ. In the context of this new hydrogen-bonding criterion, a value of λ=0 is equivalent to the original water model’s condition of σ_ij=σ_ji. The range of acceptable orientation pairs differs between bulk water (λ_b) and interfacial water (λ_h). A pair of hydrogen-bonding water molecules is designated as interfacial when either member of the pair is adjacent to one or more protein monomers. A smaller range of hydrogen-bonding orientations for interfacial water molecules (λ_h<λ_b) is the origin in the model of the lower entropy of the interfacial hydrogen bonds in hydrophobic hydration. The tolerances λ_b and λ_h directly affect the fraction of orientation pairs suitable for bonding. The total number of possible orientation pairs is q² for a pair of bonding arms on adjacent water molecules. If one of the bonding arms has any orientation σ_ij out of q possible orientations, then its partner must have an orientation from σ_ij−λ_h to σ_ij+λ_h to satisfy the interfacial hydrogen-bonding criteria. Thus, there are (2λ_h+1)q possible interfacial bonding pair orientations out of q² total orientation pairs. Because λ_h<λ_b, a smaller fraction of orientation pairs will satisfy the bonding criteria interfacially ((2λ_h+1)/q), than in the bulk ((2λ_b+1)/q).

Along with the entropic cost of interfacial hydrogen bonds, there is an enthalpic bonus. Since those hydrogen bonds formed in reduced orientational entropy configurations around hydrophobic solutes less frequently sample the more distorted and weaker bonding structures present in bulk 49,50, they are given an additional enthalpic contribution (J_H) on top of the base hydrogen-bond strength J. Thus, the complete Hamiltonian of our model is

(2)

(3)

Our model is related to that of Marques et al. 43, using the same representation of a protein in water as a hydrophobic homopolymer on a two-dimensional square lattice. However, we use a different Hamiltonian for the protein-water interactions, incorporating the entropic and enthalpic characteristics of water hydrogen bonding on a molecular level into a two-tiered framework. In contrast, as modeled by Marques et al. 43, the protein has a repulsive interaction with the solvent given by

(4)

We used a modified version of the Wang-Landau method to estimate the density of states, Ω 51. In the conventional Wang-Landau approach, the simulation performs a random walk in energy (U) with probability proportional to the reciprocal of the density of states 1/Ω(U). The random walk is performed within the range of attainable energies in the model. The density of states is not known a priori, but is determined on-the-fly during the simulation. The simulation is initialized at a random configuration with the density-of-states estimator set at Ω(U)=1 for all energy levels. Trial moves from an old configuration (o) to a new configuration (n) at energy levels U_o and U_n, respectively, are accepted with probability

(5)

In our simulations, we employ a slightly different strategy for determining the density of states that is adapted specifically to our model 43. Instead of a random walk in energy space, we perform a random walk in the two variables in our system Hamiltonian: N_HB and N_XHB. Any state specified by these variables corresponds to a specific system energy and volume given by Eqs. (2). The outcome of these simulations is the density of states, Ω(N_HB, N_XHB), which can then be translated to Ω(U, V). The advantage of performing a random walk in N_HB and N_XHB, as opposed to U and V, is that the parameters J, J_H, and Δυ need not be assigned values during the simulation. This allows us to gather the system thermodynamics for any parameter set from a single simulation.

Here we consider the histogram of visited states to be flat if every bin of h(N_HB, N_XHB) is at least 80% of the average histogram 〈h(N_HB, N_XHB)〉. The simulation ends when the modification factor is <exp(10⁻⁸). To improve sampling of rare protein configurations we also used several biased protein trial moves in addition to conventional translational and orientational trial moves. An explanation and derivation of the acceptance criteria for these trial moves is given in Appendix A .

Due to the large number of system states sampled even for a relatively small protein (∼85,000 for 20 monomers), the density of states was subdivided into smaller overlapping regions for expediency. Even with a small subset of phase space to sample, the system states around and including the protein native state required that the simulation be run in parallel on multiple processors. The processors independently performed random walks in the same subset of phase space and periodically communicated their estimates for Ω and h with each other, simultaneously testing whether the criteria for histogram flatness had been satisfied. The complete simulation of a 20-monomer protein required nearly 6000 CPU hours. The simulation of proteins up to 50 monomers is possible with a recently developed method which decouples the protein and water contributions to the density of states. This technique will be explained in detail in a forthcoming publication.

Extracting protein properties from the simulation data requires converting the density of states into more useful metrics. Since the simulation involves fluctuations in both internal energy and volume, it is natural to reweight the density of states to represent the isobaric-isothermal ensemble. Given a pressure P and temperature T, the probability of a state, i, specified by N_HB and N_XHB, is

(6)

(7)

(8)

Fig. 2 shows the phase diagram of a 17-mer protein, with contours showing constant probabilities of occupying the folded state. The model protein exhibits heat-, cold-, and pressure-unfolding and the dome of stable temperatures and pressures has the same qualitative shape as observed experimentally. The folded protein is stabilized by minimizing the exposed surface area of hydrophobic monomers, thereby limiting the number of interfacial water molecules forced to pay an entropic cost for hydrogen-bonding around the protein. This corresponds to the maximally compact conformation shown in Figure 3a. Upon increasing temperature, the protein gradually unfolds and exposes more hydrophobic monomers to the solvent. The thermal energy is sufficient to overcome the entropic cost of forming additional interfacial hydrogen bonds around the exposed monomers. The thermally denatured protein is an ensemble of conformations, with several examples shown in Figure 3b. At lower temperatures and high pressures the enthalpic benefit of additional interfacial hydrogen bonds drives the protein to unfold into a fully extended state, shown in Figure 3c. This configuration maximizes the number of water-protein contacts and allows for water to form the maximum number of interfacial hydrogen bonds.

Display large version of this figure

Figure 2

The phase diagram of a 17-mer homopolymer protein for J_H/J=0.2, Δυ/υ₀=0.348, λ_h=0, λ_b=1, and q=30. The inner line demarcates the region within which the probability of observing the folded state is 87.5% or greater. The other lines similarly demarcate the regions within which the folded probabilities are >75%, 62.5%, 50% (bold), 37.5%, 25%, and 12.5% (outermost).

Display large version of this figure

Figure 3

Representative conformations for a 17-mer model protein in the (a) native state, (b) thermally denatured ensemble of states, and (c) high-pressure cold-unfolded state.

The importance of the entropic penalty and the enthalpic bonus in cold and pressure denaturation is illustrated in Figure 4 and Figure 5. When λ_b=0, there is no entropic penalty (since λ_b=λ_h=0 in this case) and one sees that the protein is not susceptible to cold denaturation. The entropic penalty was applied by increasing λ_b while holding λ_h=0. Increasing the entropic penalty increases the slope of the cold denaturation transition line as well as the maximum pressure at which the protein is stable at any nonzero temperature. This effect will be explained quantitatively below. Fig. 5 shows that there is no cold unfolding without the enthalpic bonus. Increasing the enthalpic bonus increases the stability of the cold-denatured state at the expense of the folded state. At J_H/J=0.3, the protein shows cold denaturation even at zero pressure, while for J_H/J>0.3 the protein no longer folds at positive pressures. The effect of changing J_H demonstrates that the formation of stronger hydrogen bonds between the water molecules close to the protein is also a required driving force for cold-unfolding in this model.

Display large version of this figure

Figure 4

Contours of 50% folded probability for a 20-mer protein for varying values of the entropic penalty for forming interfacial hydrogen bonds. Simulations were performed with λ_h=0 and changing λ_b. To retain the same bulk water thermodynamics, the total number of water orientations q increases such that the fraction of bonding configurations for a pair of bonding arms on adjacent molecules (i.e.: (2λ_b+1)q/q²=(2λ_b+1)/q) is kept constant, at 0.1. The other model parameters remain constant at J_H/J=0.2 and Δυ/υ₀=0.348.

Display large version of this figure

Figure 5

Contours of 50% folded probability for a 20-mer protein for varying values of the enthalpic bonus for interfacial hydrogen bonds J_H/J. The other model parameters are Δυ/υ₀=0.348, λ_b=1, λ_h=0, and q=30.

An analytical expression derived in Appendix B yields further insight into the mechanism of cold denaturation in this model. It is shown there that the transition pressure between the native and unfolded states at zero temperature can be calculated analytically, and is given by

(9)

The effect of the entropic cost on cold denaturation can be determined from the slope of the cold-denaturation line given by Eq. (10), also derived in Appendix B ,

(10)

In contrast to cold unfolding, which is driven both by entropy and enthalpy effects, thermal denaturation at higher temperatures is an entropy-driven transition. The partial opening of the protein upon thermal denaturation breaks a few bulk hydrogen bonds (ΔN_HB<0), increasing both the protein and water disorder (ΔS>0). The simulations also show that the number of interfacial hydrogen bonds remains generally constant upon thermal denaturation (ΔN_XHB≈0). The expression for the thermal denaturation temperature at zero pressure, derived in Appendix B , can be simplified to

(11)

Display large version of this figure

Figure 6

Thermal denaturation temperature at zero pressure for proteins of various sizes. J_H/J=0.2, Δυ/υ₀=0.348, λ_h=0, λ_b=1, and q=30 were used for all proteins. The line is a guide to the eye. Irregularities in the trend of decreasing thermal denaturation temperature with protein size are due to the discrete nature of the lattice. Equation (11) connects the structure of the native and denatured states to the thermal denaturation temperature through the change in entropy upon unfolding. The degeneracy of the native, compact state does not increase regularly with protein size. For example, a 9-mer or 16-mer has a square-shaped native state, while other protein sizes have rougher surfaces.

Fig. 7 examines the entropy and volume changes upon unfolding in our model and in a phenomenological model based on experimental studies of the most typical proteins 1. While cold unfolding and high pressure thermal unfolding in our model show the same denaturation thermodynamics as experiments, it is clear that the model does not reproduce the typical phenomenological slope of the phase diagram for thermal denaturation at ambient (low) pressure. Experimentally, there typically is a point on the denaturation curve where the slope of the phase transition for thermal denaturation is nearly infinite, corresponding to a volume change upon unfolding ΔV ≈ 0. This point is lacking in our model, along with regions of positive volume change upon unfolding. The cause originates in the model, where Eq. (3) directly links the volume of the system to the number of hydrogen bonds. For our model, denaturation is always associated with a negative volume change because protein unfolding involves the breaking of hydrogen bonds. A discussion of the volumetric properties of small hydrophobic solutes based on our model and their relation to the volumetric behavior associated with protein folding is given in Appendix C . There it is shown that, for small solutes, the model also correctly mimics most of the trends with temperature and pressure seen experimentally. Nevertheless, our model does show a point on the denaturation curve where ΔS ≈ 0, allowing for a distinct separation between cold denaturation (ΔS<0) and thermal denaturation (ΔS>0).

Display large version of this figure

Figure 7

A comparison between the thermodynamics of unfolding in (a) our model and in (b) a schematic summary of experimental results. Panel b was adapted with permission from Hawley 1.

To relate the model results to experimental protein phase diagrams, we converted the temperature and pressure into dimensional quantities. A value of 23kJ/mol was used in the temperature scaling for the average strength of a hydrogen bond, J. To scale the pressure we used 18cm³/mol for the molar volume of water. The value of Δυ/υ₀ was chosen to approximate experimental values of the volume expansion of a water molecule upon forming a hydrogen bond. Using the densities of ice and liquid water of 0.917 and 1g/cm³ at 0°C, the molar volume change upon freezing is 1.63cm³/mol of water. The ratio of the heat of fusion to the heat of sublimation indicates that ∼13% of hydrogen bonds are broken upon ice melting. An individual water molecule in ice has two hydrogen bonds per molecule, therefore, on average, 0.26 hydrogen bonds per molecule are broken upon melting. Dividing the molar volume change upon freezing (1.63cm³/mol) by the fraction of hydrogen bonds per molecule formed upon freezing (0.26) gives an estimate for the volume change upon forming a mole of hydrogen bonds, or Δυ=6.27cm³/mol.

For a parameter selection of J_H/J=0.4, λ_h=0, λ_b=5, q=110, and Δυ/υ₀=0.348, the phase diagram of a model 20-mer is shown in Fig. 8. Pressure denaturation of the model protein is in the kbar region, of the same order of magnitude as the experimental results shown in Fig. 1. The temperatures for thermal and cold denaturation at ambient pressure are between 20 and 40°C below those seen experimentally. The model underestimates these temperatures because of the lack of favorable internal protein-protein interactions such as hydrogen-bonding and electrostatics. Including these forces might confer additional thermal stability to the protein.

Display large version of this figure

Figure 8

Contour of 50% folded probability for a model 20-mer with temperature and pressure converted into dimensional quantities using J=23kJ/mol and υ₀=18cm³/mol. Parameter values are J_H/J=0.4 and λ_b=5, λ_h=0, q=110, and Δυ/υ₀=0.348.

Finally, a comparison between our model predictions and those of water-implicit protein models yields an interesting observation. The thermal denaturation of proteins in the original HP model is a broad transition, occurring over a range of dimensionless temperatures of (1) 32. Our model clearly exhibits a much sharper transition between the native and denatured states, as shown in Fig. 2. The explicit inclusion of water gives rise to sharper phase transitions, in agreement with experimental observations. It remains to be seen whether sharp transitions are an artifact of the homopolymer case, or if heteropolymers would also exhibit this phenomenon.

We have developed a model for an all-hydrophobic homopolymer in water based on the thermodynamics of the solvation of hydrophobic solutes. The model protein denatures at high temperature, low temperature, and high pressure, showing many of the same denaturation characteristics observed in experiments. Our model shows sharp unfolding transitions compared to water-implicit protein models. The key new feature here is the introduction of additional restrictions on the orientation for hydrogen bonding at the protein interface compared to the bulk. The two-tiered set of hydrogen-bonding interactions for the water modulates the cold denaturation transition by balancing between the stability of the native and cold-denatured states.

Using a homopolymer to describe the protein limits our study to the effects of protein size and of the water-protein interaction. The model also shows a decrease in the protein configurational degrees of freedom upon cold unfolding, an unphysical feature that we plan to improve upon. Calculations are in progress on a heteropolymer model, with both hydrophobic and polar monomers, as in the HP model 32. Incorporating details of the chemistry of proteins into the model will allow us to investigate how sequence selection and composition can affect the range of native state stability. We also plan to investigate three-dimensional systems, as well as the impact of correlation among the orientations of hydrogen bonds on each water molecule.

Appendix A

Biased protein trial moves

In addition to simple one-monomer and two-monomer translation trial moves, we used biased protein monomer translation moves to improve the sampling of rare protein configurations. In these biased trial moves, one or two monomers of the protein are translated without disrupting the hydrogen bonding structure of the solvent around it. To keep the number of hydrogen bonds constant, we bias the selection of water orientations for the water molecule which translates into the site formerly occupied by the protein monomer. This type of bias allows the protein to effectively explore configurations when the solvent is highly structured. To balance the improved sampling of biased trial moves with the shorter computation time of unbiased trial moves, 50% of protein translation trial moves are unbiased and 50% are biased.

The transition matrix that determines the probability to perform a one-monomer trial move from o to n is

(12)

where p_i is the probability of selecting monomer at site i for translation, p_j is the probability of selecting site j for the monomer to translate into, p_w is the probability of translating water w at site j back to site i, and p_σ is the probability of giving water molecule w the set of bonding arm orientations {σ} ({σ}={σ_i1, σ_i2, σ_i3, σ_i4}). In the biased trial move, all of these probabilities are symmetric, except p_σ, as we bias the selection of these orientations to keep the number of hydrogen bonds formed by water w (N_HB) constant from state o to state n.

Water w has a total of four hydrogen-bonding arms that are faced with several options in the new site depending on their neighbors. If a bonding arm points toward the protein it can have one of q orientations since the water molecule cannot hydrogen-bond with the protein. If a bonding arm points toward other water molecules the arms can either hydrogen-bond by matching their orientation to their neighbor or form nonbonding pairs. The total number of bonding arm pairs irrespective of their bonding status is denoted by N_pairs. To maintain a constant value of N_HB for water w, N_HB of the N_pairs bonding arm pairs are selected at random to be hydrogen-bonded and the orientations set accordingly. The remaining N_pairs−N_HB nonbonding pairs are given nonbonding orientations at random. If N_HB>N_pairs, the move is rejected automatically.

Based on these criteria, the probability of selecting orientations for water w is given by

(13)

The first term is the probability of selecting a set of N_HB bonding pairs out of the N_pairs total bonding arms for hydrogen-bonding. The second term is the probability of selecting appropriate hydrogen-bonding orientations for those N_HB bonding arms. The third term is the probability of selecting appropriate nonbonding orientations for the remaining pairs. The final term is the probability of selecting any orientation for the bonding arms facing protein monomers.

Detailed balance requires that the opposing fluxes between states o and n be statistically equal, satisfying the relation

(14)

where p is the equilibrium probability of a state and A is the acceptance probability of a transition. Since the equilibrium probability of a state is the inverse of its density of states, we can derive the acceptance criteria:

(15)

The principle is the same for a biased two-monomer trial move, with some additions. With two translating water molecules there are now eight bonding arms that are given biased orientations. Moreover, the translating waters can be subject to either bulk or interfacial hydrogen bonding criteria. The probability for selecting a set of orientations is

(16)

where the subscripts i and b denote interfacial and bulk, respectively. The acceptance criteria for a biased two-monomer move is

(17)

where the change in a quantity such as ΔN_pairs is N_pairs(n)−N_pairs(o).

Appendix B

Analytical approximation for two-state transitions

The probability of any state, i, specified by N_HB and N_XHB can be calculated from simulation data for any pressure and temperature by

(18)

where Δ is the isothermal-isobaric partition function and k_B is Boltzmann’s constant; Ω_i is the degeneracy of state i, while U_i and V_i are the energy and volume of state i, calculated from Eqs. (2). At a two-state transition between states 1 and 2, p₁=p₂. Equating p₁ and p₂ and taking the natural logarithm gives the relationship

(19)

where P_t and T_t are the transition pressure and temperature. Using Boltzmann’s relation S_i=k_B ln Ω_i to replace the density of states with the entropy of state i and multiplying by k_B reduces the equation to

(20)

Substituting for the energy and volume using Eqs. (2)and rearranging yields the relation

(21)

ΔN_HB, ΔN_XHB, and ΔS denote the difference in the respective property between states 1 and 2. Recall that Δυ is not the change in a property between states 1 and 2, but instead the model parameter describing the uniform expansion of the lattice upon hydrogen-bond formation.

Equation (21) explicitly links the system states specified by N_HB and N_XHB, the model parameters J, J_H, and Δυ, and the macroscopic variables T_t and P_t. An adaptation of the Clapeyron equation for this model can be derived by solving for P_t in Eq. (21) and taking the derivative with respect to temperature:

(22)

Equations (21) allow us to examine the shape of the protein phase diagram in greater detail and to investigate more deeply which mechanisms contribute to cold-, pressure-, and heat-denaturation. Two limiting cases are used to study the effects of model parameters on the thermodynamics. We derive the pressure at which a transition occurs between two ground states by solving for P_t in Eq. (21) when T_t=0:

(23)

This expression gives the transition pressure between two ground states: the cold-denatured state stable at high pressures and the native state stable at low pressures.

We can also use the transition temperature at zero pressure to study thermal denaturation:

(24)

This analysis is an approximation valid for a two-state transition without significantly populated intermediate states. The simulation results presented in this article show that this approximation is essentially exact during cold denaturation at low temperatures because both the native and cold-denatured states are well-defined and are the only highly populated states. Because thermal denaturation involves a transition into an ensemble of denatured states, the analysis is a less exact description at higher temperatures.

Appendix C

Partial molar volume calculations

As noted by Kauzmann 9 and mentioned in the Introduction, the volume change of transfer of hydrophobic solutes into water and the volume change of protein unfolding typically exhibit different dependences on pressure. Simulations of small hydrophobic solutes in water were performed to probe the model’s ability to describe volumetric properties and to determine if the discrepancies pointed out by Kauzmann are observed. The solutes examined were smaller versions of the model homopolymer: one-site, two-site, and flexible three-site hydrophobic solutes in a fully occupied lattice of water. Changing the solute concentration by replacing water molecules with solute molecules on the lattice induces changes in the system volume at any given temperature and pressure. From the slope of this trend, we can calculate the quantity (∂V/∂N₁)_T,P,N, where N₁ is the number of solute molecules, N₂ is the number of solvent molecules, and the N subscript denotes a fixed number of occupied lattice sites. Note that this quantity is not the partial molar volume of the solute, but a relationship between the two can be derived by examining the total differential for the volume,

(25)

where the solute and solvent partial molar volumes ( and respectively) appear as

(26)

The solute’s partial molar volume can be related to the quantity estimated from simulation results by taking the partial derivative of Eq. (25) with respect to N₁ while fixing the temperature, pressure, and total number of lattice sites, which yields the relation

(27)

The quantity (∂N₂/∂N₁)_T,P,N is −1 for a one-site solute, −2 for a two-site solute, and so forth. At infinite dilution Eq. (27) reduces to

(28)

for the case of a one-site solute, where is the partial molar volume of the solute at infinite dilution. υ₂ is the molar volume of the solvent, which we can calculate from pure water simulations.

The temperature and pressure behavior of the partial molar volume of our model solutes in water correspond qualitatively to experimental trends for small apolar solutes. Figure 9ace, show that the partial molar volumes of each of the model solutes increases with increasing temperature. This trend matches the experimental observations of Masterton, who observed increasing partial molar volumes with increasing temperature for methane, ethane, propane, and benzene 53. The calculations also show that partial molar volumes decrease with increasing pressure (Figure 9bdf), conforming to the trends observed experimentally for methane 54, benzene, and toluene 55.

Display large version of this figure

Figure 9

Partial molar volumes of the model hydrophobic solutes in water, in units of υ₀. Data from graphs a, c, and e are at constant P=0. Data from graphs b, d, and f are at constant T=0.108. Parameter values are J_H/J=0.2, λ_b=3, λ_h=0, q=70, and Δυ/υ₀=0.348.

As stated in the Introduction, the measured volume change upon transfer from a hydrophobic phase into water typically, although not always, exhibits a different pressure dependence than the volume change upon protein unfolding. The volume change upon hydrophobic transfer for a model solute can be estimated as the difference between the partial molar volume of the solute in water and the molar volume of the pure solute,

(29)

The model definition for the molar volume of the pure solute is N_sitesυ₀, where N_sites is simply the number of sites comprising the solute (i.e.: 1, 2, or 3 in our case). In units of υ₀, it reduces simply to the number of monomers in the solute, i.e., 1 for a one-site solute, 2 for a two-site solute, etc. We note that υ₀ is a temperature- and pressure-independent parameter, and hence our model for a hydrophobic homopolymer in water is not designed to predict the thermodynamic properties of a pure hydrophobic species.

Given this definition, we observe negative values of Δυ_transfer for each solute over the ranges of temperature and pressure shown in Fig. 9. The partial molar volume of the solutes under these conditions remains less than the molar volume of the pure solute. This behavior at low pressures matches the observed experimental transfer volume changes, where the transfer of methane from liquid hexane into water at 1 bar corresponds to a molar volume decrease of 22.7cm³/mol 5. However, at high pressures where small hydrophobic solute transfer volumes are positive 9, this model predicts negative transfer volumes.

The behavior of Δv_transfer for model solutes is consistent with the volume change upon unfolding for the model protein. As shown schematically in Figure 7a, the model protein exhibits a volume decrease upon unfolding at low and high pressures. While this does not conform to Kauzmann’s assertions about the pressure dependence of protein volume changes 9, nor to experimental observations of many proteins such as chymotrypsinogen 1, some proteins such as staphylococcal nuclease 3 show a volume decrease upon thermal denaturation at low pressures. The factors that determine whether a protein has a positive or negative volume change upon unfolding are still a matter of some debate 56. The model protein presented here captures the correct volumetric properties of a subset of real proteins.