# Origin of the Native Driving Force for Protein Folding

\[xx Preprint NCU/CCS-1998-0920; NSC-CTS-981002

## Abstract

We derive an expression with four adjustable parameters that reproduces well the $`20\times 20`$ Miyazawa-Jernigan potential matrix extracted from known protein structures. The numerical values of the parameters can be approximately computed from the surface tension of water, water-screened dipole interactions between residues and water and among residues, and average exposures of residues in folded proteins.

\]

Protein structure and design is a very important topic in life science where physics and mathematics are indispensable to its understanding . Recently Li et al.  pointed out some highly interesting and unexpected properties of Miyazawa and Jernigan’s $`20\times 20`$ potential matrix ($`M`$) for protein structure . This matrix, whose elements are statistically deduced pair-wise interaction potential energies among the twenty types of amino acids in proteins of known structure, has been widely applied to protein design and folding simulations . Li et al. noticed that $`M`$ has a highly accurate leading principal-component representation: variations of the elements of $`M`$ from their mean can be expressed in terms of only the two leading eigenvalues of $`M`$ and the eigenvector $`\stackrel{}{q}`$ of the leading eigenvalue such that

$$M_{ij}c_2q_iq_j+c_1(q_i+q_j)+c_0,$$
(1)   
where $`i`$ and $`j`$ label the 20 amino acids, and $`c_0=1.38`$, $`c_1=5.03`$ and $`c_2=7.40`$, in units of $`RT`$, the gas constant times (room) temperature.

Two features of the right-hand-side of Eq. (1) stand out: 1) Not all residue-dependent terms are genuine two-body interactions; the $`c_1`$ terms represent one-body, mean-field potential energies. 2) Both the two-body $`c_2`$ terms and the one-body $`c_1`$ terms depend on the same set of $`q`$’s. Numerically, because the magnitudes of the $`q`$’s are small, the $`c_1`$ terms dominate over the $`c_2`$ term. This is consistent with the widely held notion that the earliest and fastest part of a protein folding process is by and large controlled by the hydrophobicity  of the residues. Tables I and II show that indeed $`q`$ is moderately correlated with the hydrophobicities ($`\mathrm{\Delta }G`$) . The product, pairwise form of the two-body terms reminds one of dipole-dipole interaction, and this in turn would imply a connection between the one-body terms and the dipole moments of the residues. Tables I and II also show a noticeable correlation between $`q`$ and the dipole moments ($`Q`$) of the side-chains of the residues . In the rest of the paper we will derive an expression for the MJ matrix in terms of an average “bare” residual solvation energy (for a hypothetical residue with vanishing dipole), interactions between the dipole moments of the residues and water molecules, and the degree of exposure to water (expressed as its complement, the burial factor) of a residue in a folded protein. We show that except for the burial factor of the residues the other three adjustable parameters appearing in the expression all have clear physical meanings with numerical values that can be computed approximately. The average burial factors for hydrophobic and hydrophilic residues that emerge from our analysis of the MJ matrix are 0.8 and 0.2, respectively (they are related and should approximately sum to 1). In this paper energy will be given in units of $`RT=0.60\mathrm{kcal}/\mathrm{mol}=4.2\times 10^{21}J`$ and dipole moments will be given in Debyes ($`D`$).

Dipole-dipole interaction. The interaction in vacuum between two electric dipoles $`\stackrel{}{Q}_i`$ and $`\stackrel{}{Q}_j`$ separated by $`\stackrel{}{R}_{ij}=\widehat{n}R_{ij}`$ is $`V_{ij}=(\stackrel{}{Q}_i\stackrel{}{Q}_j3(\widehat{n}\stackrel{}{Q}_i)(\widehat{n}\stackrel{}{Q}_j))/(4\pi ϵ_0R_{ij}^3)`$. If the carriers of the dipoles are relatively unconstrained we expect attraction and $`|\mu _r||Q_i||Q_j|V_{ij}0`$, where $`|\mu _r|=D^2/2\pi ϵ_0R_{ij}^3`$. In what follows, $`Q_i`$, $`i=1,\mathrm{},20`$ is the dipole moment of the $`i^{th}`$ side-chain, and $`Q_w`$ is the dipole moment of a water molecule. For residue-residue interaction, taking the inter-side-chain distance to be $`R_{ij}R_06.5\AA `$ , and recalling that an electron-positron pair separated by one $`\AA `$ is equal to 4.8 $`D`$, we have $`|\mu _r|0.172`$ $`(RT)`$, which may be viewed as a maximum value for the coupling since in a real setting it is expected to be weakened owing to the presence of water molecules.

One-body terms. Let $`E_0`$ be the average bare surface-dependent solvation energy of a residue in water when the residue-water dipole interaction is not taken into account; $`N_w`$ be the average number of water molecules in contact with a residue; $`\mu _w`$ be the average effective dipole-dipole coupling between the $`i^{th}`$ residue and a water molecule. Then, with residue-water interaction energy included and possible dependence of $`E_0`$, $`\mu _w`$ and $`N_w`$ on $`i`$ ignored, the residue-water interaction energy is $`E_i=\mu _wQ_iQ_wN_w+E_0\mu _wQ_i^{}Q_wN_w`$, where for convenience we write $`Q_i^{}Q_i+Q_0`$ and $`Q_0=E_0/(\mu _wQ_wN_w)`$. A hydrophobic (hydrophilic) residue would have $`E_i>0`$ ($`E_i<0`$). If $`N_i`$ is the number of the type $`i^{th}`$ residues in a peptide, then the energy of an unfolded peptide in water is $`U=_iN_iE_i`$. Suppose that after folding $`\mathrm{\Delta }N_i`$ fewer $`i^{th}`$ residues are exposed to water. Then the binding energy of the folded relative to the unfolded state is $`\mathrm{\Delta }U=_i\mathrm{\Delta }N_iE_i`$. The negative sign means that in folding, the peptide will maximize (minimize) those $`\mathrm{\Delta }N_i`$ whose $`E_i`$ are the most positive (negative), subject to the constraint of polymeric nature of the peptide.

Relation between $`q`$ and $`Q`$. Equating $`\mathrm{\Delta }U`$ with the binding energy obtained from Eq. (1) by summing the one-body terms over all pairs we have

$$\mathrm{\Delta }Uc_1N_c\underset{i}{}N_iq_i^{}=\mu _wQ_wN_w\underset{i}{}Q_i^{}\mathrm{\Delta }N_i$$
(2)   
where $`q_i^{}q_iq_0`$, $`q_0`$ is a constant and $`N_c`$ is the average number of contacts a residue has in a folded state. Matching the $`i`$-dependent terms we have

$$c_1q_i^{}\xi _iQ_i^{},\xi _i=\mu _w(\mathrm{\Delta }N_i/N_i)(N_wQ_w/N_c).$$
(3)   
Because in a folded protein proportionally more hydrophobic ($`h`$) residues than polar ($`p`$) residues will be hidden from water, one expects $`\mathrm{\Delta }N_i/N_i`$, hence $`\xi _i`$, to have a strong residual dependence. To minimize the number of parameters we allow $`\xi _i`$ to have only two values: $`\xi _h`$ and $`\xi _p`$, and have them determined by separate linear fits to $`q`$’s belonging to hydrophobic and hydrophilic residues, respectively. Excluded in the fits are residues whose hydrophobicities are ambivalent  \- Tyr, Ala, Gly, Thr, Ser and Pro. Demanding that the two fits have the same intercepts we obtain

$$q_0=0.055,Q_0=2.9;\xi _h=0.56,\xi _p=0.14$$
(4)   
The linear correlation between $`q`$ and $`\xi Q^{}`$ over the complete set of 20 residues - following  and , the first eight amino acids in Table I are taken to be hydrophobic - is 0.949, which is dramatically better than the correlation between $`q`$ and $`Q`$, see Fig. 1(a) and Table II.

The burial factor. Since on average the numbers of hydrophobic and polar residues in a protein are approximately equal and about half of all residues are buried in the core, we have $`N_hN_p`$, $`\mathrm{\Delta }(N_h+N_p)/(N_h+N_p)1/2`$ and hence $`\mathrm{\Delta }N_p/N_p1\mathrm{\Delta }N_h/N_h`$. From the ratios of the two $`\xi `$’s we thus deduce the burial factors for hydrophobic and polar residues, respectively, to be

$`\mathrm{\Delta }N_h/N_h0.80,\mathrm{\Delta }N_p/N_p0.20.`$ (5)   
That is, our analysis of the MJ matrix suggests that on average four times as many hydrophobic residues are buried in the core than are polar residues.

Two-body terms. We define the true two-body part of the MJ matrix to be the matrix minus the one-body and constant part of Eq. (1): $`M_{ij}c_0c_1(q_i+q_j)`$. This two-body part is again well approximated by $`c_2^{}q_iq_j`$, $`c_2^{}=10.7`$, with which it has a linear correlation of 0.832. When $`c_2^{}q_iq_j`$ is re-expressed in terms of $`Q^{}`$ using Eq. (3) the shift $`q_0`$ induces an additional one-body term such that

$$M_{ij}C_2\xi _i\xi _jQ_i^{}Q_j^{}+c_1^{}(q_i+q_j)+\mathrm{const}.$$
(6)   
Where $`C_2=c_2^{}/c_1^2=0.423`$ and $`c_1^{}=c_1c_2^{}q_0=5.62`$. The linear correlation between $`M_{ij}c_1^{}(q_i+q_j)`$ and $`\xi _i\xi _jQ_i^{}Q_j^{}`$ is 0.681, see Fig. 1(b). Given that the dipole moments and $`\xi _h`$ and $`\xi _p`$ are predetermined, the first term on the right-hand-side of Eq. (6) is a one free parameter ($`C_2`$) fit to 210 pieces of “noise” in the MJ matrix. The mediocre quality of the correlation nevertheless suggests that the two-body term cannot be explained by dipole interactions alone; interactions depending on charge and polarizability may need to be included. The inclusion of such terms may cause the two-body term to deviate from having the simple $`qq`$ form suggested by in Eq. 1. Owing to its relative small magnitude such a deviation should be tolerable to the original MJ matrix.

MJ matrix in terms of $`Q^{}`$. Re-expressing the one-body term in Eq. (6) in terms of $`Q^{}`$ and rationalizing notations by writing $`\mu _{ij}=C_2\xi _i\xi _j`$ and $`\xi _i^{}=\xi _ic_1^{}/c_1`$ we finally have

$$M_{ij}\mu _{ij}Q_i^{}Q_j^{}+\xi _j^{}Q_j^{}+\xi _i^{}Q_i^{}+\text{const.}$$
(7)   
where $`\mu _{hh}=0.13`$, $`\mu _{hp}=0.032`$, $`\mu _{pp}=0.0078`$, $`\xi _h^{}=0.63`$ and $`\xi _p^{}=0.15`$. The two sides of the equation has a linear correlation of 0.922, see Fig. 1(c). Since $`Q_i`$ is either zero or positive, the negative values of $`\mu _{ij}`$ imply that the dipoles mostly succeed in causing the residues to lower their energies. That is, even in a folded state the the residues appear to be sufficiently unrestricted to find optimum orientations. To the extent that the dipole moments of the side-chains are not free parameters, the expression on the right-hand-side is a four parameter fit - $`C_2`$, $`E_0`$, $`\mathrm{\Delta }N_h/N_h`$ and $`\mu _w`$ (see below) - to the complete MJ matrix.

Residue-residue dipole coupling. By definition $`\mu _{ij}(\mathrm{\Delta }N_i/N_i)(\mathrm{\Delta }N_j/N_j)`$. With $`\mathrm{\Delta }N_i/N_i`$ describing the percentage of buried residues in a folded protein, the inequalities $`|\mu _{pp}|<|\mu _{hp}|<|\mu _{hh}|<|\mu _r|`$ correctly take into account the dielectric property of water: the coupling between residues shielded from water is stronger than that between residues that are not. The magnitude of the weighed average of the residue-residue coupling, $`\overline{\mu }_{ij}=(7\mu _{pp}+6\mu _{hp}+7\mu _{hh})/20=0.041`$, is about four times less than the bare coupling strength of $`|\mu _r|=0.172`$.

Water-residue coupling. We can obtain the effective water-residue coupling from the relation $`\xi _i=\mu _w(\mathrm{\Delta }N_i/N_i)(N_wQ_w/N_c)`$ given earlier. Using the value 6.5 $`\AA `$ for the average effective diameter of a residue and the value 2 $`\AA `$ for the diameter of a water molecule, we estimate that a residue may have a maximum of 12 residue contacts and 57 water molecule contacts. In practice the number of contacts is encumbered by the presence of the peptide backbone and geometric constraints, such that in fact $`N_c7\text{[3]}`$. We therefore scale $`N_w`$ down to $`35`$. With $`Q_w=1.85D`$, we deduce from Eqs. (5) and (7) that $`\mu _w0.076(RT)`$. The negative sign of $`\mu _w`$ is consistent with the notion that the presence of dipole in a residue reduces its hydrophobicity. Taking the average water-residue distance to be 4.25 $`\AA `$ we expect the bare water-residue coupling to be $`(6.5/4.25)^3=3.5`$ times stronger than the bare residue-residue coupling. However, in an unfolded state the residues are completely exposed to water. We therefore expect the approximate relations $`|\mu _{pp}|<|\mu _w|/3.5|\mu _{hp}||\overline{\mu }_{ij}|<|\mu _r|`$, which are satisfied.

Solvation energy, surface tension and hydrophobicity. With $`\mu _w`$ and $`Q_0`$ extracted from the data we now find the bare solvation energy to be $`E_0=\mu _wQ_0Q_wN_w=14.6RT`$. Although hydration is an exceedingly complex process and is not fully understood, the effective surface tension of water, or surface free energy cost to water forced to sit against a hydrophobic surface has been estimated to be $`\sigma =40`$ erg/cm<sup>2</sup> . For a residue of diameter $`R_0`$ the free energy cost is $`W=4\pi (R_0/2)^2\sigma =13RT`$, which is reasonably close to the value of $`E_0`$. The fact that a good fit to the MJ matrix demands that $`E_0`$ enters $`\mathrm{\Delta }U`$ in Eq. 2 multiplied by $`\mathrm{\Delta }N_i`$ is indication that $`E_0`$ needs to be surface energy. When the water-residue dipole interaction energy is included, the total solvation energies $`E_i`$ of the residues then delineate into groups with distinct hydrophobicities, with the seven most hydrophobic (hydrophilic) having an average solvation energy of $`13.2RT`$ ($`9.3RT`$).

Very recently Keskin et al.  re-analyzed the MJ matrix and derived the approximation (for ease of discussion the $`W_i^{}`$ used here has an additional negative sign relative to that in ): $`M_{ij}\mathrm{\Delta }W_{ij}^{}+W_i^{}+W_j^{}+\mathrm{const}.`$, where the one-body term $`W^{}`$ is essentially defined as the mean-field of $`M_{ij}`$ and $`\mathrm{\Delta }W_{ij}^{}`$ is a four parameter fit to $`M_{ij}`$ minus its mean-field. The analysis confirms the dominance of the one-body term in the MJ matrix. The overall fit to the MJ matrix, with a correlation of 0.99, is excellent and the fit to the two-body part is about the same as that given by the dipole picture: the correlation between $`M_{ij}W_i^{}W_j^{}`$ and $`\mathrm{\Delta }W_{ij}^{}`$ is 0.67. Not surprisingly $`W^{}`$ and $`q`$ are closely related. The expression $`\eta c_1q+1.16`$, with scale factor $`\eta =1.17`$, reproduces $`W^{}`$ with a linear correlation of 0.997. The value of $`\eta `$ is mostly explained by the fact that the mean-field calculated from the right-hand-side of Eq. (1) is $`1.22c_1(q_i+q_j)`$. Incidentally, $`\eta c_1=5.89`$ is very close to the value of the renormalized coefficient $`c_1^{}=5.62`$ given in Eq. (6).

In Table I are listed values for $`q`$, $`Q`$, $`W^{}`$, $`\xi Q^{}`$, and hydropathy scales $`\mathrm{\Delta }G`$ (in units of $`RT`$) corrected for self-solvation for the side-chains of the twenty amino acids . Recall that $`\xi `$ contains the burial factor (see Eq. 4) and $`Q^{}`$ is $`Q`$ shifted by an amount proportional to $`E_0`$ (see Eq. 3). The pairwise linear correlation of the entries in Table I are given in column 2 of Table II. The correlation between $`\xi Q^{}`$ and $`W^{}`$ (and $`q`$) is very significantly better than that between $`Q`$ and $`W^{}`$ (and $`q`$). The linear relations connecting the solvation energy with $`\xi Q^{}`$, $`W^{}`$ and $`q`$: $`E_i(\mathrm{\Delta }N_i/N_i)/N_c=\xi _iQ_i^{}c_1(q_iq_0)(W_i^{}W_0^{})/\eta `$, where $`W_0^{}=0.71`$ is a shift, highlight the importance of taking into account the burial factor of a residue in a folded protein when interpreting the one-body terms of the MJ matrix.

The hydropathy scales shown in Table I are derived for side-chains in model peptides rather than in proteins. They include the effect of self-solvation that reduces the hydropathies of the polar side-chains , but does not include the effect of burial factor. This probably explains why, as seen in Table II, the $`\mathrm{\Delta }Gq`$, $`\mathrm{\Delta }GW^{}`$, $`\mathrm{\Delta }GQ`$ and $`\mathrm{\Delta }G\xi Q^{}`$ correlations are of similar quality.

The $`q`$ and $`W^{}`$ values of proline suggest it to be polar, while its $`Q`$, $`\xi Q^{}`$ and $`\mathrm{\Delta }G`$ values say it is ambivalent or even hydrophobic. The third column in Table II shows that the correlations listed either remain unchanged or improve when proline is excluded from the linear fit. The ambiguous hydrophobicity of this residue may be related to the fact that is has a looping structure.

We summarize our interpretation of Eq. (1) being a good approximation of the MJ matrix as follows. The one-body part, or hydrophobicity (or hydropathy) energy, is made up of two parts: free energy cost to water to accommodate the residue surface, and attractive dipole interaction between residue and water. Because polar residues have large dipole moments, hydrophobic residues have small or no moments and ambivalent residues have something in between, the hydropathic/hydrophobic energy is strongly attractive, weakly attractive and strongly repulsive for polar, ambivalent and hydrophobic residues, respectively. Residue-residue dipole interactions accounts for a sizable portion, but not all, of the two-body part. Aside from using the given dipole moments for the residues and having two burial factors, one each for the hydrophobic and polar residues, no residue-dependent adjustments were made in deriving Eq. (7), our rendition of Eq. (1). That is, we have not attempted a detailed fit of the MJ matrix. The correlation between the dipoles of the residues and $`q`$ becomes unequivocal and the strengths of the dipole couplings extracted from the MJ matrix become reasonable only when the burial factors are included in the formulation. That the factor is important reveals the dynamical nature of protein folding: strengths of interactions change as the folding progresses. Protein folding is a very complicated process that depends on many details and the MJ matrix does not tell its whole story. It does however contain the most basic structural information at the molecular level of those proteins whose structures are known. The success of the present analysis in understanding the main features of the MJ matrix gives us confidence that the model used here may provide a starting point for building a true potential suitable for use in a molecular dynamical description of early folding of protein in water.

The authors thank C. Tang, G.M. Crippen, Y. Duan, M. Wortis, D.C.Y. Lu and P.G. Luan for discussions and the referee for pointing them to reference  and for useful suggestions. HCL thanks the Physics Department of Simon Fraser University for hospitality in the Summer of 1998 during which part of the paper was written. This work is partly supported by grant NSC88-M-2112-008-009 from National Science Council (ROC).