Thermodynamic Molecular Switch in Sequence-Specific Hydrophobic Interaction: Two Computational Models Compared

We have shown in our published work the existence of a thermodynamic switch in biological systems wherein a change of sign in ΔCp°(T)reaction leads to a true negative minimum in the Gibbs free energy change of reaction, and hence, a maximum in the related Keq. We have examined 35 pair-wise, sequence-specific hydrophobic interactions over the temperature range of 273–333 K, based on data reported by Nemethy and Scheraga in 1962. A closer look at a single example, the pair-wise hydrophobic interaction of leucine-isoleucine, will demonstrate the significant differences when the data are analyzed using the Nemethy-Scheraga model or treated by the Planck-Benzinger methodology which we have developed. The change in inherent chemical bond energy at 0 K, ΔH°(T0) is 7.53 kcal mol compared with 2.4 kcal mol, while ‹ts› is 365 K as compared with 355 K, for the Nemethy-Scheraga and Planck-Benzinger model, respectively. At ‹tm›, the thermal agitation energy is about five times greater than ΔH°(T0) in the Planck-Benzinger model, that is 465 K compared to 497 K in the Nemethy-Scheraga model. The results imply that the negative Gibbs free energy minimum at a well-defined ‹ts›, where TΔS° = 0 at about 355 K, has its origin in the sequence-specific hydrophobic interactions, which are highly dependent on details of molecular structure. The Nemethy-Scheraga model shows no evidence of the thermodynamic molecular switch that we have found to be a universal feature of biological interactions. The Planck-Benzinger method is the best known for evaluating the innate temperature-invariant enthalpy, ΔH°(T0), and provides for better understanding of the heat of reaction for biological molecules.


INTRODUCTION
In 1971, T.H. Benzinger [19] proposed a thermal work function to take into account both Boltzmann statistical energy effects and energies of quantum-mechanical bonds. While the latter are usually not altered significantly in macromolecular reactions, it was Benzinger's conjecture that the large-scale and long-range changes of conformation that accompany protein folding or assembly might generate significant energy differences due to the cumulative alteration of their numerous covalent bond structures.
Our 1988 [9] and 1994 [10] analyses of available and highly precise data on protein folding in several ribonuclease systems, as well as six self-associating protein systems over the temperature range 220-360 K, enabled us to show that the "chemical bond" component of the overall enthalpy was of significant magnitude, and that it was independent of temperature. This analysis demonstrated the validity of the central concept.
Our subsequent studies [1,2,3,4,5,6,7,9,10,11,12,13,14,15,16] have shown unequivocally that each macromolecular process has a temperature-independent component of the enthalpy in addition to the traditionally recognized heat capacity of thermodynamics. Specifically, the temperatureindependent component of the chemical bond enthalpy is the inherent strength of the bond as measured at 0 K.
In biological interactions, the innate thermodynamic quantities [9,10,11,12,13,14,15,16], in particular the innate temperature-invariant enthalpy, represent differences in chemical bonding energy of products minus reactants, and thus control the heat of reaction at 0 K. The net difference between the temperature-dependent heat capacities of reaction [∆Cpº(T)], taken as products minus reactants, dictates how behaves at high temperatures.
o reaction ∆H We believe that we have uncovered a type of physical behavior characteristic of living systems that may prove to be remarkably general. It is, of course, known that living systems can live and operate optimally only at a sharply defined temperature, or over a limited temperature range at best. We assert that this implies that basic biochemical macromolecular interactions exhibit a well-defined negative free energy minimum (that is to say, favorable) as a function of temperature. Such a situation is not common in simple chemical systems, where a monotonic change of ∆Gº, ∆Hº, and K eq over an experimental temperature range is typical.
We have found that the critical factor driving the pair-wise hydrophobic interaction of leucine-isoleucine, and 35 similar dipeptides over the temperature range of 273-333 K reported by Nemethy and Scheraga in 1962[8] is a temperature-dependent heat capacity change of reaction, ∆Cpº(T) reaction , which is positive at low temperature but switches to a negative value at a temperature well below the ambient range. This change of sign of the critically important ∆Cpº(T) reaction (product minus reactants) has such significant consequences that we refer to it as a "thermodynamic molecular switch". It determines the behavior patterns of the Gibbs free energy change, and hence a change in the equilibrium constant, K eq , and/or spontaneity. Note that ∆Gº(T) reaction = -RT ln K eq , so that the roles of ∆Gº and K eq are rigidly coupled. The subsequent, mathematically predictable changes in ∆Hº, ∆Sº, ∆Wº, and ∆Gº which arise as a result of this thermodynamic molecular switch are demonstrated in this pair-wise, sequence-specific hydrophobic interaction [5,6,7].
Based on Chun's development of the Planck-Benzinger methodology, the change in inherent chemical bond energy at 0 K, ∆Hº(T 0 ), is 3.0 kcal mol -1 for Leu-Ile, about three times smaller than Nemethy and Scheraga's approach [8]. At <T m >, the thermal agitation energy, In all biological interactions, ∆H o (T) and ∆S o (T) are positive at low temperature. As reaction temperature increases, both ∆H o (T) and ∆S o (T) become negative, creating negative Gibbs free energy minimum: that is process 1 of the chart goes to process 2, creating process 3 (see Table 1). The change of sign in ∆Cp o (T) reaction leads to a true negative minimum in the Gibbs free energy of reaction, that is ∆Cp o (+) → ∆Cp o (-), designated as a thermodynamic molecular switch. It determines the behavior patterns of the Gibbs free energy change, and hence a change in the equilibrium constant, K eq , and/or spontaneity. is 12.8 kcal mol -1 , about five times greater than ∆Hº(T 0 ). This pair-wise, sequence-specific hydrophobic interaction is highly similar in its thermodynamic behavior to that of other biological systems, except that the negative Gibbs free energy change minimum at <T s > occurs at a considerably higher temperature, 355 K compared to about 300 K. The melting temperature, <T m >, is also high, 465 K compared to 343 K in a biological system. The implication is that the negative Gibbs free energy minimum at a well-defined <T s > has it origin in the hydrophobic interactions, which are highly dependent on details of molecular structure [6,7].
We have shown in our unpublished work the existence of a thermodynamic molecular switch in the interactions of 35 dipeptides over the temperature range of 273-333 K as initially reported by Nemethy and Scheraga [8] and which we have subsequently analyzed by the Planck-Benzinger methodology. In the hydrophobic interaction of each of 35 dipeptide pairs, a change of sign in ∆Cpº(T) reaction leads to true negative minimum in the Gibbs free energy of reaction, ∆Gº(T) reaction and hence a maximum in the related K eq , as shown in Fig. 1 and Table 1. Indeed, all interacting biological systems we have examined using the Planck-Benzinger methodology have shown such a thermodynamic molecular switch, suggesting that its existence may be universal [1,2,3,4,5,6,7].
A closer look at a single example, the pair-wise sequence-specific hydrophobic interaction of leucine-isoleucine, will demonstrate the significant differences when the Nemethy-Scheraga data are treated by the Planck-Benzinger methodology.

THE GIAUQUE FUNCTION AND PLANCK-BENZINGER THERMAL WORK FUNCTION
One form of the free energy function, the Giauque function [20,21,22,23,24], has been extensively used in chemistry and physics.  An equivalent formulation has recently found application in biochemical literature as the Planck-Benzinger thermal work function [1,2,3,4,5,6,7,9,10,11,12,13,14,15,16], where the application to a given situation is quite different. Here the Planck-Benzinger thermal work function, ∆W(T) = ∆H(T 0 )−∆G(T) [19], represents the strictly thermal components of any intra-or intermolecular bonding term, that is, energy other than the inherent difference of the 0 K portion of the interaction energy. The latter is the only energy term in constant pressure processes at absolute 0 K. Thus ∆W(T) expresses completely the thermal energy difference of the process involved. Application of the Planck-Benzinger thermal work function permits the separation of 0 K energy differences and energy differences associated with heat capacity integrals for a fuller understanding of reaction energies.

THE PLANCK-BENZINGER APPROACH
In order to analyze the thermodynamic processes operating in a pair-wise hydrophobic interaction such as leucine-isoleucine, it is necessary to extrapolate the thermodynamic parameters over a much broader temperature range. The enthalpy, entropy, and heat capacity terms are evaluated as partial derivatives of the Gibbs free energy function defined by Helmholtz-Kelvin's expression [25,26]. In continuing studies on dozens of interacting protein systems, it has been shown in our laboratory that the third-order polynomial function provides a good fit in the temperature range accessible in biochemical systems. In fact, it is shown to be correct in the low-temperature limit. The rationale for selecting the linear and nonlinear third-order(T 3 model) polynomial functions for ∆G o (T) = α + βT 2 + γT 3 (macromolecular interaction) and ∆H (T) = α + βT 3 + e γT (protein unfolding) are found in the fundamentals of relevant quantum theory [25]. The constraints for this fitted model of the Gibbs free energy change as a function of temperature are these: the enthalpy and Gibbs free energy must intersect at 0 K at a positive value and the slope must be zero where ∆G o (T 0 ) = ∆H o (T 0 ) at 0 K as T∆S o approaches zero. The cubic term is necessitated by quantum mechanical considerations, and the presence of the quadratic term is dictated by the data rather than theory [1,2,3,4,5,8,9,10,11,12,13,14,15].

DETERMINATION OF THE GIBBS FREE ENERGY CHANGE AS A FUNCTION OF TEMPERATURE
Nemethy and Scheraga's theoretical treatment of hydrophobic interaction [8,28] based on a statistical mechanical theory of the thermodynamic properties of liquid water [29] and of aqueous solutions of hydrocarbon [30] led to determination of the maximum strength between isolated hydrophobic side chains. ∆Gº(T) data for a pair-wise hydrophobic interaction of maximum strength were computed from Nemethy and Scheraga's Table II [8] based on T 2 model of ∆Gº = a + bT + cT 2 over a temperature range of 273-333 K, knowing the polynomial coefficients a, b, c [8,28]. Other thermodynamic parameters were are follows: Values for the Gibbs free energy change as a function of temperature were computed for the sequence-specific hydrophobic interactions of leucine-isoleucine and compared with values determined using our T 3 model of ∆Gº(T) = α + βT 2 + γT 3 [5].

COMPUTING HYDROPHOBIC INTERACTION BY THE PLANCK-BENZINGER METHOD
In order to extrapolate down to 0 K, it is necessary to consider the normal solution states of molecules. Here the 0 K limit would presumably refer to the glassy condition [1,2,3,4,5,6,7,13,14,15], that is a condition with all thermal agitation frozen out, but retaining the general physical chemical properties of solution -since a pure crystalline form of macromolecules is rarely encountered in practice.
The approach that we follow requires exact determination of K eq for the relevant biochemical processes as a function of absolute temperature. Of critical importance, however, is the use of precise, correctly formulated expressions for as a function of temperature as shown in Eq. 3. In this treatment, the Gibbs free energy data, as shown in Fig. 2A, were fitted to a threeterm linear polynomial function in the 273-333 K temperature range, the range in which computations have been conducted as shown in Table 2.
Once evaluated as shown in Fig. 2A (see Table 2), the coefficients α, β, and γ were fitted to other thermodynamic parameters. ∆Hº(T), ∆Cpº(T), T∆Sº(T), and ∆Wº(T) are defined by Helmholtz-Kelvin expression as follows:  Table II [8] in the temperature range 273-333 K, using the general linear model (T 3 model) procedure of statistical analysis of IMSL subroutine. The solid line represents fitted data. F = 0.0001, thus the goodness of fit of the experimental data was 99.9% or better in each case [31,32]. Note: Compiled using the general linear model procedure analysis of IMSL subroutine adapted for use in KaleidaGraph 3.5.1. Values for these ∆H o (T 0 ) vary by less than 0.01%. F = 0.0001, thus the goodness of fit of the computed data was 99.99% or better in each case. Each data point between 0 and 350 K was evaluated with extrapolation of F-statistics [31,32].
To extrapolate down to 0 K, it is necessary to consider normal solution states of macromolecules. Here the 0 K limit would presumably refer to the glassy condition [12], as described earlier.
Data for the Gibbs free energies were fitted to a model and the derived quantities for other thermodynamic parameters were calculated and plotted using the International Mathematical Subroutine Library (IMSL) software for linear and nonlinear polynomial regression analysis. Each equation was interactively executed in steps of 1 K, and the values plotted and overlaid for each set of experimental conditions. This IMSL subroutine was adapted for use in KaleidaGraph 3.5.1 and each data point was evaluated with extrapolation of F-statistics in an IBM personal computer [31,32], as shown in Figs. 2B-D.
A built-in restriction in the extrapolation procedure is that the values for ∆Hº(T) and ∆Gº(T) determined from the polynomial functions intersect at 0 K with zero slope on a thermodynamic plot, thus obeying Planck's definition of the Nernst heat theorem [17,25]. By our definition, the value of ∆Hº(T 0 ) will be positive. Other polynomial functions failed to meet all three restrictions of ∆Hº(T) and ∆Gº(T) intersecting at 0 K with zero slope, and ∆Hº(T 0 ) being positive and thus were discarded.
It is clear that a temperature-dependent model simpler than the third-order Gibbs polynomial model (T 3 model) cannot be used at low temperature, and has been found to be unacceptable at room temperature; therefore it would be reasonable to apply a more complex model in the intervening temperature region only if the facts demand a more complex fit. In fact, the facts do not require a different function; the model as described has been found to have both strong correlative power and very good predictive power.
The fitted thermal data [∆Gº(T), ∆Hº(T), ∆Wº(T), T∆Sº(T), and ∆Cpº(T)] are reasonable not only over the measured experimental range (near room temperature) but also in the low temperature limit. The fitting curves do nothing strange in the experimentally inaccessible region, but rather smoothly approach the low temperature limit. The thermal dependency presented here is the best (and essentially the only) approach known.  In the thermodynamic methods based on our development of Planck-Benzinger methodology [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15], computed values must be in agreement with experimentally obtained K eq values [and therefore with values of ∆Gº(T) of reaction] of biological systems over a temperature range of 273−343 K.  Fig. 2B, over the temperature range of 300-400 K, with the magnitude of the y-axis reduced to 4 to -4 kcal mol -1 K -1 .

WHY IS HELMHOLTZ-KELVIN'S EXPRESSION CONSIDERED TO BE A CONTINUOUS FUNCTION?
It is true that a change of phase for a specific chemical substance causes an infinite discontinuity in Cp, a finite discontinuity in S and H, and an abrupt change in the slope of G vs. T. Similar relevant data for phase changes in biological reactions are generally lacking; in our work the reference condition at 0 K is taken as a glassy solid. In this case, all thermodynamic functions are continuous down to 0 K.
Our evaluation of the innate temperature-invariant enthalpy for hydrogen-bonded water, in equilibrium with nonhydrogen bonded water molecules, is based on Helmholtz free energy data reported by Nemethy and Scheraga [29]. The entropy of the system appears to remain independent of temperature, suggesting that there is no significant temperature-dependent difference in the degree of orientation between unbound and hydrogen-bound water molecules in equilibrium in the system. A similar conclusion could be reached based on the dielectric relaxation of water as a function of temperature, as reported by Collie et al. [33]. This implies that there is a nonzero entropy difference for the transformation between water monomer and n-mer at 0 K (unpublished results [13]). As with most small molecules, the thermal agitation energy is minimal. ∆W(T) = T∆S(T), ∆H(T) = ∆H(T 0 ) over the entire temperature range from 0 K to the temperature of interest. As already noted, the assumed reference condition is glassy ice at 0 K, which probably does have a small ∆S(T 0 ). However, T∆Sº(T 0 ) = 0, in any event. In the case of ∆Gº(T) of formation, phase transition is indeed important. When dealing with ∆Gº(T) of reaction, no phase transition is taking place. In this case, all thermodynamic functions are continuous.  Table 3). The change in inherent chemical bond energy at 0 K, ∆H o (T 0 ), is 2.4 for Leu-Ile. In this pair, as with others we have examined, the value of ∆H o (T 0 ) decreases as the length of the hydrophobic side chain decreases. <T s > remains constant at 355 K and <T m > is about 470 K in each system. It is apparent that the strength and stability of the hydrophobic interaction is determined by the packing density of the side chains. At <T m >, the thermal agitation energy is about five times greater than ∆H o (T 0 ). Additionally, the thermal agitation energy for the same series of pair-wise sequence-specific hydrophobic interaction, evaluated at <T m >, decreases in the same order, that is as the length of the side chain decreases.   Table 3). As shown in Table 1, in this hydrophobic interaction, at low temperature, ∆H o (T) and ∆S o (T) are both positive (entropy-driven process), becoming negative as temperature increases (enthalpydriven process), whereas ∆G o (T) changes from positive to negative then reaches a negative value of maximum magnitude at <T s >, and finally becomes positive as temperature increases (see Figs. 2B-D). That is, process 1 goes to process 2, creating cooperative enthalpy-entropy compensation between <T h > and <T m >, as shown in Table 1 and Fig. 1.

NEMETHY AND SCHERAGA'S METHOD OF ANALYSIS
Fitting the T 2 model of ∆Gº = a + bT + cT 2 , for the Gibbs free energy as a function of temperature of Nemethy and Scheraga [8,28] and our T 3 model of ∆Gº(T) = α + βT 2 + γT 3 , based the Planck-Benzinger methodology, it is apparent from Fig. 3A that both models fit the theoretically computed ∆G o (T) data equally well over the temperature range of 273-333 K. As shown in Fig.   3B, the values for ∆H o (T) and ∆G o (T) determined from a linear T 2 polynomial function do not intersect at 0 K with zero slope and thus are inconsistent with Planck's definition of the Nernst heat theorem [17]. The value of ∆H o (T 0 ), shown in Table 3, was found to be 7.53 kcal mol -1 , five times greater then ∆H o (T 0 ) obtained by the Planck-Benzinger approach. The thermal agitation energy was -8.45 kcal mol -1 , slightly lower than that obtained by our approach (see Fig. 3B).  Table II [8] in the temperature range 273-333 K, using the general linear polynomial model (T 2 model: ∆G o (T) = a + bT + cT 2 ) procedure of statistical analysis of IMSL subroutine. The solid line represents fitted data. F = 0.0001, thus the goodness of fit of the experimental data was 99.99% or better in each case.
In Nemethy and Scheraga's analysis of thermodynamic parameters, fit of the data for heat capacity vs. temperature when extrapolated to 0 K is not fully adequate, as shown in Figs. 3C-D, therefore an unacceptable prediction of entropy is to be expected. Furthermore, no thermodynamic molecular switch is observed, rather ∆Cp  Table 4). FIGURE 3B. Thermodynamic plot of Nemethy-Scheraga T 2 model for Leu-Ile hydrophobic interaction. Each data point between 0 and 500 K was evaluated with extrapolation of F-statistics [31,32]. The solid line represents fitted data, F = 0.0001, thus the goodness of fit of the experimental data was 99.99% or better in each case. The values for ∆H o (T) and ∆G o (T) determined from a linear T 2 polynomial function do not intercept at 0 K with zero slope and thus are inconsistent with Planck's definition of the Nernst heat theorem [16]. Fig. 3B, over temperature range of 0-500 K, with the magnitude of the y-axis reduced to 0.03 to -0.04 kcal mol -1 . No thermodynamic molecular switch is observed.

FIGURE 3D.
A close-up view of a portion of Leu-Ile hydrophobic interaction as shown in Fig. 3B, over the temperature range of 300-400 K, with the magnitude of the y-axis reduced to 4 to -4 kcal mol -1 K -1 . In each case the respective value at a given temperature is determined in a straightforward way using the following expression:

TABLE 4 Thermodynamic Molecular Switch in ∆Cp o (T) and ∆S
Note that the value of the enthalpy change at 0 K, ∆H o (T 0 ), is important, but distinct and separate from the thermal agitation term. Many authors have entirely ignored this, particularly when dealing with biological systems.  (Planck-Benzinger) Leu-Ile With proper rearrangement, this equation also yields the Giauque and Planck-Benzinger thermal work functions [1,2,3,4,5,6,7,9,10,11,12,13,14,15,16] Table 1). In this thermodynamic switch unique to and characteristic of hydrophobic interaction, ∆Cp o (T) changes from positive to negative at <T Cp >, shown in Fig. 2C. As seen in Table 4, the Gibbs free energy change switches sign at <T h > and <T m >, where ∆G o (T) = 0, and reaches a negative minimum at <T s >. The <T s > value is 355 K at T∆S o (T) = 0 (see Table 1

THERMODYNAMIC MOLECULAR SWITCH IN SEQUENCE-SPECIFIC HYDROPHOBIC INTERACTIONS
For the pair-wise, sequence-specific hydrophobic interaction of Leu-Ile we have examined,  Table 5). The significance of the negative Gibbs free energy minimum in equilibria of sequencespecific dipeptides and on a wider scope, of biological systems, is that one is dealing with a true condition of stability, which is the maximum K eq of the reaction. It is no surprise that the temperature of maximum stability should be found over a broad temperature range in biological systems, since such systems are known to exist optimally from arctic temperatures to several hundred degrees Centigrade in fumerals of the ocean floor.
These data demonstrate that the critical factor is a temperature-dependent heat of reaction, ∆Cp o (T), which is positive at low temperature but switches to a negative value at <T Cp > = 180 K, the ambient temperature range. Note that this thermodynamic switch in the sign of ∆Cp o (T) reaction determines the behavior patterns of the Gibbs free energy change and hence a change in the equilibrium constant, K eq , and/or spontaneity, as shown in Figs It is clear that nonlinear temperature dependence of ∆Cp o (T) reaction is the essence of the thermodynamic molecular switch. In the present instance, the hydrophobic interactions seem to be the best candidates for such behavior. At present the ultimate range of applicability of the thermodynamic molecular switch concept is not known.
In the pair-wise sequence-specific hydrophobic interaction of Leu-Ile, the value of ∆H o (T 0 ) is 2.40 kcal mol -1 . It is apparent from results using the Planck-Benzinger approach that the strength and stability of the hydrophobic interaction is determined by the packing density of the side chains. At <T m >, the thermal agitation is about five times greater than ∆H o (T 0 ) in these systems. Additionally, the thermal agitation energy for the same series, evaluated at <T m >, decreases in the same order, which is as the length of the side chain decreases.
This pair-wise, sequence-specific dipeptide hydrophobic interaction is highly similar in its thermodynamic behavior to that of other biological systems, except that the negative Gibbs free energy change minimum at <T s > occurs at a considerably higher temperature, 355 K compared to about 300 K. The melting temperature, <T m >, is also high, 470 K compared to 343 K in a biological system. The implication is that the negative Gibbs free energy minimum at a welldefined <T s > has its origin in the sequence-specific hydrophobic interactions, which are highly dependent on details of molecular structure. Since detailed statistical mechanical computations on our system have not been attempted, no further clarification can be provided at this time.
In addition to the specific dipeptide interaction previously described here, we have shown in our unpublished work the existence of a thermodynamic molecular switch in the interactions of 35 dipeptides wherein a change of sign in ∆Cp o (T) reaction leads to true negative minimum in the Gibbs free energy of reaction and hence a maximum in the related K eq . Indeed, all interacting biological systems we have examined using the Planck-Benzinger methodology have shown such a thermodynamic molecular switch, suggesting that its existence may be universal [3,4,5,6,7,8,9,10]. All interacting biological systems that we have thus far examined using the Planck-Benzinger approach point to the universality of this thermodynamic switch.