text:
# Warm and cold Denaturation in the Phase Diagram of a Protein Lattice Model.

Under physiological condition, a protein adopts a unique three dimensional structure which is totally encoded in its amino acid sequence . Because, its biological function is closely related to this native structure, an important challenge in molecular modeling is to predict the native structure of a protein, given its amino acid sequence. It is now commonly assumed that the solvent plays a very important role in the folding process of proteins towards their native state . Its contribution to the conformational free energy should, therefore, be introduced in the computation of the energy in simulations of biological molecules .

Many proteins are in their native form between $`10^{}\mathrm{C}`$ and $`40^{}\mathrm{C}`$ . At higher temperatures, a reversible, warm denaturation occurs. For these proteins, a cold denaturation also occurs at low temperatures. For example, Privalov  showed experimentally that, under medium pH condition, the evolution of the heat capacity of myoglobin as a function of the temperature exhibits two peaks that are the signatures of two modifications in the order of the chain: it only folds at medium temperatures, and both a cold and a warm denaturation occurs.

A widely used class of models for studying theoretically protein folding is the lattice representation . Up until now, solvent effects were taken into account in these models using a parameter independent of the temperature, introduced in the intrachain calculations and favoring compact conformations. Dinner et al.  used the random energy model (REM)  to compute phase diagrams of protein-like chains. With this solvent model, they found three different regions: a native state, a globular state and a random coil state. Their phase diagrams, however, do not exhibit any sign characteristic of a cold denaturation. Several more sophisticated solvent models have been already proposed in the literature , but they have not yet been able to simulate a cold denaturation in a protein model.

Two theoretical analyses of disordered phases at low temperature have been proposed recently. Hansen et al.  used a zipper model equivalent to the 1-D Ising model, and obtained warm and cold transitions. De Los Rios and Caldarelli  proposed a Potts model that allows the observation of two peaks in the heat capacity profiles of short chains but the warm unfolding transition disappears as the length of the chain increases. Moreover, both models are homopolymer models, and, as such, do not exhibit regions in which the chain folds into a unique, compact, native, structure.

In this letter, we present a simple solvation model for REM, based on a study of the thermodynamic and the structural properties of the aqueous environment of a protein. We find a phase diagram which exhibits a native and a denatured regions with both warm and cold denaturation transitions. Our results are qualitatively in good agreements with the experimental findings of Privalov . We apply our solvation model to a 16-monomers chain constrained on a two-dimensional lattice for which 802075 different self avoiding walk conformations non equivalent by symmetry can be easily enumerated. We introduce a parameter $`B_s`$ in the model (more details are given below) which is related to the quality and the pH of the solvent and acts as a control parameter. At small $`B_s`$, the protein has a strong tendency to fold, while at large $`B_s`$ solvation becomes more effective and the protein tends to unfold.

The main results of our investigation are summarized in fig.1. Fig.1(a), (b), (c) and (d) shows the behavior of the specific heat $`C`$ and of the order parameter $`Q`$. Both quantities are plotted versus temperature (fig.1(a) and (b)) and versus temperature and solvation parameter (fig.1(c) and (d)). The two peaks on the heat capacity curves shown on fig.1(a) for $`B_s=12.0`$ and $`B_s=13.0`$ are the signatures of the modification of the order of the chain. A warm denaturation occurs at the temperature $`T_w`$ corresponding to the higher temperature heat capacity peak. The temperature of the other peak, when occurs, is the temperature of cold denaturation $`T_c`$. The sharp decay of the order parameter indicates that the transitions are first order type. The fact that, $`Q=1`$ between $`T_c`$ and $`T_w`$, shows that the chain is in a unique native compact structure. $`Q0`$ at temperature below than $`T_c`$ or above than $`T_w`$ then a huge number of extended conformations of the chain are relevant. One must note the particularity of the cold denaturation transition : the chain is in a disordered state at low temperature and becomes ordered in a unique structure as the temperature increases.

No peaks occurs, for values of $`B_s>11.2`$ (fig.1(c)) and the chain is always denatured (fig.1(d)). The peak corresponding to the cold denaturation do not appear for $`B_s<14.4`$ and, only the warm denaturation still occurs. The curves $`T_c(B_s)`$ and $`T_w(B_s)`$ are reported on the ($`B_s`$, $`T`$) diagram of the system (fig.1(e)). One sees that for $`(B_s,T)=(11.2,0.58)`$ the two curves intersect. At this particular point a naturation at constant temperature occurs by modifying the solvent condition without increasing the heat capacity. This transition seems occur without latent heat. Phase diagram for myoglobin obtained from result of Privalov are shown on fig.1(f). The good agreement between the calorimetric experiments and our theoretical results strongly supports our model of solvation that we introduce now.

The mechanism driving proteins towards compact structures is strongly coupled to the hydrophobic effect . The latter arises from a loss of solvent entropy when a hydrophobic monomer is transfered from the interior of the protein to the vicinity of the aqueous solvent (fig.2). In addition, intrachain interactions play the ultimate role for finding the unique, native structure among all the compact structures involving internal organization. Let us note $`E_{\mathrm{intr}}^{(m)}`$ the intrachain energy of the chain in the conformation $`m`$ and $`E_{\mathrm{solv}}^{(mm^{})}`$ the solvation energy resulting of the interactions between the solvent molecules in the conformation $`m^{}`$ and between solvent molecules and the chain units. The partition function of the system in thermal equilibrium reads :

$$Z(T)=\underset{m\mathrm{\Omega }}{}\underset{m^{}\mathrm{\Omega }^{}(m)}{}\mathrm{exp}\left(\frac{E_{\mathrm{intr}}^{(m)}+E_{\mathrm{solv}}^{(mm^{})}}{T}\right)$$
(1)
where $`\mathrm{\Omega }`$ is the conformational space of the chain and one must note that the conformational space of the solvent $`\mathrm{\Omega }^{}(m)`$ depends on the chain conformation $`m`$. The free energy of solvation is written :

$$F_{\mathrm{solv}}^{(m)}(T)=T\mathrm{ln}\underset{m^{}\mathrm{\Omega }^{}(m)}{}\mathrm{exp}\left(\frac{E_{\mathrm{solv}}^{(mm^{})}}{T}\right)$$
(2)
Then, it comes :

$$Z(T)=\underset{m\mathrm{\Omega }}{}\mathrm{exp}\left(\frac{F_{\mathrm{tot}}^{(m)}(T)}{T}\right)$$
(3)
where the total free energy of the conformation $`m`$ is simply : $`F_{\mathrm{tot}}^{(m)}(T)=E_{\mathrm{intr}}^{(m)}+F_{\mathrm{solv}}^{(m)}(T)`$

The chain is represented by a string of $`N`$ beads constrainted on a square lattice. The intrachain energy is given by the classical expression :

$$E_{\mathrm{intr}}^{(m)}=\underset{i=1}{\overset{N}{}}\underset{j>i}{\overset{N}{}}B_{ij}\mathrm{\Delta }_{ij}^{(m)}=\frac{1}{2}\underset{i=1}{\overset{N}{}}\underset{ji}{\overset{N}{}}B_{ij}\mathrm{\Delta }_{ij}^{(m)}$$
(4)
where $`\mathrm{\Delta }_{ij}^{(m)}`$ equals 1 if the monomer $`i`$ and $`j`$ are first neighbors on the lattice and the symmetric couplings $`B_{ij}=B_{ji}`$ are chosen at random in a gaussian distribution with a standard deviation equals 1. In contrast with previous works  the mean of the gaussian law is taken equals to 0 in order to generate attractive and repulsive interactions between the monomers.

Each empty site of the lattice is filled up with a group of solvent molecules which interacts with its four first neighbor sites on which are located either monomers of the chain or other groups of solvent molecules. Then, the free energy of solvation is written as the sum of interactions involving at least one solvent sites (solvent-monomer interactions and solvent-solvent interactions) :

$$F_{\mathrm{solv}}^{(m)}(T)=\underset{i=1}{\overset{N}{}}n_i^{(m)}f_i(T)+n_s^{(m)}f_s(T)$$
(5)
where $`n_i^{(m)}`$ is the number of solvent sites first neighbors of the $`i^{th}`$ monomer on the lattice and $`n_s^{(m)}`$ is the total number of first neighbors interactions between any pairs of solvent sites. Both quantities depend on the conformation of the chain. $`f_i(T)`$ is the contribution to the free energy of the interaction between a group of solvent molecules and the monomer $`i`$ and $`f_s(T)`$ is the free energy of interaction between two solvent sites. They are given by :

$$f_i(T)=T\mathrm{ln}\underset{j=1}{\overset{N_s}{}}\mathrm{exp}(B_i^{(j)}/T)\mathrm{and}f_s(T)=T\mathrm{ln}\underset{j^{}=1}{\overset{N_s}{}}\mathrm{exp}(B_s^{(j^{})}/T)$$
(6)
where $`B_i^{(j)}`$ is the energy of the bond between the $`i^{th}`$ monomer and a group of solvent molecules in the conformational state $`j`$ and the state $`j^{}`$ of the bond between two solvent sites has for energy $`B_s^{(j^{})}`$. The two summations run over the conformational states of the system composed of the solvent molecules of one group. We assume that, the number of states $`N_s`$ which can adopt the system of solvent molecules is the same in both cases, but, the solvation energy density functions of $`B_i`$ and $`B_s`$ are different (fig.2). We give below the arguments for choosing these functions.

We assume that the density energy function of the interaction between the solvent and a monomer is well represented by a gaussian law. The non-degenerated energy minimum of the $`N_s`$ states corresponds to the clatracle conformation (fig.2b) and at sufficiently low temperature, only this state is relevant. The mean of the gaussian law is taken as the solvation energy reference and is chosen equals to 0. The $`N_s`$ states, i.e. the $`N_s`$ values of $`B_i`$, are chosen at random in the gaussian law with a standard deviation equals 2. The value of $`B_i^{\mathrm{min}}=\mathrm{min}_jB_i^{(j)}`$ is specific for each amino acid and is the energy of the solute surrounded by the clatracle conformation. It gives an insight of the hydrophobicity of the monomer $`i`$.

On the other hand, one assumes that in pure solvent, the water molecules can evolve very freely and then, the $`N_s`$ conformational states of a bond between two sites of solvent molecules have the same energy (fig.2a), let says $`B_s`$. $`f_s(T)`$ can then be rewritten : $`f_s(T)=B_sT\mathrm{ln}N_s`$ where the $`B_s`$ parameter is related to the pH of solution as follow : when pH of the solution equals 7, the water molecules form a large number of hydrogen bonds, then $`B_s`$ is very negative. As pH of the solvent is lowered from 7, more and more hydrogen bonds between water molecules are broken and then the solvation energy $`B_s`$ increases but is still negative. Then, the larger $`B_s`$, the lower the pH of the solution.

Let us now show that, the total free energy of conformation can be rewritten as the sum of effective pairwise interactions between monomers. The total number of lattice links, $`n_{\mathrm{tot}}`$ can be written as the sum of number of covalent bonds, intrachain contacts, monomer-solvent sites interactions and solvent-solvent sites contacts : $`n_{\mathrm{tot}}=(N1)+_{i=1}^N_{j>i}^N\mathrm{\Delta }_{ij}^{(m)}+_{i=1}^Nn_i^{(m)}+n_s^{(m)}`$ Let us note, $`\alpha _i`$ is the number of first neighbors the $`i^{th}`$ monomer except the next and previous ones ($`\alpha _1=\alpha _N=3`$ and $`\alpha _i=2`$ for $`2iN1`$). One has : $`n_i^{(m)}=\alpha _i_{ji}\mathrm{\Delta }_{ij}^{(m)}`$ Then, following eqs. 4, 5 and the free energy of conformation reads :

$$F_{\mathrm{tot}}^{(m)}(B_s,T)=F_{\mathrm{tot}}^{(ext)}(B_s,T)+\underset{i=1}{\overset{N}{}}\underset{ji}{\overset{N}{}}\left(\frac{1}{2}B_{ij}f_i(T)+\frac{1}{2}f_s(B_s,T)\right)\mathrm{\Delta }_{ij}^{(m)}$$
(7)
where $`F_{\mathrm{tot}}^{(ext)}(B_s,T)=(n_{\mathrm{tot}}3N1)f_s(T)+_i\alpha _if_i(T)`$ is the free energy of the extended chains. For these chains, the intrachain energy vanishes and the total free energy is only the contribution of the free energy of solvation. As, $`_i_{ji}f_i(T)\mathrm{\Delta }_{ij}^{(m)}=_i_{ji}\frac{1}{2}(f_i(T)+f_j(T))\mathrm{\Delta }_{ij}^{(m)}`$, eq.7 can be rewritten :

$$F_{\mathrm{tot}}^{(m)}(B_s,T)=F_{\mathrm{tot}}^{(ext)}(B_s,T)+\underset{i=1}{\overset{N}{}}\underset{j>i}{\overset{N}{}}B_{ij}^{\mathrm{eff}}(B_s,T)\mathrm{\Delta }_{ij}^{(m)}$$
(8)
where $`B_{ij}^{\mathrm{eff}}(B_s,T)`$ are the effective couplings which now, dependent on the temperature :

$$B_{ij}^{\mathrm{eff}}(B_s,T)=B_{ij}f_i(T)f_j(T)+f_s(B_s,T)$$
(9)
One must note that the free energy of conformation is rewritten as the sum of temperature-dependent pairwise interactions between monomers.

This temperature dependence of the couplings, that we study now, induces the cold denaturation phenomenon. In eq.13, the solvation contribution of the $`i^{\mathrm{th}}`$ monomer to the couplings $`B_{ij}^{\mathrm{eff}}(B_s,T)`$ is $`\delta _i(T)=f_i(T)+f_s(B_s,T)/2`$. The derivatives of $`\delta _i(T)`$ reads, $`\delta _i(T)/T=s_is_s/2`$ where $`s_s=\mathrm{ln}N_s`$ is the entropy of the group of neat solvent and $`s_i(T)`$ is the entropy of the solvent around the solute $`i`$. One has $`s_i(0)=0`$ and $`s_i(T)\underset{T\mathrm{}}{}\mathrm{ln}N_s`$. It comes : $`\delta _i(0)/T=\mathrm{ln}(N/2)`$ and $`\delta _i(T)/T\underset{T\mathrm{}}{}\mathrm{ln}(N/2)`$. Then, there is a value of $`T_i`$ for which $`\delta _i(T_i)/T=0`$ and for which $`\delta _i(T)`$ has a maximum value. And, as $`\delta _i(0)=B_i^{\mathrm{min}}+B_s/2`$, for some values of $`B_s`$, the contribution $`\delta _i(T)`$ is positive for small and large value of $`T`$, inducing on averaged, repulsive interactions of monomer $`i`$ and, in contrast, $`\delta _i(T)`$ has negative values for medium values of the temperature and compact structures are then favored.

For each value of ($`B_s,T`$) an order parameter, is computed and define the state of the chain :

$$Q(B_s,T)=\frac{1}{N_{\mathrm{max}}}\underset{m_1\mathrm{\Omega }}{}\underset{m_2\mathrm{\Omega }}{}N_{\mathrm{com}}^{(m_1m_2)}P_{eq}^{(m_1)}P_{eq}^{(m_2)}$$
(10)
where the two summations are taken over all the chain conformations and $`N_{\mathrm{com}}^{(m_1m_2)}`$ is the number of common contacts between the conformations $`m_1`$ and $`m_2`$, $`N_{\mathrm{max}}`$ is the number of contacts of the more maximally compact conformation (here $`N_{\mathrm{max}}=9`$) and $`P_{eq}^{(m)}(B_s,T)=\mathrm{exp}[F_{\mathrm{tot}}^{(m)}(B_s,T)/T]/Z(B_s,T))`$ is the equilibrium probability of occurence of the conformation $`m`$. Hence, when only one structure $`m`$ is relevant, one has $`P_{\mathrm{eq}}^{(m)}=1`$ and then $`Q(B_s,T)=N^{(m)}/N_{\mathrm{max}}`$, with $`N^{(m)}`$ is the number of contacts of the conformation $`m`$. If the conformation $`m`$ is more maximally compact, one has $`Q(B_s,T)=1`$. On the other hand, if a large number of structure have a non-null probability of occurence, one has : $`Q(B_s,T)0`$.

The $`B_{ij}`$ values are those of the sequence R given in  after centering the mean of these values on 0. In order to fully solvate the more extended chain, one supposes that at least $`2N+2`$ groups of solvent molecules surround the polymer. The thermal averaged of the quantities of interest are computed for each couple of values of $`B_s`$ and $`T`$ such as $`B_s`$ vary from -5.00 to -9.00 with a step of 0.05 and $`T`$ varying from 0.00 to 2.00 with a step of 0.02. One generates sets of $`N_s=10000`$ values of $`B_i^{(j)}`$ for each amino acid to compute the partial free energies of solvation $`f_i(T)`$, following eq.6. Phase diagrams computed with $`N_s=2000,5000,10000,20000`$ remain qualitatively unchanged. The boundary between the two regions is only shifted along the $`B_s`$ axis.

In conclusion, we have presented a new model of solvation based on thermodynamic and structural properties of the aqueous solvent around protein. We showed, that the interactions involving in this solvation model can be written as temperature-dependent pairwise interaction between monomers including, the temperature and a solvent parameter which can be related to the pH of the solution. The model has been applied to a two-dimensional chain and we obtained a phase diagram that exhibits warm and cold denaturations and which is qualitatively in good agreement with experimental results. It must also be noted that our solvation model is easily applicable to longer chains on a three-dimensional lattice.

Acknowledgments to Dragi Karevski for helpful discussions to Enrico Carlon and Christophe Chipot for critical reading of the manuscript and to Centre Charles Hermite, Vandoeuvre-lès-Nancy, France for allocation of computer time.