Logo of biophysjGuide for AuthorsAbout this journalExplore this journalBiophysical Journal
Biophys J. 2013 Sep 17; 105(6): 1482–1490.
PMCID: PMC3785879
PMID: 24048000

Solvation Free Energy of the Peptide Group: Its Model Dependence and Implications for the Additive-Transfer Free-Energy Model of Protein Stability

Dheeraj S. Tomar,†▵ D. Asthagiri,†▵∗ and Valéry Weber‡
Author information Article notes Copyright and License information Disclaimer

Associated Data

Supplementary Materials

Abstract

The group-additive decomposition of the unfolding free energy of a protein in an osmolyte solution relative to that in water poses a fundamental paradox: whereas the decomposition describes the experimental results rather well, theory suggests that a group-additive decomposition of free energies is, in general, not valid. In a step toward resolving this paradox, here we study the peptide-group transfer free energy. We calculate the vacuum-to-solvent (solvation) free energies of (Gly)n and cyclic diglycine (cGG) and analyze the data according to experimental protocol. The solvation free energies of (Gly)n are linear in n, suggesting group additivity. However, the slope interpreted as the free energy of a peptide unit differs from that for cGG scaled by a factor of half, emphasizing the context dependence of solvation. However, the water-to-osmolyte transfer free energies of the peptide unit are relatively independent of the peptide model, as observed experimentally. To understand these observations, a way to assess the contribution to the solvation free energy of solvent-mediated correlation between distinct groups is developed. We show that linearity of solvation free energy with n is a consequence of uniformity of the correlation contributions, with apparent group-additive behavior in the water-to-osmolyte transfer arising due to their cancellation. Implications for inferring molecular mechanisms of solvent effects on protein stability on the basis of the group-additive transfer model are suggested.

Introduction

The thermodynamics of protein unfolding in the presence of aqueous osmolytes (small organic cosolutes) is of fundamental interest in understanding the forces stabilizing the folded protein (for example, see (1–9)) and in the broader quest to understand how biological systems adapt to environmental stresses by changing the solvent properties (10–13). Experimental investigations on the molecular role of the osmolytes on the thermodynamics of protein unfolding are usually anchored by

ΔGNU=ΔGNU(0)+∑iαiΔgtr,i,
(1)

where ΔGNU and ΔGNU(0) are the experimentally accessible unfolding free energies of the protein in the osmolyte solution and in water, respectively. (For a 1 M osmolyte solution, ΔGNU−ΔGNU(0) is the m-value (2); here, we follow the sign convention of Auton and Bolen (14)). Drawing upon the seminal studies by Tanford and co-workers (1,15), ΔGNU−ΔGNU(0) is partitioned into the water-to-osmolyte-solution transfer free energy, Δgtr,i, contribution due to group i, potentially allowing one to understand the role of the solvent at the level of each individual group. Typically, the groups include the peptide unit and the side chains of the amino acid residues, and the Δgtr,i contributions are obtained from chemically apposite model compounds. The factor αi accounts for the fractional change in the degree of solvent exposure of group i in going from the native (N) to the unfolded (U) state.

The pioneering experimental studies by Bolen and co-workers interpreted using Eq. 1 leads to an important insight: both conformation-protecting and denaturing osmolytes exert their influence primarily by changing the solubility of the peptide backbone (3,7), an identification with significant consequences to our understanding of protein folding (16). However, these results pose a paradox: whereas the group-additive decomposition in Eq. 1 is, in general, not valid (to be discussed below), ΔGNU−ΔGNU(0) calculated using Eq. 1, together with reasonable approximations of the unfolded state of the protein, describes the experimental ΔGNU−ΔGNU(0) rather well (17,18). Toward resolving this paradox, here we examine the vacuum-to-solvent and water-to-osmolyte-solution transfer free energies of the peptide group and its dependence on the choice of model compounds.

Glycyl peptides, (Gly)n, are a common model for the peptide unit. The peptides can be blocked n-acetyl-glyn-methyl amides (as in this study) or zwitterionic (15). Using blocked glycyl peptides, and by carefully minimizing peptide-peptide interactions in peptide solubility, Auton and Bolen (14) have sought the peptide-group transfer free energy, Δgtr,p, that is “independent of the model compound and the choice of concentration scale.” In their studies, Δgtr,p was obtained by appropriately combining the transfer free energy of chains of various lengths n, and, somewhat more robustly, by equating Δgtr,p to the slope of the transfer free energy with respect to n, the so-called constant-increment method. Model independence was demonstrated by showing that the transfer free energies based on (Gly)n agree with those from cGG, the cyclic-diglycine molecule. This concordance, although pleasing, is puzzling as well, for in cGG, the CO and NH of the peptide are cis and the molecule has a net zero dipole moment, whereas in the (Gly)n system, the CO and NH are trans, and the peptide dipole moment is nonnegligible. Thus, either the conformation of the peptide is unimportant in Δgtr (Eq. 1) or there are other factors that lead to this result, or a combination of both.

Implicit in Eq. 1 is the important assumption that the same unfolded-state conformation is present in both water and the aqueous osmolyte solution. In practice, limiting models of the unfolded state—an extended protein (more rigorously a random coil in a good solvent) on one end, and, on the other end, a structure obtained by combining excised fragments of folded proteins that retain intramolecular interactions representative of compact denatured states—are considered and an attempt is made to bracket the m-value (17,18). However, conformational preferences can be different in different phases, due to both the intramolecular interaction in a given conformation and how that conformation is solvated in the medium, as can be directly seen within the multistate picture of hydration (19,20). (Among the manifestations of this effect is the phenomenon of self-solvation (21,22), an effect that reveals limitations of additivity.) It is of course remarkable that despite such obvious approximations, the transfer free energy model (Eq. 1) appears to describe the experiments as well as it does (17,18). Here, paralleling the experimentalist’s approach, we consider only peptides in the fully extended conformation and dissect how and why apparently additive behavior emerges for this scenario that is likely favorable for additivity.

Important past simulation examinations of group additivity in peptide hydration appear to come to opposite conclusions. Using a continuum dielectric model of the solvent, Avbelj and Baldwin (23,24) argue against group additivity. For example, they find that in an extended (Ala)9 chain, the electrostatic contribution to the hydration free energy of the peptide unit, as well as its interaction with other peptide units, depends on the location of the chosen peptide unit. However, from all-atom simulations of (Gly)n, Hu et al. (25) find that both the vacuum-to-water and the water-to-osmolyte-solution transfer free energies are reasonably linear with n. Since the condition for the use of the constant-increment approach seems to be satisfied, these authors concluded in favor of group additivity and obtained the group contribution from the slope of the free energy versus n.

However, as we show below, the seemingly contradictory conclusions from the earlier simulation studies can be reconciled; and this reconciliation depends on appreciating that an independent group-additive contribution to the hydration free energy is not a consequence even when the conditions for use of the constant-increment approach are satisfied. Here we use all-atom simulations and the quasichemical organization of the potential distribution theorem (19,20,26) to examine the solvation free energies of acetyl-(Gly)n-methyl amide peptides and of cGG. A virtue of the quasichemical formulation is that it makes transparent the role of correlated fluctuations of the binding energies of two groups on the molecule and its role in the thermodynamics of hydration, issues that are of central interest in understanding group additivity.

We show that even for an idealized solute that is incapable of making any close (near-neighbor) contacts with the solvent, the group-solvent binding energies between neighboring groups are correlated and make a nonnegligible contribution to the net solvation free energy. This implies that identifying a model independent, group contribution to the solvation free energy is, in principle, not possible, even for this idealized solute. The situation for a real solute is bound to be considerably more complicated.

Group additive behavior in the solvation free energy can be suggested if solvent-mediated correlation between different groups are similar for a series of model compounds, such as (Gly)n, but the identified group-additive free energy will necessarily depend on the chosen model. The existence of correlations is consistent with the continuum dielectric study (23,24) and its uniformity is what leads to a linear dependence of solvation free energy on n, consistent with the earlier all-atom study (25). It is important to note that such correlation effects arise at different energy scales: within the quasichemical formulation, it is found for the idealized solute that is incapable of near-neighbor contacts with the solvent, a collection of spheres that only repels the solvent, and for the physical solute. However, group additivity in the water-to-osmolyte-solution transfer free energy can arise if these correlation contributions cancel, but in that case, care is needed in inferring mechanistic conclusions about solvent effects on stability from such an apparently group-additive model.

Methods

The simulation procedure closely follows our earlier work (27). The peptides are modeled in the extended configuration with the long axis aligned with the diagonal of the simulation cell and the center of the peptide placed at the center of the simulation cell. (Initial configurations were energy-minimized with restraints to keep the peptide extended.) The peptide atoms are held fixed throughout the simulation. The solvent was modeled by the TIP3P (28,29) model and the CHARMM (30) force field with correction terms for dihedral angles (31) was used for the peptide. A total of 2006 TIP3P molecules solvated the peptide. Parameters for urea and trimethylamine N-oxide (TMAO) were obtained from Weerasinghe and Smith (32) and Kast et al. (33), respectively. A total of 449 urea molecules (for a molar concentration of ∼8 M) and 195 TMAO molecules (for a molar concentration of ∼4 M) were used.

The simulations were performed at a temperature of 298.15 K and a pressure of 1 bar, as described earlier (27,34). The equations of motion was propagated using the Verlet algorithm with a time step of 2.0 fs. SHAKE (35) was used to constrain the geometry of water molecules. The Lennard-Jones interactions were terminated at 10.43 Å by smoothly switching to zero starting at 9.43 Å. Electrostatic interactions were treated with the particle mesh Ewald method with a grid spacing of 0.5 Å.

We applied atom-centered fields to carve out a molecular cavity in the liquid, in contrast to the spherical cavity used in our earlier study (27). The functional form of the field is given by Eq. 4b in that study (27). To build the field to its eventual range of λ = 5 Å, we progressively applied it in unit Å increments and computed the work done in applying the field using Gauss-Legendre quadrature. In Section S.II of the Supporting Material, we present further details of this procedure and how the conditional binding energy distribution, P(ε|ϕλ), was calculated.

cGG was built and optimized using the Gaussian (G09) quantum chemistry package (36). For consistency with the (Gly)n simulations, the partial charges and Lennard-Jones interaction parameters were obtained from the backbone atoms of the CHARMM forcefield.

Theory

The excess chemical potential, μex, of a solute in the solution is that part of the Gibbs free energy that would vanish if the interaction between the solute and solvent were to vanish. Formally (20,37),

βμex = ln〈eβε〉, 
(2)

where 〈…〉 denotes ensemble averaging with respect to P(ε), the probability density distribution of the binding (or interaction) energy, ε, of the solute with the solvent. As usual, β = 1/kBT, where T is the temperature and kB the Boltzmann constant.

Following earlier work (27,34), to calculate μex from Eq. 2, we regularize P(ε) by introducing an auxiliary constraint, a field ϕλ that pushes the solvent molecules away from the solute’s surface to a range λ. (The envelope defined by λ also specifies the inner shell.) This construct has the virtue of tempering the solute-solvent interaction, and, for solvent pushed far enough (typically evacuating the first hydration shell is sufficient), the distribution of binding energies is Gaussian. Formally, with the introduction of the field,

βμex=lnx0[ϕλ]︸localchemistry−lnp0[ϕλ]︸packing+βμex[P(ε|ϕλ)]︸long−range.
(3)

−ln x0[ϕλ] is the free energy required to apply the field in the solute-solvent system: it reflects the strength of the solute interaction with the solvent in the inner shell. −ln p[ϕλ] is the free energy required to apply the field in the neat solvent system: it reflects the intrinsic properties of the solvent. For ϕλ modeling a hard exclusion of solvent −ln p0[ϕλ], is precisely the hydrophobic contribution to hydration (38). βμex[P(ε|ϕλ)] is the contribution to βμex from long-range solute-solvent interactions. In molecular dynamics simulations, we calculate −ln x0[ϕλ] or −ln p0[ϕλ] simply by the work required to apply ϕλ. Fig. 1 presents a schematic of the decomposition of μex according to Eq. 3.

An external file that holds a picture, illustration, etc. Object name is gr1.jpg

Schematic depicting the quasichemical organization (Eq. 3) of Eq. 2.

The chemistry and packing contributions defined in Eq. 3 also have an alternate physical interpretation that proves helpful in appreciating the correlation effects to be defined below.

For example, −lnx0[ϕλ]=ln(1+∑i≥1Ki[ϕλ]ρwi), where i water molecules (from the bulk at a density ρw) are sequestered within the inner shell and the equilibrium constant for forming the solute plus i water cluster is Ki (19,20,26,39). A similar clustering expansion, with an equilibrium constant K˜i, obtained in the absence of the solute, applies to the packing contribution.

For ϕλ excluding solvent from the first hydration (or inner shell), the conditional binding energy distribution, P(ε|ϕλ), can be well described by a Gaussian of mean 〈ε|ϕλ〉 and variance 〈δε2|ϕλ〉 (27,34,40), and the long-range contribution is then given by

μex[P(ε|ϕλ)]=〈ε|ϕλ〉+β2〈δε2|ϕλ〉.
(4)

Correlation effects in the long-range contribution

Now, consider decomposing ε into contributions εi due to various groups i = 1, …, n, comprising the solute under consideration. For a pairwise additive force field, as is used in this study, such a decomposition can be unambiguously made. For the conditioned solute, even the individual binding energy distributions, P(εi|ϕλ), are Gaussian distributed, but in general, εi is correlated with εj(j ≠ i). When two conditioned solutes are spatially adjacent, the presence of one solute can influence the binding energy distribution of the other; physically this happens because the solvent has to jointly satisfy binding with both the solutes. In this correlated case, P(ε = ∑iεi) is Gaussian distributed with a mean ∑i〈εi|ϕλ〉 and a variance ∑i〈δεi2|ϕλ〉+2∑i>jρij〈δεi2|ϕλ〉〈δεj2|ϕλ〉, where ρij is the correlation coefficient (41). The long-range contribution to the free energy is then given by

μex[P(ε|ϕλ)]=∑iμex[P(εi|ϕλ)]+β∑i>jρij〈δεi2|ϕλ〉〈δεj2|ϕλ〉,
(5)

where μex[P(εi|ϕλ)] is described by Eq. 4. The second summation in Eq. 5 can be rewritten as a sum over all nearest-neighbor pairs, ∑(i,i+1), the next-nearest pairs, ∑(i,i+2), etc. From the summation arranged in this fashion, we can then identify the effect of correlations at various spatial length scales to the free energy μex[P(ε|ϕλ)]. Note that this formulation precisely identifies the contributions solely due to the individual groups, namely the quantities μex[P(εi|ϕλ)]; we call this the self-contribution of the group i. For ease of presentation, when we speak of, say, an (i, i + 2) correlation, we mean the correlation between the binding energies of groups i and i + 2, respectively, and the solvent.

Correlation effects in packing and chemistry

To identify correlation contributions in ln p0, we consider a group i and the rest of the peptide forming the background, iback. (For clarity, we do not display the dependence on the field, ϕλ.) Then, the probability that the inner shell is evacuated around i and iback is, by Bayes’ rule, p0(i · iback) = p0(i) × p0(iback|i) = p0(i)p0(iback)pc, where pc is defined as the correlation contribution. (Of course, p0(i · iback) is the quantity p0 in Eq. 3.) p0(i) and p0(iback) are obtained using the same procedure as used for p0. pc is then evaluated from these quantities.

In contrast to Eq. 5, pc does not admit a simple algebraic expression. Nevertheless, the physical meaning of pc can be appreciated as follows. Consider two initially well-separated spherical cavities, assumed identical for simplicity. Further, for simplicity, we assume that only j water molecules populate the cavity with an equilibrium constant K˜j; without loss of generality, j can be the dominant coordination state of the cavity (39). When these two j water clusters are now brought together, we can expect interference between the j water clusters. For example, when the two cavities just touch each other, the equilibrium constant for forming the 2j water cluster, in the absence of any interference, is just K˜j2. However, physically we expect the equilibrium constant for distributing 2j water molecules in the combined volume of the two cavities to be different from the noninterfering case. The correlation contribution is just a measure of the deviation from this noninterfering case. Of course, the interference effects will be greatly amplified when the cavity volumes overlap. Thus, within the clustering perspective, the additivity approximation amounts to assuming no cluster interference.

A decomposition similar to the one noted above also applies to x0, with the physical meaning of xc paralleling the physical meaning of pc. (As a technical aside, note that when calculating x0(i), the rest of the chain remains physically coupled with the solvent and vice versa for x0(iback).)

Results and Discussion

Solvation of the conditioned solute

We first consider the solvation of the conditioned solute, as this helps us focus on solvation effects uncluttered by the complexities of near-neighbor solute-solvent interactions.

Hydration of (Gly)n

The solvation free energy of the conditioned solute, μnex[P(ε|ϕλ)], depends linearly on n for all solutions considered here. For the analysis below, we exclusively focus on the vacuum-to-water transfer (Fig. 2).

An external file that holds a picture, illustration, etc. Object name is gr2.jpg

Long-range contribution to the free energy of blocked (Gly)n modeled in the extended conformation in water. Open circles are the simulation results. Solid circles are based on using the average values of the self (i) and (i, i + 1) correlation contributions from the (Gly)7 chain (Table 1) to reconstruct the free energy for all other n. (See text for further details.)

For understanding the correlation contributions to the slope, we focus on the internal groups of the (Gly)n chain, as these are the ones changing with n. For the blocked (Gly)7 model, group 0 is the methyl, group 1 is the CONHCH2 group formed between the acetyl group and the N-terminus of the protein, and group 8 is the terminal CONHCH3 group. The remaining six (6) CONHCH2 groups are termed the internal groups. For the (Gly)3 model, per this convention, there are two internal groups. In Table 1, we present the average contribution due to various orders of correlation between these internal groups.

Table 1

Average values of correlation contributions of various orders per internal peptide unit

Water
Urea
TMAO
(Gly)3(Gly)7(Gly)7(Gly)7
(i)−1.58−1.53−1.97−1.53
(i, i + 1)−0.78−0.79−0.79−0.79
(i, i + 2)—0.250.250.25
(i, i + 3)—−0.11−0.11−0.11
Total−2.36−2.18−2.62−2.18

(i) indicates that only the contribution of the group with the solvent is included (the first term on the right in Eq. 5); (i, i + 1) indicates the first-neighbor correlation contribution, βρi,i+1〈δεi2|ϕλ〉〈δεi+12|ϕλ〉, and so on. For reference, note that the slope of the μnex[P(ε|ϕλ)] versus n curve for water (Fig. 2) is −2.23 kcal/mol. All values are in kcal/mol.

Notice that the self-contribution ((i), Table 1) is fairly different from the slope of the μnex[P(ε|ϕλ)] versus n curve. As Eq. 5 shows, this term is also the one that can be rigorously identified as a contribution solely due to the group. Progressively including contributions from (i, i + 1), (i, i + 2), etc., correlations, we find a sum that is reasonably close to the slope of μnex[P(ε|ϕλ)] versus n. (The slight discrepancy between the sum computed in Table 1 and the slope arises because the linear fit is not perfect.) Observe that the contribution from various orders of correlation are fairly similar for (Gly)7 and (Gly)3. Likewise, the correlation of the end groups with the internal groups are also fairly similar for these two models (data not shown), implying insensitivity to chain length for the correlations involving long-range interactions.

The identification of the self and correlation contributions depends on how (Gly)n is partitioned, whereas μnex[P(ε|ϕλ)] is itself invariant to the chosen partitioning scheme. If instead of decomposing the chain into single peptide units i, adjacent peptide units are treated as one group, we find that as i increases, the sum of the self-contributions itself progressively approaches μnex[P(ε|ϕλ)] (Section S.V in the Supporting Material). This intuitively reasonable result emphasizes the length-scale dependence of (approximate) additivity and the subtleties involved in identifying a group for a group contribution scheme.

It proves insightful to consider how well the average values of various orders of correlation for (Gly)7 describe the free energy for all other chain lengths. To this end, we take the average value of the (i, i + 1) correlation contribution from Table 1, the self-contributions for groups 0, 1, and 8, and the (i, i + 1) contribution between groups 0 and 1 and between groups 7 and 8, and use Eq. 5 to compute the free energy for all n. These values are referred to as reconstructed in Fig. 2. (The correlation coefficients, ρij, for (Gly)7 are given in Table S2.) The good agreement for all n, including n = 1 (which is all end-groups in our notation), reveals the underlying uniformity of these self-, (i), and nearest-neighbor, (i, i + 1), correlations in this model system.

From Table 1, we find that ignoring all higher-order correlations will lead to an error in the net free energy of ∼0.65 × 7 ≈ 4.6 kcal/mol; ignoring only the (i, i + 3) correlations still leads to an error of ∼0.7 kcal/mol. To put these numbers in perspective, the m-value, which is the principal target of the group additive decomposition (Eq. 1), is usually only a few kcal/mol (17) for modest-sized proteins. Borrowing ideas from Dill (42), for the group-additive decomposition of m to be satisfactory, a rough estimate for the allowable error/group for a modest ∼100-amino-acid protein is ∼10 cal/mol/group. (The experimental error for the peptide group transfer contribution is roughly of this magnitude (14).) As the simulation results suggest, ignoring long-range correlation contributions itself can easily lead to much larger errors.

Long-range correlations in aqueous osmolyte solutions

Table 1 compares the average values of the various orders of correlation in the solvation of (Gly)7 in different solvents. Remarkably, we found that all orders of correlation excluding the self-contribution are identical. Physically this reflects the electrostatic basis of the long-range correlations, and these, being sensitive to the dielectric constant (already high for water itself), do not change much across osmolyte solutions considered here. Thus, for the idealized solute, the transfer free energy from water to the osmolyte solution can be entirely determined by the self-contribution, which is also the contribution that obeys group additivity.

Solvation of the physical solute

Fig. 3 shows that μnex versus n for blocked (Gly)n obeys a linear dependence. Similar linearity also holds for the chemistry, packing, and long-range contributions (Eq. 3) individually. (The reasons for the linearity of the chemistry and packing contributions are considered below.) According to the constant-increment method (14), we consider the slope of these curves as the contribution of an individual group to the free energy. These values are collected in Table 2 together with results for cGG. (Following established experimental approach (14), the value for cGG is scaled by half to obtain the value for one CH2CONH group.)

An external file that holds a picture, illustration, etc. Object name is gr3.jpg

The solvation free energy (Eq. 3) versus n for blocked (Gly)n in the extended conformation. TMAO and urea are aqueous solutions. Straight-line fits are shown by a solid line for water and a dashed line for urea. The fit for TMAO is indistinguishable from that for water within statistical uncertainties and hence is not shown.

Table 2

Peptide-group transfer free energies from vacuum to solvent obtained from the slope of μnex versus n

WaterUreaTMAO
(Gly)n−5.0−5.4−5.0
−0.40.0
cGG/2−6.2−6.6−6.2
−0.40.0

Values for cGG have been scaled by 1/2. Below each line for the model system studied, we present the free- energy values for transferring from water to the solution under study. At the specified osmolyte concentration, using the linear extrapolation approximation (2,25), the experimental estimate of the transfer free energies for urea (8 M) and TMAO (4 M) are−0.34 kcal/mol and 0.35 kcal/mol, respectively. All values are in kcal/mol. Standard error (1 σ) for water-to-osmolyte transfer is ∼0.1 kcal/mol (see Supporting Material).

As found experimentally (14), the simulated water-to-osmolyte transfer free energies are found to be indifferent to the peptide model used. For water-to-aqueous-urea transfer free energy, the agreement between experiments and simulations is also satisfactory: the urea concentration is ∼8 M, and assuming a linear dependence of transfer free energy on osmolyte concentration (2,25), we find that for 1 M urea solution, the transfer free energy is −50 ± 13 cal/mol versus the experimental estimate of −43 ± 10 cal/mol (14). We find a net-zero transfer free energy to aqueous TMAO solution (4 M), in contrast to the experimental estimate of 87 cal/mol (14), a discrepancy that is likely due to inadequacy in the force-field model for TMAO (25,42).

From Table 2, we note that the good agreement between (Gly)n and cGG in water-to-aqueous-osmolyte transfer free energy masks the rather poor agreement in transfer free energies from vacuum to the respective solution. Although it can be argued that water-to-osmolyte-solution transfer is the most relevant transfer experimentally, from the perspective of a physical theory, the vacuum-to-solution transfer quantities have the virtue of highlighting the role of intergroup correlations transparently.

To further illuminate the model dependence of the transfer free energies (Table 2), we consider how these values are used in modeling the m-value. As Eq. 1 indicates, the transfer free energies are scaled by the fractional solvent exposure of the group relative to that in the model compound (17,18). Using commonly used atomic radii (43–45) for calculating solvent-accessible surface area (SASA), we find that the SASA of the CONHCH2 group in cGG is ∼1.5 times that in the (Gly)n model. (The SASA of the peptide in cGG is obtained by dividing the SASA for the entire molecule by 2. For the (Gly)n model, the change in SASA with n gives the SASA/peptide.) Thus, relative to (Gly)n, the transfer free energy/unit area of the model compound is ∼33% smaller in magnitude for the peptide from cGG. In the context of the m-value, using the peptide from cGG as a model can lead to both quantitative and qualitative errors.

Correlation effects in the solvation of the physical solute

Packing contributions

Table 3 lists the correlation contribution to p0 for two different chain lengths and two different solvents. (Recall that in calculating p0, solute-solvent interactions are turned off.) In each case, we focus on the peptide group from the middle of the chain.

Table 3

Correlation contributions in the packing contribution

Solventnp0p0 (i)p0 (iback)pc
Water462.022.760.5−21.2
681.122.978.8−20.6
Urea465.924.564.2−22.8
685.724.683.2−22.1

p0, p0(i), and p0(iback) values are reported in energy units by taking the logarithm and multiplying by −kBT. For each peptide of length n, the group i refers to the central peptide in the chain. For example, for n = 6 (there being seven peptide groups), i = 4. For water, the standard error on pc is 0.3 kcal/mol and for urea it is 0.4 kcal/mol. pc averaged over n = 4 and n = 6 cases for water is −20.9 ± 0.2, and for urea it is −22.5 ± 0.3. All values are in kcal/mol.

Consider the n = 6 case for water. We find that the work done in evacuating the inner shell that will accommodate i alone (22.9 kcal/mol) plus the work for evacuating the inner shell to accommodate the background, iback, alone (78.8 kcal/mol) is ∼21 kcal/mol more than would be required if these domains were evacuated simultaneously. This is physically reasonable, because evacuating the shell for i should aid in evacuating the shell around the neighbors: that is, water molecules populating the inner shell around the peptide interact cooperatively.

Comparing n = 4 and n = 6 cases for water, we find that within statistical uncertainties the correlation contributions are similar. This uniformity of correlation explains, as it did for the long-range contributions, why a linear dependence of the packing contribution versus n is obtained. Also note that within the partitioning of local and background followed here, the change in free energy is entirely due to the variation in the length of the background chain. Finally, the trends noted for water also hold for urea. However, the correlation contributions differ between the solvents.

The analysis of the correlation effects in the packing contribution shows that even for a collection of cavities, nonadditive, solvent-specific behavior should be expected.

Chemistry contributions

Table 4 lists the correlation contribution to x0 corresponding to the cases considered in Table 3. Notice also that the net local contribution, x0(i), due to a defined group is invariant to the chain length, as it should be, provided the end effects are small. Further, as noted above, the change in free energy is entirely due to the variation in the length of the background chain. Just as we found for the correlation effects in packing, the correlation contributions are uniform across chain lengths for a given solvent but differ between solvents.

Table 4

Correlation contributions in the chemistry contribution

solventnx0x0 (i)x0 (iback)xc
Water4−74.4−19.4−67.912.9
6−98.3−19.7−91.713.1
Urea4−78.9−21.1−72.114.3
6−104.2−21.1−97.714.6

x0, x0(i), and x0(iback) values are reported in energy units by taking the logarithm and multiplying by +kBT. For water, the standard error on xc is 0.3 kcal/mol and for urea it is 0.4 kcal/mol. The average xc is 13.0 ± 0.2 for water and 14.5 ± 0.3 for urea. All values are in kcal/mol.

Cancellation of correlation effects in water-to-osmolyte-solution transfer

Table 5 combines the results for the correlation contribution to the net chemistry-plus-packing contribution. Within the statistical uncertainty of the calculation, the correlation contribution to the net short-range effects is the same for both urea and water, just as we found for the long-range contributions (Table 1). Thus, in the water-to-aqueous-osmolyte transfer these correlation contributions will cancel, leaving only a group-additive contribution.

Table 5

Correlation contributions to the net chemistry plus packing contribution

Solventnx0p0x0(i)p0(i)x0(iback)p0(iback)xcpc
Water4−12.43.3−7.4−8.3
6−17.23.2−12.9−7.5
Urea4−13.03.4−7.9−8.5
6−18.53.5−14.5−7.5

The x0/p0, x0(i)/p0(i), and x0(iback)/p0(iback) values are reported in energy units as in Table 4. The standard error on kBT ln xc/pc is ∼0.6 kcal/mol. The average xc/pc is −7.9 ± 0.4 for water and −8.0 ± 0.4 for urea. Within the noted uncertainties, the correlation contribution is the same for both urea and water. All values are in kcal/mol.

The uniformity of correlation effects seen in the long-range (Table 1 and Fig. 2), packing (Table 3), and chemistry (Table 4) contributions explains why the net μnex is linear in n (Fig. 3) for the extended (Gly)n chains considered here. However, as the above analysis makes clear, linearity of μnex versus n does not mean we can define a group-additive contribution ignoring the context in which the group is found. The above analysis also helps reconcile the linear μnex versus n observed in an earlier all-atom simulation (25) that was interpreted as supporting group additivity with the continuum dielectric calculations that questioned additivity (23,24).

Concluding Discussions

Group additivity is a widely used approximation in chemistry. For a small-molecular solute in the gas phase, a system that is entirely characterized by strong, short-range interactions, group additivity may be an acceptable approximation. However, in the treatment of a many-body system, such as a protein in a solvent, characterized by many different scales of energies, group-additive ideas must be considered with sufficient care.

Group-additive models have been often used in describing protein folding, and cautionary notes based on scaling of errors with chain length have been presented as well (46). Efforts to explore this issue by all-atom simulation are beginning to appear. Notably, the failure of additivity was demonstrated in a recent exhaustive study of the hydration of 15 amino acids and their analogs (47). As in that earlier work, here we have used all-atom simulations to examine the meaning of a group-additive transfer free energy of a peptide unit. Distinguishingly, here we have developed a framework to understand how the solvation thermodynamics of a defined group depends on its neighbors. Such a framework, it is hoped, will help improve our understanding of both the successes and limitations of additivity and may perhaps show the way to more effective models.

Based on the aforementioned framework, we find that even for an idealized solute with no near-neighbor solute-solvent interaction, the net solvation free energy of the solute comprises contributions due to the correlated interaction of the solvent with distinct groups in the solute. As is intuitively reasonable, the contribution of individual group-solvent interaction to the net free energy is the most dominant. However, the binding energy of a group i with the solvent is correlated with the binding energy of its neighbors i ± 1, i ± 2, etc. These correlated fluctuations can either raise or lower the free energy of the solute, and, it should be noted, can persist even for spatially distant groups. For the linear (Gly)7, correlations persist up to (i, i + 3). Consideration of the length-scale dependence of additivity shows that identifying separate free energy contributions due to polar and nonpolar groups, a commonly used construct in protein biophysics, is itself of limited validity, especially if those groups are adjacent spatially.

For the physical solute, we find that water molecules that populate the inner hydration shell behave in a cooperative manner, emphasizing that it is not possible to treat a defined group as independent of its neighbors. A similar behavior, perhaps even longer-range of correlations, can be expected for a topographically complicated object and for objects with formal charges, such as charged amino acid residues. On this basis, it is difficult to rationalize the common practice of subtracting the transfer free energy of glycine from that for the target amino acid to isolate a contribution due to the side chain alone.

For the extended-(Gly)n chains studied here, the correlation contributions to the short-range and long-range contributions to the free energy are similar in both water and the aqueous osmolyte solutions studied here. Thus, in the water-to-aqueous-osmolyte transfer, these correlation contributions largely cancel, resulting in an apparently group-additive behavior. We suspect that cancellation of correlations may partly explain the success of the group-additive molecular transfer model in describing the experimental m-value (17,18). Inclusion of conformational flexibility, not considered here, will likely weaken the solvent-mediated inter-residue correlations, further aiding the group-additive decomposition of net water-to-aqueous-osmolyte transfer free energy. On the other hand, differential weighting of conformations in different solvents and presence of charged residues will all render group additivity less tenable. Rigorously evaluating such cases in realistic proteins is necessarily left for future studies. The analysis here does indicate, especially given that m-values are typically only a few kcal/mol, that care may be needed in drawing mechanistic conclusions about molecular forces stabilizing a protein based on a group-additive decomposition of the m-value. Indeed, the case of urea-induced protein denaturation appears to support this cautionary note: whereas the group-additive approach suggests that urea denatures a protein largely because of favorable interactions with the peptide backbone (48), simulation suggests that denaturation is driven by promiscuous urea-protein interactions mediated by (largely nonspecific) dispersion interactions (49).

Acknowledgments

We thank Robert Baldwin, George Rose, Michael Paulaitis, Safir Merchant, and Purushottam Dixit for helpful discussions. We thank the two anonymous reviewers for their constructive comments and criticisms.

This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE- AC02-05CH11231.

Supporting Material

Document S1. Figs. S1–S3, Tables S1 and S2, and Supporting Methods:
Document S2. Article plus Supporting Material:

References

1. Tanford C. Protein denaturation. C. Theoretical models for the mechanism of denaturation. Adv. Protein Chem. 1970;24:1–95. [PubMed] [Google Scholar]
2. Greene R.F., Jr., Pace C.N. Urea and guanidine hydrochloride denaturation of ribonuclease, lysozyme, α-chymotrypsin, and β-lactoglobulin. J. Biol. Chem. 1974;249:5388–5393. [PubMed] [Google Scholar]
3. Bolen D.W., Baskakov I.V. The osmophobic effect: natural selection of a thermodynamic force in protein folding. J. Mol. Biol. 2001;310:955–963. [PubMed] [Google Scholar]
4. Athawale M.V., Dordick J.S., Garde S. Osmolyte trimethylamine-N-oxide does not affect the strength of hydrophobic interactions: origin of osmolyte compatibility. Biophys. J. 2005;89:858–866. [PMC free article] [PubMed] [Google Scholar]
5. Beck D.A.C., Bennion B.J., Daggett V. Simulations of macromolecules in protective and denaturing osmolytes: properties of mixed solvent systems and their effect on water and protein structure and dynamics. In: Häussinger D., Sies H., editors. Methods in Enzymology: Osmosensing and Osmosignalling. Elsevier; London: 2007. pp. 373–396. [PubMed] [Google Scholar]
6. O’Brien E.P., Dima R.I., Thirumalai D. Interactions between hydrophobic and ionic solutes in aqueous guanidinium chloride and urea solutions: lessons for protein denaturation mechanism. J. Am. Chem. Soc. 2007;129:7346–7353. [PubMed] [Google Scholar]
7. Bolen D.W., Rose G.D. Structure and energetics of the hydrogen-bonded backbone in protein folding. Annu. Rev. Biochem. 2008;77:339–362. [PubMed] [Google Scholar]
8. Kokubo H., Hu C.Y., Pettitt B.M. Peptide conformational preferences in osmolyte solutions: transfer free energies of decaalanine. J. Am. Chem. Soc. 2011;133:1849–1858. [PMC free article] [PubMed] [Google Scholar]
9. Guinn E.J., Pegram L.M., Record M.T., Jr. Quantifying why urea is a protein denaturant, whereas glycine betaine is a protein stabilizer. Proc. Natl. Acad. Sci. USA. 2011;108:16932–16937. [PMC free article] [PubMed] [Google Scholar]
10. Yancey P.H., Clark M.E., Somero G.N. Living with water stress: evolution of osmolyte systems. Science. 1982;217:1214–1222. [PubMed] [Google Scholar]
11. Hochachka P.W., Somero G.N. Oxford University Press; New York: 2002. Biochemical Adaptation. Mechanism and Process in Physiological Evolution; pp. 217–289. [Google Scholar]
12. Yancey P.H. Organic osmolytes as compatible, metabolic and counteracting cytoprotectants in high osmolarity and other stresses. J. Exp. Biol. 2005;208:2819–2830. [PubMed] [Google Scholar]
13. Bandyopadhyay A., Saxena K., Chakraborty K. Chemical chaperones assist intracellular folding to buffer mutational variations. Nat. Chem. Biol. 2012;8:238–245. [PMC free article] [PubMed] [Google Scholar]
14. Auton M., Bolen D.W. Additive transfer free energies of the peptide backbone unit that are independent of the model compound and the choice of concentration scale. Biochemistry. 2004;43:1329–1342. [PubMed] [Google Scholar]
15. Nozaki Y., Tanford C. The solubility of amino acids and related compounds in aqueous urea solutions. J. Biol. Chem. 1963;238:4074–4081. [PubMed] [Google Scholar]
16. Rose G.D., Fleming P.J., Maritan A. A backbone-based theory of protein folding. Proc. Natl. Acad. Sci. USA. 2006;103:16623–16633. [PMC free article] [PubMed] [Google Scholar]
17. Auton M., Rösgen J., Bolen D.W. Osmolyte effects on protein stability and solubility: a balancing act between backbone and side-chains. Biophys. Chem. 2011;159:90–99. [PMC free article] [PubMed] [Google Scholar]
18. Auton M., Bolen D.W. Predicting the energetics of osmolyte-induced protein folding/unfolding. Proc. Natl. Acad. Sci. USA. 2005;102:15065–15068. [PMC free article] [PubMed] [Google Scholar]
19. Pratt L.R., Asthagiri D. Free Energy Calculations: Theory and Applications in Chemistry and Biology. Springer; Berlin: 2007. Potential distribution methods and free energy models of molecular solutions; pp. 323–351. [Google Scholar]
20. Beck T.L., Paulaitis M.E., Pratt L.R. Cambridge University Press; Cambridge, United Kingdom: 2006. The Potential Distribution Theorem and Models of Molecular Solutions. [Google Scholar]
21. Roseman M.A. Hydrophilicity of polar amino acid side-chains is markedly reduced by flanking peptide bonds. J. Mol. Biol. 1988;200:513–522. [PubMed] [Google Scholar]
22. König G., Boresch S. Hydration free energies of amino acids: why side chain analog data are not enough. J. Phys. Chem. B. 2009;113:8967–8974. [PubMed] [Google Scholar]
23. Avbelj F., Baldwin R.L. Limited validity of group additivity for the folding energetics of the peptide group. Proteins. 2006;63:283–289. [PubMed] [Google Scholar]
24. Avbelj F., Baldwin R.L. Origin of the change in solvation enthalpy of the peptide group when neighboring peptide groups are added. Proc. Natl. Acad. Sci. USA. 2009;106:3137–3141. [PMC free article] [PubMed] [Google Scholar]
25. Hu C.Y., Kokubo H., Pettitt B.M. Backbone additivity in the transfer model of protein solvation. Protein Sci. 2010;19:1011–1022. [PMC free article] [PubMed] [Google Scholar]
26. Paulaitis M.E., Pratt L.R. Hydration theory for molecular biophysics. Adv. Protein Chem. 2002;62:283–310. [PubMed] [Google Scholar]
27. Weber V., Asthagiri D. Regularizing binding energy distributions and the hydration free energy of protein cytochrome c from all-atom simulations. J. Chem. Theory Comput. 2012;8:3409–3415. [PubMed] [Google Scholar]
28. Jorgensen W., Chandrasekhar J., Klein M.L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983;79:926–935. [Google Scholar]
29. Neria E., Fischer S., Karplus M. Simulation of activation free energies in molecular systems. J. Chem. Phys. 1996;105:1902–1921. [Google Scholar]
30. MacKerell A.D., Jr., Bashford D., Karplus M. All-atom empirical potential for molecular modeling and dynamics studies of proteins. J. Phys. Chem. B. 1998;102:3586–3616. [PubMed] [Google Scholar]
31. Mackerell A.D., Jr., Feig M., Brooks C.L., 3rd Extending the treatment of backbone energetics in protein force fields: limitations of gas-phase quantum mechanics in reproducing protein conformational distributions in molecular dynamics simulations. J. Comput. Chem. 2004;25:1400–1415. [PubMed] [Google Scholar]
32. Weerasinghe S., Smith P.E. A Kirkwood-Buff derived force field for mixtures of urea and water. J. Phys. Chem. B. 2003;107:3891–3898. [Google Scholar]
33. Kast K.M., Brickmann J., Berry R.S. Binary phases of aliphatic N-oxides and water: force field development and molecular dynamics simulation. J. Phys. Chem. A. 2003;107:5342–5351. [Google Scholar]
34. Weber V., Merchant S., Asthagiri D. Communication: regularizing binding energy distributions and thermodynamics of hydration: theory and application to water modeled with classical and ab initio simulations. J. Chem. Phys. 2011;135:181101. [PubMed] [Google Scholar]
35. Ryckaert J.P., Ciccotti G., Berendsen H.J.C. Numerical integration of the Cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 1977;23:327–341. [Google Scholar]
36. Frisch M.J., Trucks G.W., Fox D.J. Gaussian, Inc.; Wallingford CT: 2009. Gaussian 09, Revision A.1. [Google Scholar]
37. Widom B. Potential-distribution theory and the statistical mechanics of fluids. J. Phys. Chem. 1982;86:869–872. [Google Scholar]
38. Pratt L.R., Pohorille A. Theory of hydrophobicity: transient cavities in molecular liquids. Proc. Natl. Acad. Sci. USA. 1992;89:2995–2999. [PMC free article] [PubMed] [Google Scholar]
39. Merchant S., Asthagiri D. Thermodynamically dominant hydration structures of aqueous ions. J. Chem. Phys. 2009;130:195102. [PubMed] [Google Scholar]
40. Asthagiri D., Ashbaugh H.S., Pratt L.R. Non-van der Waals treatment of the hydrophobic solubilities of CF4. J. Am. Chem. Soc. 2007;129:10133–10140. [PubMed] [Google Scholar]
41. Utiramerur S., Paulaitis M.E. Cooperative hydrophobic/hydrophilic interactions in the hydration of dimethyl ether. J. Chem. Phys. 2010;132:155102. [PubMed] [Google Scholar]
42. Dill K.A. Additivity principles in biochemistry. J. Biol. Chem. 1997;272:701–704. [PubMed] [Google Scholar]
43. Canchi D.R., Jayasimha P., Garcia A.E. Molecular mechanism for the preferential exclusion of TMAO from protein surfaces. J. Phys. Chem. B. 2012;116:12095–12104. [PMC free article] [PubMed] [Google Scholar]
44. Bondi A. van der Waals volumes and radii. J. Phys. Chem. 1964;68:441–451. [Google Scholar]
45. Lee B., Richards F.M. The interpretation of protein structures: estimation of static accessibility. J. Mol. Biol. 1971;55:379–400. [PubMed] [Google Scholar]
46. Sanner M.F., Olson A.J., Spehner J.C. Reduced surface: an efficient way to compute molecular surfaces. Biopolymers. 1996;38:305–320. [PubMed] [Google Scholar]
47. König G., Bruckner S., Boresch S. Absolute hydration free energies of blocked amino acids: implications for protein solvation and stability. Biophys. J. 2013;104:453–462. [PMC free article] [PubMed] [Google Scholar]
48. Auton M., Holthauzen L.M.F., Bolen D.W. Anatomy of energetic changes accompanying urea-induced protein denaturation. Proc. Natl. Acad. Sci. USA. 2007;104:15317–15322. [PMC free article] [PubMed] [Google Scholar]
49. Hua L., Zhou R., Berne B.J. Urea denaturation by stronger dispersion interactions with proteins than water implies a 2-stage unfolding. Proc. Natl. Acad. Sci. USA. 2008;105:16928–16933. [PMC free article] [PubMed] [Google Scholar]
Articles from Biophysical Journal are provided here courtesy of The Biophysical Society