Through the human genome project we now know that a human cell can synthesize about 20,000 to 25,000 different proteins (1). Proteins are an important class of biological macromolecules present in all biological organisms, and constitute high proportion of the dry mass of all cells (2). Most of the biological processes in all cells are executed by proteins. The amino acid sequence of a protein contains all information needed for adopting its three-dimensional structure. However, misfolding does occur, even though help from other molecules, such as chaperons, for correct and fast in vivo folding are in place (3–5).

Denaturation studies are very useful for investigating the thermodynamic properties of proteins. Transition from native to denatured states can be brought about by changing the properties of protein's environment. In general, this can be done by increasing the temperature, adding chemical denaturants or changing the pH.

Urea and guanidinium ion (used in the form of guanidinium chloride-GdnHCI) favor the denatured state by increasing the solubility of the unfolded chain in an aqueous solution. In comparison to temperature denaturation, chemical denaturation is often a reversible process. This is possible since the hydrophobic groups of the unfolded chain are shielded by the denaturants, which prevent aggregation.

The unfolding free energy (ΔG _U ) depends linearly on the denaturant concentration as: Eq.1 Δ G U = Δ G U H 2 O - m U [ Denaturant ]

Where is Δ G U = Δ G U H 2 O the free energy of unfolding in the absence of denaturant and m _U denotes the dependency of free energy on denaturant concentration (i.e. m _eq ) (6). A good linearity is observed at high denaturant concentrations and Δ G U = Δ G U H 2 O is obtained by extrapolation to the zero concentration of denaturant. Δ G U = Δ G U H 2 O values calculated from guanidinium chloride and urea denaturation are in very good agreement (7) which gives this relation some further credibility.

One of the major challenges in the field of protein science is to predict the stability and function of proteins from their primary structures. To accomplish this task, efficient algorithms are needed to relate the structure to stability. The availability of about 77,000 protein structures in Protein Data Bank (PDB) (8) and a great deal of experimental works on the thermodynamic stability of proteins have provided a wealth of information which can be used for the development of empirical functions that relate thermodynamic and structural parameters.

The success of such approach in developing structure-based methods to predict various thermodynamic parameters that define the Gibbs energy, i.e., the enthalpy, entropy and heat capacity changes, has been shown previously (9–13). In the process of unfolding, the major contribution to the enthalpy change arises from the disruption of intramolecular interactions such as van der Waals and hydrogen bonds and also solvation of the interacting groups. Therefore, the change in solvent accessible surface area (ΔSASA) upon unfolding has been used as a mean for predicting the ΔH as presented below: Eq.2 Δ H = ∑ i α i ( ρ ¯ ) ⋅ Δ SAS A i

Where ΔSASA _i is the change in SASA of atom i upon unfolding, and α i ( ρ ¯ ) is a coefficient that depends on the atom type and the average packing density of that atom within the protein (14).

The heat capacity change (ΔC _p ) in protein unfolding largely arises from changes in the hydration of groups that are buried in the native form away from the surrounding aqueous environment. ΔCp is correlated to the changes in SASA upon unfolding, as shown in the following equation: Eq.3 Δ C p = ∑ i a i ⋅ Δ SAS A i

Where a _i is the contribution of atom i per unit area and ΔSASA _i is as defined above. Using both equations, good correlations were obtained between experimental and calculated ΔH and ΔCp values (14).

The aim of current study is to develop empirical models to account for the effect of crosslinks on ΔSASA and hence on thermodynamic parameters (i.e., m _eq and ΔCp) of protein unfolding based on computational approach.

The changes in solvent accessible surface area (ΔSASA) upon unfolding, as determined by the differences in solvent accessibilities of native form (calculated from the crystal structure) and denatured form (modeled by an extended polypeptide chain) are given for a set of 45 proteins in Table 1. The table also shows m _eq values from denaturation experiments, ΔCp of unfolding, number of residues as well as crosslinks present in each of these proteins taken from the compilation made by Myers et al (10).

Figures 1A and 1B demonstrate dependencies that exist between the denaturants m _eq values and the changes in the solvent accessible surface area upon unfolding. There are significant linear correlations in both cases, with the correlation coefficient (R) values of 0.85 and 0.87 for GdnHCl and urea, respectively. The slopes of the linear regression lines are 0.25 and 0.17 cal/ (mol.M.Å ² ) for GdnHCl and urea, respectively, indicating the stronger denaturing effects of GdnHCl.

Denaturation heat capacity changes (ΔCp) were also correlated with the ΔSASA strongly with the correlation coefficient of 0.97 as shown in Figure 1C. The same linear correlations between m _eq values and ΔSASA have been shown previously by Myers et al (10). ΔSASA has been also related linearly to ΔCp by others (20, 21).

The main purpose of this study is to reevaluate the effect of crosslinks on ΔSASA and also predict the m _eq and ΔCp of unfolding based on protein sequence information. These latter two parameters are amongst the important criterion indicative of the stability of proteins. Therefore, prediction or any improvement in the prediction of these values has significant theoretical and practical applications.

The presence of crosslinks such as disulfide bonds and heme groups in a protein (as shown in Table 2) will result in a more compact unfolded state, thus reducing the solvent accessibility of the unfolded polypeptide chain. To compensate for the effects of crosslinks, Myers et al (10) employed the results of different empirical methods (22) to estimate the magnitude of the reduction of solvent accessible surface area (ΔSASA) per disulfide bond. The reduction of ΔSASA per crosslink was estimated to be about 900Å ² .

In the current study, to find out more about the effect of crosslinking through theoretical and computational methods, the different unfolded models were generated for crosslink-containing proteins while the crosslinks were preserved or removed in the unfolded states generated by instantaneous unfolding method based on assigning standard bond length and angle values. Then the SASA values were calculated for the generated unfolded structural models (Table 2).

To quantitatively indicate the effect of crosslinks on ΔSASA upon unfolding a new term called Crosslinking Factor (CLF) was introduced (CLF was described in Materials and Methods section.) Effectively, CLF is a measure of reduction in the SASA of unfolded protein as a consequence of presence of a single crosslink, such as disulfide bond, and calculated to be equal to 918.5 Å ² . This value is the average of crosslinking numbers calculated for 16 crosslink-containing proteins listed in Table 2 for which the m _eq and ΔCp values were available.

In five proteins listed in the table, cross-links are formed via ligation of central ion atom of heme groups by sulfur or nitrogen atoms of the side chains of the interacting residues. The average of crosslinking numbers for these proteins (759.0 Å ² ) is smaller than the average of the numbers (991.0 Å ² ) for the remaining proteins where the crosslinks are formed by disulfide bounds. However, the difference is not statistically significant (p-value>0.05). None of these values are statistically different from the calculated CLF value of 918.5.

Based on the above findings, the ΔSASA values were corrected for the effect of cross-links on the solvent accessibility of the unfolded state by taking 918.5 Å ² per crosslink off the ΔSASA (called ΔSASA _corrected) and then the corrected values were re-correlated to the m _eq and ΔCp values. Linear correlation coefficients improved to 0.90, 0.88 and 0.99 for GdnHCL and urea m _eq as well as ΔCp values, respectively as shown in Figure 2.

The extent of increase in SASA upon unfolding of a protein highly depends on the number of residues (i.e. protein size) and the constraints present in the unfolded state. The unfolded state of a protein is populated by an ensemble consisting huge number of conformationally distinct species. The presence of structural constraints limits the conformational space available to be explored by the protein polypeptide chain. Our analyses, in agreement with the results of others (10), show that the amount of area buried in each protein correlates very strongly (R=0.99) with the number of residues in each protein (Eq. 9). The strong correlation between ΔSASA and the number of residues, makes it possible to estimate the thermodynamic parameters using equations 10 to 12.

Where k denotes the number of residues for a given protein. These equations provide means to predict m _eq and ΔCp directly based on the primary structure information. The results of experimental studies are in close agreement with the results of our theoretical calculations which indicate the important thermodynamic parameters can be predicted using ΔSASA upon unfolding and taking into account the presence of crosslinks in the protein.

In a different approach to estimate SASA of unfolded state of proteins, we have used MD to simulate the unfolding behavior of proteins in denaturing condition, as stated in Materials and Methods section. Four of the proteins in our dataset (listed in Table 3) were subjected to MD simulations for 10 ns at 500 °K while inserted in a solvation box filled by a mixture of water and urea molecules. As can be seen from the table, the maximum SASA values for the unfolded conformations of proteins obtained by MD are smaller than that achieved by non-simulation method. Consequently, the ΔSASA values are also relatively smaller. Eq.9 Δ SASA = 90.01 × k - 804.55 , R = 0.99 ; N = 45 Eq.10 m eq ( GdnHCl ) = 22.55 × k - 0.362 ( n × CLF ) + 835.25 , R = 0.89 ; N = 34 Eq.11 m e q ( U r e a ) = 17.55 × k - 0.173 ( n × C L F ) - 517.92 , R = 0.92 , N = 34 Eq.12 Δ C p = 18.48 × k - 0.163 ( n × C L F ) - 327.208 , R = 0.99 , N = 25

Figure 3 shows the snapshots of conformational changes during unfolding simulation of IgG binding domain of protein G (IBPG) which has 56 residues with no crosslink. As time evolves, both tertiary and secondary structures of IBPG are lost and at the same time its SASA increases. The maximum SASA achieved during 10 ns is 5772 Å ² which is less than that estimated for fully extended conformation (8143 Å ² ).

The presence of crosslinks in the unfolded state will result in a more compact unfolded form and the higher the number of crosslinks, the more pronounced is this effect. For example, as shown in Figure 4, the unfolded conformation of lysozyme (hen egg white), a 129-residue protein with four disulfide bonds, retained more globular shape at the end of MD simulation, although it loses the elements of secondary structures.

As shown in Table 3, the SASAs of the investigated proteins increased at the end of MD simulation. However, the extent of this increase is bigger for the protein with no crosslink. For example 1PGB which is a 56-residue protein without any crosslink showed 54% increase in SASA upon unfolding using MD method. However, applying the same unfolding condition on 1AAL, a protein with almost equal size (i.e. 58 residues) and two disulfide bond has led to only 40% increase in SASA. In all studied cases, the maximum SASA for unfolded conformations achieved by MD are smaller than that of instantaneous method.

Analyses of MD trajectories showed that the RMSD differences for C atoms increases as time evolves approaching high values in the range of ∼14–19 Å for the studied proteins during the simulation. The rate of RMSD increase was dramatically fast for 1PGB and 2TRX, with no and one crosslink, respectively. However, the rate was gradual in the case of 1AAl and 1AKI with two and four crosslinks, respectively.

Although MD simulation under the condition used in this study can unfold the proteins and also demonstrates the effect of crosslink, however, using this method the SASA of unfolded conformations never reached to the SASA values of the unfolded conformations obtained by instantaneously decomposing the protein native structure just by taking into consideration to preserve the standard bond lengths, bond angles and other standard chemical structure geometries. This could be due to insufficient simulation time or entrapment of protein in an ensemble of conformations in a local minimum of energy landscape. However, to find out more about these issues and draw more sensible conclusion, further computational experiments such as hydrodynamic simulation are required.

Crosslinks such as disulfide bonds and heme groups have profound effect on the conformational flexibility and SASA values of unfolded state and hence influence the stability of the proteins by decreasing the entropy gain as well as reduction of ΔSASA upon unfolding. Studies of proteins with chemical crosslinks have shown clearly that the major effect of the crosslink on the stability results from a decrease in the conformational entropy of the unfolded molecule (23, 24).

On the other hand, attempts to increase the stability of proteins through introducing disulfide bonds suggest that the structural restraints in the native state due to the crosslinks may also make an important contribution to the net effect of the crosslinks on the stability (25). Furthermore, inspection of the model structures of Micro-myoglobin (Mb) revealed a role for heme in stabilizing the folded state (26).

Doig and Williams (27) investigated the effect of disulfide crosslinks on hydrophobicity derived stability of proteins. Based on data obtained from solvent transfer experiment, they calculated the non-polar ΔSASA to be 590 and 690 Å ² per disulfide bond according to free energy of hydration and ΔC _p measurements, respectively. Taking into account that the fraction of total area buried which is non-polar is about 0.70, these values correspond to a reduction in the total area change of 850 Å ² and 990 Å ² per disulfide. Using solvent perturbation difference spectroscopy, Pace et al demonstrated that the solvent accessibility of the aromatic residues (Tyr and Trp) in three studied proteins (lysozyme, RNase A and RNase T1) was changed upon unfolding (22). Myers et al, used these experimental data to estimate an approximate average value of 900 Å ² reduction in ΔSASA _unfolded per disulfide, assuming a universal change in accessibility across all residue types (10).

The results of these experimental methods have been averaged and used by Myers et al to compensate for the effects of crosslinks on ΔSASA. However, it may suffer from an over simplification by using only the changes in accessibility of just two amino acids and extrapolating these changes to the total area. It should also be mentioned that these results have been concluded from very limited number of experiments performed on only three globular proteins (22). Eq.13 m e q ( G d n H C l ) = 835.25 - 325.49 × n + 22.55 × k R = 0.889 ; N = 34 Eq.14 m e q ( U r e a ) = - 517.92 - 155.85 × n + 17.55 × k R = 0.878 , N = 34 Eq.15 Δ C p = - 327.21 - 149.49 × n + 18.48 × k R = 0.990 , N = 25

One of the shortcomings of using either CLF, introduced in this work, or experimentally derived value of 900, introduced by Myers et al to compensate for the effects of crosslinks on ΔSASA and hence estimation of m _eq and ΔCp values is the scarcity of the data used. The correction value close to 900 Å ² (proposed by Myers and here as CLF) can be justified by fitting equations 12 to 14 presented in Myers et al where the disulfide bond corrections that maximize the fits are all close to 900 Å ² . However, there is no need to use a correction factors such as 900 Å ² proposed by Myers or CLF to account for the effects of crosslinks on m _eq or ΔCp. Although we believe the correction factor most likely is close to 900 Å ² , but it is not a magic number and any other value close to that can be used to do the correction and then draw empirical equations to relate m _eq or ΔCp to the corrected ΔSASA (or to the combination of number of amino acids and CLF as we used in here).

The coefficients in the final mathematical equations will be adjusted to balance out any changes in the value of correction factor. In a situation where the ultimate aim is to be able to predict the thermodynamic parameters as precise as possible, one may decide to use different structural descriptors to derive empirical equations for the prediction purposes.

To this end we have furthered our investigation by trying to develop different empirical equations to predict m _eq and ΔCp values. We have examined the effects of different structural properties such as number of amino acids, number of crosslinks, size of loops representing the total number of amino acids involved in the loops formed by crosslinks, and the central position of the regions in the loop area on the prediction of the thermo-dynamic properties. The best statistical Multiple Linear Regression (MLR) models were achieved using variables representing total number of amino acids and the number of crosslinks.

In order to test the predictive power of these models, the Leave One Out (LOO) cross validation method was used. The mean absolute percentage error of predictions (MAPE) of m _eq (GdnHcl), m _eq (Urea) and ΔCp values for all proteins listed in Table 1 (or proteins with crosslinks listed in Table 4) based on Myers’ models are 19.7 (21.5), 22.9 (31.8) and 13.8 (17.9). The corresponding MAPEs using models presented in equations 13, 14 and 15 are 22.3 (20.5), 23.6 (31.1) and 13.5 (16.8), respectively.

The results show that both methods are not statistically different in predicting the evaluated thermodynamic parameters, either for all data points (proteins in Table 1) or for the proteins with crosslinks (i.e. values indicated inside the brackets), and simple MLR equations based on limited number of structural descriptors, i.e. number of residues and number of crosslinks, are able to perform equally well. In fact equations 13 to 15 are identical to equations 10 to 12 and the only difference is the way to represent the effect of crosslink on the parameter of interest. For example in equation 13 the coefficient of variable n (i.e. number of crosslinks) equals to the coefficient of the second term on the right hand side of the equation 10 multiplied by the value of CLF. These equations are also too close to the equations proposed by Myers et al (eqs. 12 to 14 in reference 10). For instance, the coefficient of variable n in equation 14 above (i.e. 155.58) is very close to the 139.3 calculated by multiplying 0.14 and 995 in equation 13 in Myers’ study (10).

In summary, it can be concluded that the proposed relationships represent valuable tools for predicting thermodynamic parameters of protein folding using the primary sequence information. The proposed crosslinking factor (CLF; which shows the effect of a single crosslink on ΔSASA upon unfolding) of 918.5 Å ² obtained based on computational simulation is very close to the previously published experimentally derived value of 900 Å ² . Such a correction factor can be used to estimate the ΔSASA upon unfolding which in turn can be used for the prediction of thermodynamic parameters such as m _eq and ΔCp. For the prediction of these parameters, one may also use number of amino acids (k) and number of crosslinks (n) without need to any kind of correction factor.

Although the correction factor for the effect of crosslink on ΔSASA is a quantitative value describing a fundamental property in protein folding, however, for the prediction purposes, the use of more simple properties taken from the primary structure of proteins gives as well accurate results. In addition, the current work demonstrates an example where theory is capable of reproducing the results obtained from experimental works.