Prediction of hot spot residues at protein-protein interfaces by combining machine learning and energy-based methods

Prediction of hot spot residues

We discuss first the binary classification problem, i.e. the problem of predicting if a residue is a hot spot (ΔΔG ≥ 2 kcal/mol) or not (ΔΔG < 2 kcal/mol). Support Vector Machines (SVMs) are widely used tools in computational biology and well suited for the task (see [20, 21] and references therein). As input features we consider the 8 terms associated with side-chain inter-molecular and environment inter-molecular energies (the addition of side-chain intra-molecular energies in this case does not improve the prediction accuracy). A summary of the results is reported in Table 1a according to various performance measures. The precision P is the fraction of true hot spots among the set of residues predicted to be hot spots; the recall R is the fraction of correctly identified hot spots relative to all those present in the data set; the F 1 score is a weighted average of the precision and recall; the Matthews Correlation Coefficient (MCC) is a commonly used measure of the quality of binary classifications (see Methods section for more details).

Results quoted in Table 1a have been obtained through a 16-fold cross-validation strategy. We have verified that they are robust with respect to the cross-validation scheme employed, i.e. by using a more stringent 12-fold cross-validation the results do no change appreciably (see Additional file 2: Supplemental Table S1). As can be expected they tend to be slightly worse but the difference does not appear to be significant. A further control is to apply the models on their training sets and check if statistically the predictions are much better than on the test set. We have verified that the performances on the training sets are only marginally more accurate but still comparable to those on the test set (see Additional file 2: Supplemental Table S1). These considerations about the 12-fold cross-validation and the performances on the training sets suggest that there is no over-fitting of the data in our analysis and that the SVM model generalises to unseen data as quoted in Table 1a.

One concern related to the data set is that it contains a large number of immunoglobin domains and this might introduce a bias in the predictions. Ideally, to verify the extent of this problem, one should group all the immunoglobin domains in the same fold family when performing cross-validation and check if this affects significantly the performances. In practice this is not feasible, as the data set is relatively small and clustering the immunoglobins together would result in one fold dominating the database. However a first indication that there is no such bias towards predicting mutations on immunoglobins more accurately comes from the results of the 12-fold cross-validation. In this case, most (but not all) of the immunoglobin containing complexes are clustered together and the results do not change significantly. A second, a posteriori check is to analyse classification predictions separately for complexes containing an immunoglobin domain and those without. If there were a bias, we would expect a better performance on the former compared to the latter class but instead we obtain P = 0.55 ± 0.03, R = 0.60 ± 0.05, F 1 = 0.57 ± 0.03, MCC = 0.40 ± 0.05 and P = 0.52 ± 0.03, R = 0.66 ± 0.02, F 1 = 0.58 ± 0.02, MCC = 0.46 ± 0.03 respectively.

The limited contributions from the electrostatic terms can be ascribed to their weak correlations to the target outputs, i.e. to the observed ΔΔG values (see Additional file 2: Supplemental Table S2). The non essentiality of the environment van der Waals and hydrogen bond potentials seems instead to derive from the fairly high correlation existing between these two terms and the environment desolvation potentials (see Additional file 2: Supplemental Table S2). It suggests that if one of these two terms is missing the latter can effectively substitute for it (possibly with some readjustment from the other remaining terms as well). It is important to underline that these results do not imply that the mentioned 4 terms play no role and can be altogether ignored. Indeed omission of pairs of features can lead to a significant decrease in the performance.

The importance of the side-chain van der Waals term agrees with the observation that hot spot atoms form good packing interactions with the partner proteins [23–25]. The side-chain hydrogen bond term was found to provide a major contribution also in [11]. The role of the desolvation potential seems to support the O-ring hypothesis and the importance of shielding the interactions from the solvent. Exclusion of the solvent leads to a lower effective dielectric thereby increasing the strength of an interaction. In this respect it is somewhat surprising that electrostatics does not emerge as a key component in our model. Although a similar result was found in ref [11], it is possible that a better description of electrostatic effects is required, for example either by solving the Poisson-Boltzmann equation or through the generalised Born model.

The above considerations indicate that although some terms are more important than others there is no single feature that makes a dominant contribution. Rather, it seems it is the balanced combination of terms in the SVM model that allows the detection of hot spots. This possibly provides a justification of why this prediction problem is hard. It also support the claim that there is no simple patterns of hydrophobicity, shape or charge that can be used to identify hot spots [7].

To further examine the reliability and usefulness of our approach, we have compared our predictions with the predictions of the Robetta server on the same set of mutations. The server returns an estimated value for ΔΔG based on an all-atom free energy function. It can then be turned into a binary classifier by labelling a residue as a predicted hot spot if ΔΔG _calc ≥ x _th kcal/mol, where x _th is some threshold value. In the following we denote with such a classifier. Given our definition of hot spots, a natural choice is x _th = 2. We have however also explored other values for x _th because Robetta in its original implementation defines hot spots using a threshold of 1 kcal/mol and is therefore not optimised for our hot spot definition. We find indeed that more accurate predictions are obtained by setting x _th = 1.8. Results are reported in Table 1a. By comparing them to those for our SVM approach, it can be deduced that the latter yields a substantial improvement. We emphasise though that the comparison is not entirely fair: the Robetta method has been designed to predict the actual value of ΔΔG and as we will discuss below this a considerably more difficult problem than binary classification.

As mentioned above, in the original Robetta paper [11] a threshold of 1 kcal/mol was used to define experimental and predicted hot spots. With this definition, our data set is more balanced with respect to the number of positive and negative examples, respectively 165 (hot spots) and 184 (non hot spots). We have trained and tested our method with this hot spot definition as well and we report a summary of the results in Table 1b. Interestingly, in this case we obtain the best predictions by including the side-chain intra-molecular energy terms as well. Compared to the Robetta server (i.e. Robetta₁ in our notation), our method appears to perform marginally better but the difference might not be statistically significant and it is certainly less pronounced than in the case of a 2 kcal/mol threshold (see Table 1).

Recently, a machine learning approach to predict hot spot residues has been presented [13]. It is based on decision trees and trained on features such as geometrical shape and biochemical properties (e.g. atomic contacts, hydrogen bonds and salt bridges). The predictive performance of the method has been estimated to be P = 0.49, R = 0.58 and F 1 = 0.53 (F 1 = 0.55 if the method is combined with the Robetta server). Comparing methods on the basis of quoted results is problematic as data sets and cross validation strategies differ. The data sets of alanine mutations used in [13] and in our investigation are not identical but they do overlap substantially (indeed very similar performance scores are obtained by applying Robetta₂ to the two data sets, see Table 3 in [13] and Table 1a in this paper). The higher performance scores we obtain might therefore reflect a genuine improvement in hot spots prediction accuracy. The reason behind this improvement possibly lies in the inclusion of the environment energy terms in our method. Indeed results for SVM-sc in Table 2 appear comparable to those reported in [13]. Besides SVMs, we have also tested Gaussian Processes (GP) models for the classification task but found results with overall lower accuracy (see Table 1a).

After our manuscript was submitted and while still under review, two new studies on the same problem of hot spot prediction at protein-protein interfaces have been published [26, 27]. It is not straightforward to assess how our method performs compared to them. We note for example that on data sets assembled from ASEdb the reported F1 score for the methods in [26] and [27] are respectively F 1 = 0.65 and F 1 = 0.57. However, on the same two data sets Robetta₂ achieves respectively F 1 = 0.55 and F 1 = 0.59, which are both substantially higher than the value it obtains on our data set (F 1 = 0.49). This suggests that our method and those in [26] and [27] can not be compared on the basis of the quoted results alone.

We have experimented with introducing unlabelled data in the training set, a problem that is often refereed to as semi-supervised learning. Unlabelled data can sometimes improve a classifier by providing a more reliable decision boundary, for example by requiring that it lies (in the input features space) in region of low density. Transductive SVM have been developed to work in a semi-supervised learning setting and take advantage of the information content in unannotated data. In our case unlabelled data correspond to interface residues for which ΔΔG is not known. As reported in Table 1a, we find TSVMs return marginally improved predictions although the difference might not be entirely significant. One possible explanation for this improved performance is that the inclusion of unlabelled data makes the training sets in the cross-validation procedia more balanced and statistically more uniform.

The most accurate predictions are obtained for mutations involving Lys. This is a positively charged, polar amino acid capable of forming up to 3 hydrogen bonds with the amino group at the end of its side-chain; its side chain contains 4 carbon atoms which confer to it a partial hydrophobic character. Lys is among the more frequent amino acids in hot spots (see Table 6). As an "extreme" example of a mutation involving Lys we show in Figure 2 the trypsin-trypsin inhibitor complex. The residue Lys15 from the inhibitor is a hot spot having the largest measured ΔΔG in our data set (ΔΔG = 10 kcal/mol). Our correct positive prediction is driven mainly by the side-chain van der Waals and hydrogen bond energies. Arginine too is a positively charged, polar amino acid and its side-chain can form up to 5 hydrogen bonds. Predictions however are not as accurate as for Lys with a fair number of false positives. By inspecting these false predictions individually, it emerges that often the positive scores are driven by a large favourable side chain hydrogen bond energy. It is possible that in fact this term is over-estimated as the loss of some of the hydrogen bonds upon alanine mutation might be (at least partially) mitigated by, e.g., the inclusion of water molecules at the interface.

Mutated amino acid	In data set		Hot spots (ΔΔG≥ 2 kcal/mol)			Enrichment in hot spots
	(Number)	(%)	(Number)	(%) (a)	(%) (b)
Arg	33	9.46	7	21.21	8.64	0.91
Asn	22	6.30	6	27.27	7.41	1.18
Asp	29	8.31	9	31.03	11.11	1.34
Cys	1	0.29	0	0.00	0.00	0.00
Gln	21	6.02	2	9.52	2.47	0.41
Glu	31	8.88	5	16.13	6.17	0.69
His	13	3.72	1	7.69	1.23	0.33
Ile	15	4.30	4	26.67	4.94	1.15
Leu	10	2.87	1	10.00	1.23	0.43
Lys	32	9.17	11	34.38	13.58	1.48
Met	2	0.57	0	0.00	0.00	0.00
Phe	11	3.15	2	18.18	2.47	0.78
Ser	28	8.02	1	3.57	1.23	0.15
Thr	24	6.88	1	4.17	1.23	0.18
Trp	23	6.59	9	39.13	11.11	1.69
Tyr	44	12.61	20	45.45	24.69	1.96
Val	10	2.87	2	20.00	2.47	0.86

The number and percentage of amino acids in our data set are shown. The number and percentage of hot spots for each amino acid type is also reported. For a given amino acid type, (%)^(a)is the percentage of hot spots with respect to the residues of that type in the data set (i.e. entry in column 4 divided by the corresponding entry in column 2); (%)^(b)is the percentage of hot spots of that type with respect to all hot spots in the data set (i.e. entry in column 4 divided by the sum of all entries in column 4). Enrichment in hot spots is calculated as the ratio of the frequency of a given residue type in hot spots (column 6) over the frequency of the same amino acid type in the whole data set (column 3). This table should be compared with Table 2 in [19]. Note that proline and glycine are not included in our data set.

Tryptophan and tyrosine are also common in hot spots [19] (see Table 6). They are both aromatic, hydrophobic amino acid (Trp is more hydrophobic than Tyr as the latter contains a hydroxyl group) and are both capable of forming one hydrogen bond (see Figure 3 for two examples of Tyr hot spots). In general, our model appear to predict fairly accurately Trp and Tyr hot spots (see Table 5). In this respect it is worth noting that 3 false positive and 3 false negative predictions respectively for Tyr and Trp are found in the same complex, formed by the Interferon-γ Receptor and Antibody A6 (pdb code 1JRH). As discussed in [28] this is a difficult case, with the individual monomers likely to undergo significant conformational changes upon binding. Indeed predictions for this complex are among the less accurate we obtain (see Additional file 2: Supplemental Table S4). It has previously been noted that tyrosine is much more likely to be found in hot spots than phenylalanine, despite the two residues being similar and with nearly identical volumes. Presumably, this is due to the ability of tyrosine (and not phenylalanine) to hydrogen bond. It is reassuring that our model is able to reproduce this empirical observation fairly well (note that the identity of the mutated amino acid is not an explicit input feature).

For some amino acid types, our model has some clear limitations. For example, it performs poorly on mutations involving Glu despite it does reasonably well on Asp mutations. It has been observed that Asp is found more frequently in hot spots than Glu and this might be related to differences in side-chain conformational entropy [19]. This possibly also explains the low accuracy of our predictions, given that our method does not consider any entropic term. Another difficult case is isoleucine. Ile is an aliphatic, hydrophobic amino acid and similarly to Leu our method predicts no hot spots for this residue. In fact whereas Leu is rarely found in hot spots, Ile is actually enriched. Our model fails to distinguish between these two residues. In theory, one could think of exploiting the amino acid identity as input feature too and of building a model for each different amino acid. In practice, at present this is not feasible as there are not enough mutational data available.

Prediction of ΔΔG values

We turn now to the more complicated regression problem, i.e. the prediction of the actual value of ΔΔG induced by an alanine substitution. We have found beneficial in this case to include the intra-molecular energy terms as well. The results are summarised in Table 7. Both SVM regression and GP models compare favourably to the the Robetta server, although the difference is probably only marginal. Figure 4 shows the scatter plot of the predicted versus the observed ΔΔG, both for SVM regression and for the Robetta server. SVM regression returns a lower root mean square error and an higher correlation coefficient. On the other hand, the regression line in Figure 4 is closer to the ideal case (i.e. slope and intercept equal to one and zero respectively) for the Robetta predictions. Notice that for some choices of hyper-parameters SVMs do return solutions with slope and intercept of the regression line closer to the Robetta ones. For these solutions the root mean square error and correlation coefficient are also similar to those derived from Robetta.

Method	RMSE	r	Precision	Recall	F1 score	MCC
SVM	1.34 ± 0.05	0.54 ± 0.03	0.62 ± 0.05	0.44 ± 0.05	0.52 ± 0.03	0.41 ± 0.04
GP	1.36	0.52	0.58	0.47	0.52	0.39
Robetta	1.52	0.47	0.52	0.47	0.49	0.35
LLSF	1.36	0.51	0.54	0.46	0.50	0.36
LLSF(a)	1.45	0.44	0.54	0.42	0.47	0.34
LLSF(b)	1.43	0.47	0.55	0.43	0.48	0.35

SVM: Support Vector Machine, GP: Gaussian Processes, LLSF: Linear Least Squares Fit, LLSF^(a): Linear Least Squares Fit with inter-molecular side-chain van der Waals as only input, LLSF^(b): Linear Least Squares Fit with inter-molecular side-chain van der Waals and hydrogen bond as inputs. RMSE stands for root mean squared error, r for correlation coefficient and MCC for Matthews correlation coefficient (see Methods section for definition of the performance measures).

All tested methods appear limited in the accuracy they can achieve. For example, if < ΔΔG > is the average outcome of an alanine mutation simply by predicting ΔΔG = < ΔΔG > independently of the input features, one would obtain a RMSE = 1.58 kcal/mol. Of course in this case the correlation coefficient would be zero but this simple consideration suggests that root mean square errors of the order of those both we and other methods obtain are not extraordinary. In fact our results are comparable to those from a model built with an ordinary linear least squares fit (LLSF) as reported in Table 7. Compared to LLSF, SVM regression might still be slightly more accurate possibly because it is less sensitive to the effect of outlying points.

We have investigated also "minimal" version of the LLSF models, based on only one or two energy terms as inputs. This analysis highlights which are the major contributions to the models and reveals that, similarly to the binary classification problem, the most important term is the packing of the mutated side-chain against the partner protein. A model based on just the side-chain inter-molecular van der Waals energy achieves results that are worse but still comparable to the full model (see Table 7). Adding the hydrogen bond term the accuracy improves only marginally. We suggest therefore that the simple side-chain van der Waals model could be used in future as a benchmark to gauge advancements in this regression problem.

Progress in this field will probably require the development of novel approaches. Based also on our results, it seems likely that it will involve modelling of both the mutated and unbound structures. This indeed might be the major limiting factor in our strategy. Some improvements might also be expected by a more accurate description of the different energy terms. For example a more rigorous treatment of solvent effects on electrostatics might be necessary. Additionally, terms can be introduced to explicitly account for the structural plasticity and adaptability of hot spot regions, e.g. by means of normal mode analysis or dynamical simulations. It is also possible that the addition of evolutionary information can lead to more reliable predictions. Indeed, in a recent investigation it has been highlighted that hot spots tend to establish conserved physico-chemical interactions across homologous interfaces [29]. However, in order to achieve major improvements it will be essential to rely on a larger data set of alanine mutations.

Data Sets

The data set consists of protein complexes whose structures have been solved by X-ray crystallography and for which alanine mutational data are available. Structures are obtained from the Protein Data Bank (PDB) [36]. Alanine mutation data are collected from the Alanine Scanning Energetics database (ASEdb) [37] and from previous publications [11, 13]. We have only considered protein-protein interactions involving an extended interface. Mutation data related to protein-peptide complexes have not been included. For each reported mutation we have verified the original reference and checked that the residue mapping from sequence to structure is consistent (few discrepancies have been found). Only mutations occurring at the complex interface have been retained. Interface residues are defined as those having at least one heavy atom within 5 Å of an heavy atom in the binding partner. Similarly, two residues across the interface are considered in contact if any of their heavy atoms are within 5 Å.

To ensure that the data set is sufficiently diverse and representative of protein-protein interfaces in general, we have analysed the complexes in terms of interacting domains and of location of the binding interface (i.e. which residues and residue contacts are at the interface). Domain structures can be classified according to CATH [38], which follows a hierarchical scheme. The first five levels in this hierarchy are Class, Architecture, Topology (fold), Homologous superfamily and Sequence family (clustered at a threshold of 35% identity). We have required that no two pairs of interacting domains have the same CATH numbers at the S-level. If two domain pairs have the same CATH numbers, we verify if they use the same or a different interface to interact. In the latter case, both structure are included in the data set (in practise this issue arises for only one pair of structures, PDB codes 3HFM and 1VFB).

The data set contains 20 protein complex structures. Following previous publications [19], we define hot spots as those alanine mutations for which ΔΔG ≥ 2 kcal/mol (ΔΔG is the change in binding free energy). In total the data set comprises 349 mutations, of which 81 correspond to hot spots. Some studies such as for example in the Robetta paper [11] define hot spots using instead a threshold of 1 kcal/mol. According to this definition, in our data set there would be 165 hot spots and 184 non hot spots. We have tested our approach in this case as well and report the results in Table 1b.

The list of the 20 protein structures and their interacting domains is reported in Additional file 2: Supplemental Tables S5 and S6. Note that homologous protein complexes are present in our data set but at least one of the proteins involved has less than 35% sequence similarity to its homologue. We remark that in these cases only a limited number of mutations occurs at equivalent sequence positions, as deduced from the alignment, and even removing them from the data set does not affect the overall results. As can be noted, the data set is dominated by the immunoglobin superfamily (10 structures out of 20 comprise at least one immunoglobin domain). We have verified a posteriori that this does not introduce a bias and that mutations on immunoglobins are not predicted with higher accuracy than on other proteins (see Results and Discussion section).

A detailed list of the individual mutations with their respective ΔΔG is reported in Additional file 3: Data set S1. We have analysed the distribution of amino acid types that occur in hot spots, in relation to the distribution of amino acid types in the database (Table 6). A similar analysis was performed in ref [19] and the results roughly agree with ours (e.g. tryptophan and tyrosine are among the most frequent amino acids in hot spots). We find however some notable differences, e.g. hot spots are not enriched in arginine (they are instead substantially enriched in lysine). Whereas the analysis in ref [19] was based on the whole ASEdb database, we have filtered out redundant entries and included only mutations located at a protein-protein interface for which the crystal structure of the complex is available. This might be the origin of the observed differences.

In recent publications [13, 26, 27], mutation data extracted from the Binding Interface Database (BID) [39] have been used to test proposed hot spot prediction models. In our opinion, these data have several problems (e.g. they are not associated to ΔΔG values) and it is questionable whether they are useful in assessing the predictive power of a method. Nonetheless, for completeness we have extracted from BID our own data set and tested our method on it. We report the results in Additional file 4 together with a discussion of our concerns about this data set.

Input features: energy components

As input features for the machine learning algorithms we have used basic energy terms that have been found to be important for the stability of protein complexes. These are van der Waals potential, hydrogen bonds, Coulomb electrostatics and desolvation energy. For each of the four energy components we have separately calculated 3 different energy contributions (schematised in Figure 1):

In total therefore there are 12 input features (4 × 3), although e.g. for Support Vector Machines we have used only a subset of them as this would give better or equivalent performances to the full set (see below for more details).

The energy components are calculated from the PDB structures. Only heavy atoms are considered and hydrogen atoms, if present, have been discarded. The computation of the different energy terms draws on established force-fields [11, 40, 41] but it incorporates some adjustments. It is important to remark that all terms are pairwise additive and therefore their calculation is formally equivalent. Before applying machine learning methods input features are rescaled so that they vary within similar numerical ranges. In the classification problem, each feature is standardised by subtracting its mean and dividing by its standard deviation. For regression purposes instead features are rescaled by their respective quadratic means. The latter normalisation makes it easier to search for solutions with no constant term, i.e. such that ΔΔG = 0 if all energy terms are zero.

van der Waals energy

The energy of van der Waals interactions between two atoms i and j is calculated using a "smoothed" 6 - 12 Lennard-Jones potential

(4)

where r is the distance between the two atoms and the parameters r _ij and ϵ _ij are taken from CHARMM19 force fields [40]. The parameter r₀ = 0.5 Å has the effect of widening the region of maximum affinity and to reduce the potential energy at r = 0 to a finite value. It makes the potential less sensitive to the precise position of atoms and to minor coordinate errors and local clashes. A similar option is available in the AutoDock suite [42]. We have also introduced a long distance cut-off, i.e. V _vdW (r) = 0 if r > 8 Å.

Support Vector Machine models

Support Vector Machines (SVMs) are playing an increasingly important role in the field of computational biology [20, 21]. They can be applied both to classification and regression problems. In the context of our investigation the former corresponds to predict whether a residue is a hot spot or not, the latter to estimate the actual value of ΔΔG. We have used the program package SVM^light [45], which is available at the website http://svmlight.joachims.org/. The program provides several standard kernel functions and we have experimented with linear, polynomial and Gaussian kernels, and with different combinations of input features. For classification purposes, we find the best results using a linear kernel and a set of 8 input features corresponding to the side-chain and environment inter-molecular energies. Using either a polynomial or a Gaussian kernel and a different (possibly larger) combination of input features, the performance would not improve, so we have opted for the simpler model. In the regression task, the best results are obtained with a linear kernel and by adding the van der Waals term from the side-chain intra-molecular energies to the features set (9 input features in total). The bias parameter in the kernel has been set to zero.

Standard SVM classifiers follow a supervised ("inductive") learning approach, whereby a predictive model is built on the basis of available positive and negative examples. The SVM^light package implements also a semi-supervised learning algorithm, the so-called Transductive SVMs (TSVMs), in which unlabelled data are exploited in order to develop a better discriminatory model. In our case, unlabelled data correspond to interface residues that have not been mutated to alanine and for which therefore ΔΔG is not known. The unlabelled data are extracted from the 20 protein complexes in the data set and used to integrate the original mutational data (maintaining the 16-fold partition described above for cross-validation). Optimal performance for TSVM classifier is obtained with a linear kernel and the same 8 features as for the SVM classifier.

Cross-Validation and Model Selection

SVMs have a number of tunable parameters (hyper-parameters) which should be chosen in order to achieve good generalisation performance and avoid over-fitting. For the classification task using a linear kernel, for example, two parameters, C and j, have to be set. The parameter C controls the trade off between the training error and the square norm of the weights associated to the features; the cost factor j determines how training errors on positive examples outweight errors on negative examples. An additional parameter is present in TSVMs, the fraction p of unlabelled examples to be classified into the positive class. In a regression setting the hyper-parameter ϵ controls the tolerance to errors. The optimal values for the hyper-parameters are not known in advance. In order to both evaluate fairly the performance of the methods and do not introduce biases in the choices of hyper-parameters, we have implemented a nested-loop cross-validation scheme [46] (see scheme in Additional file 2: Supplemental Figure S1) coupled with a grid search. The procedure ensures that hyper-parameters are selected without ever considering the performance on the test sets.

In nested-loop cross-validation, the data set is split into n folds (in our case n = 16 or n = 12). (n - 1) folds are used as the training set and the remaining fold is used as test set. The hyper-parameters are optimised on the training set, by applying a grid search and an internal (n - 1)-fold cross-validation strategy. The hyper-parameters that give the best cross-validated performance are selected and the associated model is then applied to the held-out test set. This process is repeated n times so that each of the original n folds is used as the test set once and predictions are obtained for each entry in data set. The procedure therefore consists of two nested cross-validation loops: an outer one for testing, an inner one for choosing hyper-parameters. In the inner cycle, we have assessed the model performance by means of the F1 score and the root-mean-square error (RMSE) for classification and regression tasks respectively (see below for more details).

For a given performance measure (e.g. the F1 score), we estimate its value f by considering the whole data set and by comparing predictions (obtained as described above) to the observed ΔΔG. The associated statistical error is instead evaluated as follows. Let {f _i }_{i = 1... n}be the set of optimised values obtained on the n training sets, as a result of the inner loops of cross-validation. We then write

(5)

Note that the more standard procedure of calculating the score and its error from the average and standard deviation on the n test sets is not applicable in our case as individual test sets are rather small and of different sizes.

We have verified that our approach does not produce over-fitting of the data by comparing results obtained from a 16-fold and from a more stringent 12-fold cross-validation. A significant better performance for the 16-fold cross-validation would suggest over-fitting. In addition we have tested the n models on their own training sets. Better results can clearly be expected in this case with respect to the test set but a large difference would possibly indicate over-fitting of the data.

Gaussian Process models

Gaussian Processes (GPs) are machine learning methods rooted in a Bayesian probabilistic framework [16]. They can be used for both regression and classification tasks. GPs assume that the data are a finite number of noisy observations generated by an unknown function. In line with Bayesian approaches [47], GPs assign a prior probability to each possible function which is then updated into a posterior probability in light of the observed data using Bayes' rule. Predictions on new data are made by a weighted average over all possible functions (models), the weights being equal to the posterior probabilities. Notice the different learning strategy between SVMs and GPs: SVMs determine one optimal model and make predictions based on it, GPs estimate a probability distribution over all possible models and predictions are model averages. One advantage of the latter approach is that it provides also confidence measures (i.e. error bars) on the predictions.

A GP is defined by a covariance function, which is the analogous of the kernel function in SVMs. Learning in GPs involves selecting a suitable covariance function for the problem at hand. As for SVMs, standard covariance functions include linear, polynomial and Gaussian functions. For both regression and classification we find the best results using a linear kernel and the full set of 12 input features. In the regression problem, the bias parameter of the kernel is set to zero. The data likelihood is assumed to be isotropic Gaussian. Notice that by using a linear kernel, a GP regression becomes equivalent to a Bayesian linear model with a Gaussian prior probability distribution over the linear coefficients. The GP classifier is computed based on the Laplace approximation [16] and using a logistic likelihood function. Intuitively, GP classification can be viewed as a generalisation of logistic regression (see e.g. [48]). For the calculations, we have used a software developed by one of us (CA).

Linear regression models

As a baseline regression method to compare results with we have built a linear model by least squares fitting (LLSF). The model has no constant term and the best fit is calculated by minimising the deviation between calculated and observed ΔΔG values. It includes all 12 energy features. We have also constructed "minimal" versions of the model based on only one or two features. Model parameters are calculated on the training sets and performances evaluated on the test sets, following similar procedia employed for SVM and GP regressions.

Measures of prediction performance

For classification, we primarily assess the prediction performances of our and related methods using the F1 score. Let TP, FP, FN refer to the number of true positives, false positives and false negative respectively. Precision (P, also called specificity) and recall (R, also called sensitivity) are defined as

(6)

The F1 score is the harmonic mean of precision and recall

(7)

Precision, recall and the F1 score are comprised between 0 and 1, the larger the value the better the performance. A random predictor on our data set would score P _ran = 0.23, R _ran = q and , where 0.23 is the fraction of hot spots in our database and 0 ≤ q ≤ 1 is the fraction of presumed hot spots in the predicted set. It follows that 0 ≤ F 1 _ran ≤ 0.37 and if the expected frequency is equal to the observed frequency of hot spots (i.e. R _ran = q = 0.23), then F 1 _ran = 0.23. For each prediction set, we have also calculated the Matthew's correlation coefficient (MCC) given by

(8)

where TN is the number of true negative and TP, FP and FN are as above. The MCC value is between -1 and 1: a perfect predictor has MCC = 1 whereas for a random predictor MCC = 0 (MCC = -1 for a perfect inverted predictor)

The performance of the regression models is evaluated using the root-mean-square error (RMSE) and the correlation coefficient (r) between the predicted (calc) and the observed (obs) ΔΔG

(9)

(10)

where sums are over the N entries in the data set and the symbols ⟨·⟩ and σ denote averages and standard deviations respectively. A regression model can be easily mapped into a classifier by associating a residue to an observed or predicted hot spot if ΔΔG ≥ 2 kcal/mol. The same performance measures used in the classification problem can then be used to assess the regression models as well.