Predicting protein thermal stability changes upon point mutations using statistical potentials: Introducing HoTMuSiC

The accurate prediction of the impact of an amino acid substitution on the thermal stability of a protein is a central issue in protein science, and is of key relevance for the rational optimization of various bioprocesses that use enzymes in unusual conditions. Here we present one of the first computational tools to predict the change in melting temperature ΔTm upon point mutations, given the protein structure and, when available, the melting temperature Tm of the wild-type protein. The key ingredients of our model structure are standard and temperature-dependent statistical potentials, which are combined with the help of an artificial neural network. The model structure was chosen on the basis of a detailed thermodynamic analysis of the system. The parameters of the model were identified on a set of more than 1,600 mutations with experimentally measured ΔTm. The performance of our method was tested using a strict 5-fold cross-validation procedure, and was found to be significantly superior to that of competing methods. We obtained a root mean square deviation between predicted and experimental ΔTm values of 4.2 °C that reduces to 2.9 °C when ten percent outliers are removed. A webserver-based tool is freely available for non-commercial use at soft.dezyme.com.

energy estimations as well as structural, sequence and dynamical features, have been developed to get a prediction of the thermodynamic stability changes upon point mutations defined through ΔΔG(T r ), the difference of folding free energy ΔG between the mutant and wild-type proteins at room temperature. Quite interestingly, some of these methods can reach a good accuracy at very low computational cost, with the sole knowledge of the wild-type structure [13][14][15][16][17][18][19][20][21][22][23] (see also 24 for a comparison of their performances). This makes the fastest among them ideal tools for stability analyses of the entire structurome. One of the major drawbacks of these methods is that the results are usually biased towards the training datasets even if strict cross validation is applied, which makes the estimation of their true performances an almost impossible task 25 .
The impact of point mutations on the thermal stability, defined through ΔT m which measures how the protein melting temperature changes upon mutations, has been much less investigated than the thermodynamic stability, as it is more intricate and requires taking into account that the amino acid interactions are temperature dependent. Therefore, there are very few in silico tools for predicting ΔT m 16,[26][27][28] . The common strategy to study the enhancement of thermal resistance is to assume the thermodynamics and thermal stabilities to be perfectly correlated (or ΔΔG(T r ) and ΔT m to be perfectly anticorrelated). Unfortunately, even if this approximation can be used in a first instance, it is not always reliable 29 . As a consequence, the predictions of ΔT m are in general less accurate because the intrinsic errors on the ΔΔG predictions have to be summed with the errors due to the imperfect correlation of the two stabilities. As an example, the anticorrelation between ΔΔG's predicted using the thermodynamic stability change predictor PoPMuSiC 22 and measured ΔT m 's is not so satisfactory and is equal to 0.51, whereas the correlation between predicted and experimental ΔΔG's is 0.63.
For all these reasons, it is necessary to design a specific computational tool for predicting ΔT m in a fast and more precise way. This is the aim of the present paper, in which we present a new, knowledge-and thermodynamics-based, method called HoTMuSiC, which is able to predict this quantity using as sole input data, the three-dimensional (3D) structure of the wild-type protein and -when available -its melting temperature T m . A very preliminary version of this method has been published in 30 . The main reasons of the success of HoTMuSiC are rooted on the one hand in the thorough physical analysis of the system which helped correct guessing the form of the model structures, and on the other hand in the use of temperature-dependent statistical potentials 22,[31][32][33] that are extracted from non-redundant datasets of protein X-ray structures of thermostable and mesostable proteins. We would like to emphasize that these are presently the only potentials that yield an estimation of the temperature dependence of the folding free energy contributions of the different amino acid interactions -albeit in an approximate, effective, manner.

Results
Theoretical analysis. The thermodynamic stability change upon mutation is measured by ΔΔG(T r ), i.e. the difference between the standard Gibbs folding free energies of the mutant (ΔG mut ) and wild-type (ΔG wild ) proteins at the reference temperature T r : Usually T r is taken as the room temperature: T r = 298 K. ΔΔG's upon point mutations can be predicted in silico using a series of tools [13][14][15][16][17][18][19][20][21][22][23] , which reach a relatively good accuracy with a standard deviation between the experimental and predicted ΔΔG's of 1 to 2 kcal/mol. Less prediction methods have been developed for the change in thermal stability upon mutations, measured by ΔT m , i.e. the difference between the melting temperature of the mutant (T m mut ) and wild-type (T m wild ) proteins: In a first approximation the two protein stabilities can be assumed to be strongly interdependent. Indeed, focusing on two-state folding transitions and assuming: (1) the mutations to be small perturbations with respect to the wild-type state; (2) the folding heat capacity ΔCP to be T-independent; (3) its variations upon mutations to vanish (ΔΔC p ∫ ΔC p mut -(C p wild )=0), and (4) similarly for the folding enthalpy ), one can derive the simple relation: r . This equation directly relates the two stabilities, however with a coefficient that depends on the thermal characteristics of the wild-type protein through its melting temperature T m wild and its folding free enthalpy ΔH m . Note that, since the enthalpy change upon folding is negative, ΔΔG and ΔT m are anticorrelated, as expected.
Unfortunately, the situation is in general less simple, especially for highly destabilizing or highly stabilizing mutations, and we have to take into account that the enthalpy and heat capacity variations do not vanish, i.e. ΔΔH m ≠ 0 ≠ ΔΔC P . This is illustrated by the fact that the correlation coefficient between experimental ΔT m 's and ΔΔG's is about − 0.7, which signals an imperfect correlation between these quantities. In this case, and still assuming ΔC P to be T-independent, Eq. (3) becomes: Scientific RepoRts | 6:23257 | DOI: 10.1038/srep23257 Thus the simple relation between the two descriptors of protein stability is lost: the proportionality coefficient can be positive or negative, and becomes also mutation-dependent in addition of being protein-dependent. It is easier to understand the meaning of the possible types of correlations between the two stabilities with a graphical representation. If the assumption of small perturbation and as a consequence Eq. (3) holds, the full protein stability curve changes upon mutation like in Fig. 1a. Typically, in this case, the wild-type has one interaction more than the mutant (or conversely). Otherwise the scenario is more similar to Fig. 1b, where one cannot say a priori which type of connection there is between ΔΔG(T r ) and ΔT m without the knowledge of additional information. It can for example occur when an interaction that is more stabilizing at room temperature is replaced by another that is more T-resistant thereby modifying the ΔH m , or when a change in the 3D structure occurs which modifies the protein's solvent accessible surface area and thus the ΔC P .
In a first approach to the prediction of ΔT m upon mutations, we consider the small perturbation approximation and thus Eq. (3) as valid and compute ΔT m as the sum of usual, temperature-independent, statistical potential contributions. In a second approach, which is expected to be more accurate for highly destabilizing or stabilizing mutations, we use ΔΔG-values calculated at different temperatures. Note however that structural modifications are more likely to occur in such case, and thus that the accuracy of the energy evaluations on the basis of the wild-type structure alone is questionable. For the purpose of estimating ΔΔG(T), we use the formalism of the temperature-dependent statistical potentials introduced in 31-33 . Construction of the model. Standard and temperature-dependent statistical potentials. The standard formalism of statistical potentials [34][35][36][37] has been fruitfully applied to a variety of analyses that range from the prediction of protein structure, stability, and protein-protein and protein-ligand binding affinities. It basically consists in deriving a potential of mean force (PMF) from frequencies of associations of structure and sequence elements in a dataset of protein X-ray structures. Under some approximations whose validity has been discussed [38][39][40] and making use of the Boltzmann law, the simplest PMF can be written as: where c and s are structure and sequence elements respectively, F represent the relative frequencies of c and/or s which are expressed as a function of the number of occurrences n, k is the Boltzmann constant and T the absolute temperature. Following 41 , higher order potentials can be constructed by considering more than two structure elements and/or sequence elements. Considering for example two sequence elements s and s′ and one structure element c, the above expression of PMF gets modified as:  2 Using Eqs (5-6) and their generalizations defined in 41 , we derived 9 different statistical potentials from a dataset of about 4,100 proteins with well-resolved 3D structure and low sequence similarity; they are listed in Table 1. They differ by the number of sequence and/or structure elements involved and by their type. Each sequence element s is an amino acid type at a given position and each structure element c is either the spatial distance d between two residues, the backbone torsion angle domain t or the solvent accessibility a of a residue. The statistical potential formalism has recently been extended to include more properly the thermal characteristics of proteins and in particular the fact that amino acid interactions are temperature-dependent 9,31-33 . Following this approach, a dataset of about 170 proteins with known melting temperature and 3D structure was used and divided into two subsets, one with only mesostable proteins (with T m less than about 65 °C), and the other with thermostable proteins (with T m higher than about 65 °C). A series of 9 different potentials were extracted from each subset, which are listed in Table 1. To limit the number of parameters to be optimized (see next section) and thus to avoid overfitting as much as possible, these twelve potentials were combined into five potentials ( Table 1).
As expected, these T-dependent potentials reflect the thermal characteristics -mesostable or thermostableof the subset from which they are derived: the former set better describes the interactions in the low temperature region while the second set better reflects the thermal properties at high temperatures. This approach has shown good performances in the prediction of the thermal resistance and of the stability curve as a function of the temperature of proteins that belong to the same homologous family 32,33 .
The structure elements defining these potentials are the same as for the temperature-independent potentials. In contrast, the size of the mesostable and thermostable protein datasets is much smaller than the dataset used for the standard statistical potentials, as it requires the knowledge of the melting temperature. To deal with the smallness of the datasets, we used several tricks that consists of corrections for sparse data and the smoothing of the potentials, following 9,31-33 .
In addition to the standard and T-dependent potentials, we also considered volume terms in the folding free energy estimation, which are defined as the volume difference ΔV between the mutant and wild-type amino acids 22,23 . In order to take into account that the creation of a cavity in the protein interior (ΔV < 0) can have a different impact on the stability compared to the addition of stress (ΔV > 0), we introduced two separate terms (ΔV − ) and (ΔV + ) defined as where θ is the Heaviside function.
Artificial neural networks and parameter identification. The above-defined potential terms were combined to predict how the melting temperature changes upon mutations, using two different model structures. The second model (called T m -HoTMuSiC) uses the T m of the wild-type, while the first (HoTMuSiC) does not. To identify the parameters involved in these combinations, we used an artificial neural network (ANN) with peculiar activation functions.
In the first approach (HoTMuSiC), we assumed that ΔT m can be written as the sum of twelve contributions, the nine energy terms ∆∆ = ∆ − ∆ ν ν ν W W W wild mut (ν = 1, … 9) computed from the standard, T-independent, statistical potentials listed in Table 1, the two volume terms and an "independent" term that only depends on the solvent accessibility. The functional form reads as: where N r is the number of residues in the protein, and the coefficients α ν (A) were chosen to be sigmoid functions of the solvent accessibility A of the mutated residues: We have chosen the activation functions to be sigmoidal since they model a smooth transition from the protein core to the surface, and since the weight of the different energy contributions have been shown to differ in these two regions 22,42 .
To identify the fifty parameters introduced in Eq. (7), we have chosen a standard feedforward ANN with one input and one output layer (schematically depicted in Fig. 2a). The cost function to be minimized is the mean square deviation between the experimental and predicted values of ΔT m for the dataset S mut that contains N mut = 1,626 mutations for which an experimental ΔT m -value is available (see Methods section): The initial values of the weights were chosen randomly. To take into account that the ANN training algorithm can get stuck in local minima, the initialization and training processes were repeated about thirty times and the solution reaching the lowest σ-value was chosen 43,44 .
We used a strict five-fold cross validation procedure with an early stopping procedure for evaluating the method's performance. More precisely, the entire set of mutations was split randomly into a training set (containing 90% of the mutations) and a test set (with the remaining10%). The training set was then further randomly split into a smaller set (with 80% of the mutations) on which the parameters were identified, and a validation set (with 10% mutations) on which the early stop threshold was determined, namely the maximum number of iterations in the gradient descent procedure before its convergence 45,46 . This early stop is necessary to avoid overfitting, as in general the network starts to learn too much from the training dataset after a certain number of iterations, with the consequence that the error in cross validation starts to raise. As a final step in the computation, the prediction error σ is independently calculated on the test set.
In the second method for predicting ΔT m , called T m -HoTMuSiC, we used quite a different approach, with as building blocks the T-dependent statistical potentials listed in Table 1 and the T m wild of the wild-type protein. More precisely, Eq. (7) describing the ΔT m functional was modified into: wild are polynomial functions of the melting temperature of the wild-type protein and its number of residues N r . Their functional form is guessed from Eq. (3) and approximated as: To identify the 67 parameters of this second method, an ANN with an input layer, a hidden layer and an output layer is used, as shown schematically in Fig. 2b. The input layer consists of three sets of neurons, one set encoding the mesostable potentials, a second one the thermostable potentials, and a third one the volume and independent terms. Three perceptrons use these three sets of input neurons and generate the three output signals of the hidden layer. These are the input of yet another perceptron, which yields the ΔT m -prediction as output. The initialization and identification procedures of all the weights and the cross validation procedure are the same as for the first method.
The final ΔT m predictions of the T m -HoTMuSiC method are defined as the mean of the two predictions given by Eqs (10) and (7):  Table 1, two volume terms and an independent term; (b) T m -HoTMuSiC method: 3-layer ANN, consisting of 3 perceptrons with sigmoid weights; the first perceptron has 5 input neurons encoding the 5 mesostable potentials listed in Table 1, the second has 5 input neurons corresponding to the 5 thermostable potentials, and the last perceptron has 3 neurons for the volume and independent terms. The outputs of these three perceptrons (Meso, Thermo, and Vol + I) are the inputs of another perceptron with polynomial weight functions.

Performance of HoTMuSiC and T m -HoTMuSiC.
The values of the root mean square deviation σ between measured and predicted ΔT m values (Eq. 9), computed in strict cross validation, are shown in Table 2. For HoTMuSiC, we obtained σ = 4.3 °C; the Pearson correlation coefficient r between experimental and predicted ΔT m 's is 0.59. The performance of the T m -HoTMuSiC version is slightly better with σ = 4.2 °C and r = 0.61. When ten percent outliers are excluded, σ decreases to 2.9 °C and r rises to 0.75. The results are plotted in Fig. 3.
The T m -HoTMuSiC version that encodes information about the melting temperature of the wild-type protein and uses T-dependent potentials thus performs slightly better than the T-independent version, but not as much as could be expected on the basis of earlier analyses 32,33 . Indeed, in principle, the mesostable potentials should describe the mutations in the mesostable proteins much better than the thermostable potentials and vice-versa. The reasons for this mitigated result could be due to the lower accuracy of the T-dependent potentials compared to the standard ones since they are extracted from smaller protein sets. Or they could be rooted in the information loss upon the reduction of the number of potentials (see Eq. (7) versus Eq. (10)), which is done to avoid overfitting issues.
Moreover, we can see from Table 3 and Fig. 4 that the σ-value for amino acid substitutions in the core (A < 15%) and partially buried positions (15% < A < 50%) is on the average larger than that of surface mutations (A > 50%). In contrast, the correlation coefficient r is higher in the core and in the partially buried region and smaller at the surface. This apparent discrepancy is in fact due to the higher variance of ΔT m in the core, which  Table 2. Scores of HoTMuSiC and T m -HoTMuSiC; σ * and * r correspond to σ and r with 10% outliers removed.   Table 3. Scores of T m -HoTMuSiC as a function of the solvent accessibility A of the mutated residues. drives the correlation and increases the value of r. In the surface region, the predictions are more accurate (lower σ) but the ΔT m variance and the correlation coefficient are lower. When normalizing σ by the standard deviation of ΔT m , we obtain values that increase from the core to the surface (see Table 3).
Comparison with other ΔT m predictors. As far as we know, only two other ΔT m predictors are described in the literature, which are strongly different from ours. The strategy of Saraboji et al. 28 consists in predicting for a mutation from wild-type residue W to mutant M the mean value of the ΔT m 's of all the analogous mutations W → M in the training dataset. A similar strategy proposed in the same work is based on the classification of the mutations in terms of the secondary structure and solvent accessibility of the wild-type residues and predicts as ΔT m the mean of the experimental ΔT m 's occurring in the suitable class in the learning set. A limitation of this approach is that not all wild-type to mutant mutations are present in the learning set due to the lack of experimental data. The second method is called AutoMute 16,26,27 and is based on residue environment scores. It proceeds by reducing protein 3D structures to ensembles of C α atoms, and applying Delaunay tessellation to identify quadruplets of nearest neighbor residues. A log-likelihood potential is constructed from the number of occurrences of the quadruplets in a dataset of 3D structures and then used as the key ingredient in the computation of ΔT m .
To make a cross-validated comparison between our results and those of these two methods, we have chosen a subset S sub of our dataset S mut consisting of the mutations that are not present in the AutoMute learning set (see 48 for a list), and trained versions of T m -HoTMuSiC and the Saraboji method on the set S mut \S sub . The performances of the three methods are reported in Table 4. T m -HoTMuSiC shows the best performance, with an improvement of about 20% and 30% with respect to AutoMute 16 and the Saraboji method 28 , respectively. The σ-values computed on the S sub set are equal to 3.7, 4.7 and 5.4 °C for T m -HoTMuSiC, AutoMute and Saraboji method, respectively. Note that the performance of the latter two methods has been evaluated on a slightly reduced dataset, since the ΔT m values could not be computed for some of the mutations.

Discussion
We developed a thermodynamics-based and knowledge-driven ΔT m prediction method that does not exploit the common assumption of a perfect correlation between thermal and thermodynamics stabilities. The basic ingredients of our approach include standard and T-dependent statistical potentials that are combined through the use of ANN's. The performance in cross validation of the two versions of our method, HoTMuSiC and T m -HotMuSiC (which requires the T m of the wild type), are quite good with a σ-value of 4.3 and 4.2 °C, respectively, which goes down to 2.9 °C upon removal of 10% outliers. They perform significantly better, by 20 to 30%, than the two other ΔT m prediction methods described in the literature.
HoTMuSiC and T m -HotMuSiC are accessible via the webserver soft.dezyme.com and are free for non-commercial use. They are extremely fast and allow the ΔT m predictions for all possible single-site mutations in a protein in a few minutes. This webserver will be presented in an application note 49 .
Our software thus yields quite accurate results, and allows rapid screening of all possible point mutations in a protein structure and identifying a subset that is likely to yield the required thermal resistance. This subset must then be analyzed further, either by using more detailed computational techniques, or by experimental means. HoTMuSiC and T m -HotMuSiC are very useful and user-friendly tools for every researcher who wishes to rationally design modified proteins with controlled characteristics.
Notwithstanding the large applicability and good accuracy of HoTMuSiC and T m -HotMuSiC, it is worth discussing their limitations and the sources of errors that affect the predictions. These are: • The wild-type and mutant structures are supposed to be identical (up to the substituted side chain) and the possible structural modifications are only encoded in the volume terms ΔV ± ; local structure rearrangements upon residue substitutions, for example in the hydrophobic core, yet depend on many more parameters such as the residue depth and the backbone flexibility [50][51][52][53][54] . • The mutation dataset is strongly unbalanced towards destabilizing mutations, which is likely to add unwanted hidden biases even if a strict cross validation procedure is applied 25 . Only when more stabilizing mutations will be experimentally characterized will we be able to completely exclude the biasing impact of this stabilizing-destabilizing asymmetry. • The experimental conditions at which the T m measurements are performed usually differ in terms of pH, ion concentration or buffer composition, which induces noise in the learning set and errors in the predictions.
To limit this problem, we chose the entries derived from experiments performed at pH as close as possible to seven, and made a weighted average of the experimental ΔT m 's of a same mutation, when available. • The T-dependent potentials suffer from the smallness of the dataset of protein structures with known T m . Different tricks have been used to limit this issue.  Table 4. Comparison between the performances of HotMuSiC and those of the two other methods evaluated in cross-validation on the dataset S sub .
• The possible parameter overfitting is like always an important concern, especially for the T m -HoTMuSiC version, in which the number of potentials is three times larger than for HoTMuSiC. To avoid overfitting, we decided to decrease the number of parameters in T m -HoTMuSiC by fixing some coefficients of the linear combination of potentials (see Table 1).
Different ways will be explored in an attempt to further improve the prediction performances of HoTMuSiC. They obviously include the enlargement of the datasets of proteins of known structure and T m , and of the mutations of known ΔT m . We will also investigate different ANN architectures, the addition of hidden layers, and the inclusion of other features such as the change in conformational flexibility upon mutation, which seems related to the thermal stability even if a quantitative connection between the two quantities on a large scale is still missing [55][56][57][58] . Finally, it could be worth analyzing the wrong predictions in view of identifying the factors that should be taken into account to make HoTMuSiC even more performing.

Methods
Set of experimentally characterized mutations. We started collecting the mutations with experimentally measured ΔT m value from the ProTherm database 59 and the literature. Each entry was then manually checked from the original literature to remove errors and select those that satisfy the following criteria: were only considered (1) mutations in monomeric proteins of known X-ray structure with resolution below 2.5 Å, (2) mutations whose experimental ΔT m was measured in absence of chemical denaturants, (3) simple two-state (un) folding transitions, and (4) single point mutations. Destabilizing or stabilizing mutations by more than 20 °C were overlooked, as they probably induce important structural modifications that our method is unable to model. Using this procedure, we obtained a set S mut of 1,626 mutations that belong to 90 proteins and have an experimental ΔT m . More information and their list can be found in 48 .