Recognizing protein–protein interfaces with empirical potentials and reduced amino acid alphabets

Fold recognition for monomeric structures

The energy function is a sum over pairs of residues less than 4.5 Å apart. Each pair of residue types has its own interaction parameter. Several parameter optimizations were conducted, using the amino acid alphabets in Fig. 1. which contain 2, 3, 4, 6, and 20 amino acid types, respectively. We considered 810 native protein structures, along with sets of decoys constructed by threading each sequence onto unrelated structures. The structures and their decoys were split into an Optimization set (615 structures) and a Test set (200 structures). The energy parameters were chosen to maximize the $\overline{Q}$ parameter (Eq. 7), summed over the proteins of the Optimization set. The quality of each energy function can be measured by the final, optimized $\overline{Q}$ value, which ranges from 0 to 100% and reflects the ability of the energy function to group near-native structures in the low-energy range. A second quality measure is the "discrimination percentage" D: the percentage of native structures in the Optimization set that are successfully predicted to be of lower energy than all their decoys. The quality measures are obtained by cross validation; i.e., we compute $\overline{Q}$ and D for the structures of the Test set, which are not used for parameter optimization. The quality of the predictions are summarized in Table 1.

Number of amino acid groups	Optimization Set		Test Set
	$\overline{Q}$(%)	D (%)	$\overline{Q}$(%)	D(%)
2	74.3	69.8	74.4	70.0
3	76.4	72.2	75.3	72.0
4	87.5	84.4	87.4	84.6
6	88.3	86.7	91.9	89.4
20	94.4	91.9	91.3	88.4

The most detailed energy function, with all twenty amino acid types, is subject to overoptimization, or overtraining: performance on the Test set is 3–4% below that on the Optimization set (Table 1). In fact, the performance on the Test set is slightly inferior (by 1%) to that of the best six-class energy function. For the simpler alphabets (2–6 types), the performance on the Test set was as good, or slightly better than on the Optimization set.

Fig. 2 shows the discrimination D for different protein size ranges, for all alphabets. For proteins of more than 320 amino acids, discrimination is 100% with either four, six, or twenty residue classes. With two and three amino acid classes, discrimination reaches 100% when the protein size reaches 250 amino acids, except for a few large proteins. The lower discrimination for small proteins is partly due to their higher number of decoys (also shown).

The amino acid groupings used here (Fig. 1) were derived from the Blosum50 similarity matrix by Levy et al [30]. Residues with similar physicochemical properties are clustered together. The parameters show a good consistency with the hierarchy of amino acid types, as discussed above. When an amino acid group is subdivided, the parameters corresponding to the subgroups are usually in qualitative agreement with the parent group. This is not due to bias in the optimization procedure, since the best parameters at each level were usually obtained from a random scan of parameter space (see Methods). Another six-class grouping ({AVLIMC}, {FWYH}, {STNQ}, {KR}, {DE}, {GP}), proposed by Shakhnovich and coworkers [31], gave significantly poorer results, with a D of 71%.

Recognizing protein–protein interfaces

Characterizing the structures

To parameterize and test the energy functions, we used 219 protein–protein complex structures, including 195 homodimers and 24 heterodimers. We constructed about 1300 decoys each (for a total of about 290,000 decoys). The decoys were constructed by a realistic (and computationally expensive) docking procedure, with flexible proteins interacting through a molecular mechanics energy function. Some of the native structures' properties are summarized in Fig. 3. The distribution of interface sizes has a peak at about 1800 Ų, with a broad tail extending up to about 6500 Ų. The propensities P _X of each amino acid type X to be found at the interface are also shown. P _X is defined as

$P_X=\log f_X^{\mathrm{int}}/f_X^{\text{SWP}},$

(1)

where $f_X^{\mathrm{int}}$ and $f_X^{\text{SWP}}$ are the proportion of amino acid type X found in the protein–protein interfaces and in the SwissProt data base, respectively [13]. A positive propensity P _X means that the proportion of X in the interfaces is higher than the overall proportion of X in SwissProt. Hydrophobic types, especially Met, Trp, Tyr and Phe are overrepresented, compared to their natural abundance in proteins. Among the polar amino acids, Arg and His are overrepresented, and Lys is underrepresented. Homodimeric interfaces (the majority of our data set) tend to be larger and more hydrophobic than heterodimers. Nevertheless, the trends in Fig. 3 agree with recent surveys of protein–protein interfaces [12, 13].

To illustrate our decoy sets, Fig. 4 shows two native structures (PDB codes 1ARO and 1BJF) along with a simplified representation of their decoy series. The series contain 1573 and 1652 decoys each. Decoys are produced by fixing one or the other of the partners, and positioning the other using a flexible docking procedure (see Methods). In Fig. 4. one partner ("A") is arbitrarily taken as a reference, and the distribution of the other ("B") around it is shown (one structure per decoy). Partner B is seen to be widely distributed around A, with a wide variety of orientations. By construction, each decoy has a native-like interface size and a limited distortion of the internal structure of both A and B. The energy distributions for these same series of structures are discussed in the next section. The decoy series for a third, typical, dimer is illustrated in Fig. 5 in a different way. For each amino acid of one member of the dimer, Fig. 5 shows the number of decoys where it participates in the protein–protein interface. We see, again, that the surface coverage of one partner by the other in the decoy sets is rather thorough.

Recognition power with different amino acid alphabets

The energy functions for interface recognition were optimized using a combination of monomeric and dimeric structures (see Methods). We considered two Optimization sets, containing respectively 110 and 109 dimeric structures and their associated decoys, along with 110 monomeric structures and their decoys. We denote them OS1 and OS2. The dimeric structures not belonging to OS1 form the Test Set 1; similarly for OS2. The amino acid subdivisions into classes are the same as for the monomeric energy functions above. The prediction results are summarized in Table 4. The optimized energy parameters are in Tables 2, 3 and in Supplementary Material (for the 20 class function: see additional file 2: Supp2.ps). To measure performance, we employ D _int , the percentage of native interfaces ranked above all their decoys, and ${\overline{Q}}_{int}$, the average of $\overline{Q}$ over all the protein–protein complexes. D _int and $\overline{Q}$ are exactly analogous to D and $\overline{Q}$ (Eq. 7), but they only include the interface recognition statistics; fold recognition statistics are not included (whereas they are taken into account during the optimization).

Number of amino acid groups	Optimization Set 1 a				Optimization Set 2 b
	Opt. c	Test. d			Opt. c	Test d
	$\overline{Q}$	$\overline{Q}$	${\overline{Q}}_{int}$	$D_{int}^e$	$\overline{Q}$	$\overline{Q}$	${\overline{Q}}_{int}$	$D_{int}^e$
2	89.4	90.7	87.5	46.7	90.1	90.1	87.2	47.3
3	93.0	92.7	89.0	50.5	93.0	92.7	89.9	56.4
4	97.9	96.7	93.5	60.7	97.9	97.2	94.6	68.2
6	98.0	97.0	93.7	63.5	98.0	97.7	95.4	69.1
20	99.1	97.6	93.9	67.3	98.7	98.2	94.8	69.1
20/Bastolla e	-	-	96.5	56.4	-	-	93.3	62.6
20/Lu e	-	-	93.6	63.6	-	-	93.5	63.5

^a Results for the parameter sets produced with OS1. ^b Results for the parameter sets produced with OS2. ^c Results on the Optimization set. ^d Results on the Test set. ^e ${\overline{Q}}_{int}$ and D _int include only the performance for interface recognition (the fold recognition statistics are left out). ^e With the Bastolla or Lu energy functions.

The optimal energy parameters with the binary alphabet are consistent with the ones derived for fold recognition (Table 2). The parameters obtained using OS1 and OS2 are almost identical. Discrimination is fair, with a D _int of 47% and a ${\overline{Q}}_{int}$ of 87% for both the OS1 and OS2 sets. The energy parameters with the ternary alphabet are consistent with the binary ones. Contacts between members of the two "hydrophobic" classes are stabilizing, for example. The discrimination power has increased: ${\overline{Q}}_{int}$ = 89% and D _int = 53% (averaging over the OS1 and OS2 results).

With the four-class alphabet, ${\overline{Q}}_{int}$ and D _int increase sharply, to 94% and 64%, respectively. Consistent with the monomeric results, contacts within and between the {LIVIMC} and {FYW} classes are the most stabilizing. Optimal parameters from OS1 and OS2 are similar. The rms deviation between the OS1 and OS2 energy matrices (with four amino acid classes) is 3.6 kcal/mol.

Performance improves only slightly with six classes: ${\overline{Q}}_{int}$ = 95% and D _int = 66%. With the complete, 20-class alphabet, there is almost no further improvement: ${\overline{Q}}_{int}$ is 94% and D _int is 68%. It appears that for the simple, residue–residue model used here, 4–6 amino acid classes are sufficient to encode the actual complexity of the model.

To illustrate further the behavior of our energy functions, we consider the structures 1AR0, 1BJF, and 1FBT, which represent three typical situations. The energy spectra for the three decoy series are shown in Fig. 6, using the energy functions with two, six, and twenty amino acid classes.

• 1AR0: With 20 classes, discrimination is successful. The native complex is very stable compared to the decoys.

• 1BJF: The native complex is not recognized as the lowest-energy structure, but it has the third-best energy and is well-separated from the rest of the energy distribution (the decoys lying higher than it).

• 1FBT: This structure illustrates the worst case. The energy function is unable to recognize the native complex at any level of alphabet complexity. The native energy is quite high. The position of the native structure in the energy spectrum is only 0.5 standard deviations below the average decoy.

To illustrate the diversity of the decoy sets, we also calculated energy spectra with a "random energy" model. Decoys are not explicitly built; rather, we count the number of interface contacts in the native dimer A:B, say N contacts; we replace them by N contacts established between N random pairs, where the first amino acid is drawn from the surface of A and the second from the surface of B. Given the rather large numbers of interface constants, the corresponding energy spectra are essentially gaussian, with a standard deviation that is easy to compute from the composition in amino acid types (not shown). The distribution mean is taken to be the same as that of our decoy sets. The 2-class random energy spectra are shown in Fig. 6 (lefthand panels; thick lines). We see that our decoy sets lead to more diverse energies than the random energy model.

Finally, for 1AR0, 1BJF, and 1FBT, Fig. 7 illustrates our ability to discriminate decoys as a function of their deviations from the native structure. For the best case, 1ARO, there is a clear energy increase as one moves away from the native structure (large dot, lower left). For 1BJF, the native structure has a low energy, while the closest decoy structures are much higher. For 1FBT, the energy/structure correlation is fair, even though the native structure is not well-discriminated; the decoys that have a lower energy correspond to very different structures with a different, competing binding mode. Energy/structure correlations will be studied in more detail in future work. Notice that the decoy sets typically contain some, but not very many structures that are very close to native. This is not the goal of our docking approach, which focuses on the structural quality of the decoys (reasonable van der Waals and electrostatic interactions), the size of the interfaces (which should not be too small), and seeks to produce diverse interfaces. For the present class of energy functions, the number and nature of interface contacts are the most important decoy properties, and we believe these are adequately sampled with our procedure; this should lead to a successful parameterization. In contrast, we systematically include near-native structures in the blind testing of the energy functions further on.

Analysis of the energy parameters and amino acid classes

The simplified amino acid alphabets used here were derived from a clustering procedure based on the Blosum50 similarity matrix ("Blosum clustering"; see Methods, Eq. 4). An alternative classification can be obtained by performing a clustering based on the optimized energy matrix, constructed with all 20 amino acid types ("energy clustering"; Eq. 5). By comparing these two classifications, we can assess the robustness and consistency of the simplified alphabets. The two methods give somewhat different clusters, as shown in Fig. 1. At the lowest level, the energy function does not give a meaningful separation into two classes. It gives the same hydrophobic cluster as the Blosum50 matrix, with a good internal coherence (Pearson Correlation Coefficient of R = 0.58; Eq. 5). But the Blosum50 polar group is split in two by the energy clustering. At the six-class level, the two methods give the exact same amino acid groupings.

Fig. 8 illustrates the best energy function for each level of complexity, using a contour representation. The functions are seen, again, to form a coherent hierarchy. We note that with three, four, and six amino acid types, the best function was obtained independently of the simpler functions. Only the 20-class function was obtained using one of the simpler functions as a starting point (the 6-class/OS2 set, shown on the lower right). Thus, the consistency between the energy functions with the largest and the smaller alphabets is largely independent of the optimization procedure.

To further check the internal consistency of our hierarchy of energy functions, we take a "top-down" approach, using the complete, 20-class function as a starting point to infer simplified functions, which can be compared to the simplified functions described above. To infer a 6-class function, for example, we simply average the appropriate elements in the complete matrix. Thus, the lines and columns corresponding to K, R, and H in the complete 20 × 20 matrix form a submatrix, which is averaged to give the {KRH}{KRH} interaction in the inferred 6-class matrix. We also compute the standard deviation over the submatrix, which measures the internal consistency of the {KRH} subgroup within the complete matrix. Results for 4 and 6 classes are given in Table 5; results for 2 and 3 classes are qualitatively similar and are not shown. The simplified functions inferred in this way are in general agreement with the functions described above, obtained by direct optimization of $\overline{Q}$.

	H α	H β'	H γ	P α	P β	G
	-2.03(0.07)	0.13(0.06)	-1.70(0.05)	0.79(0.04)	0.54(0.04)	-0.13(0.03)	H α
		0.01(0.05)	-0.43(0.05)	0.16(0.05)	-0.22(0.05)	0.11(0.04)	H β'
H α	-2.03(0.07)		-0.65(0.02)	0.37(0.05)	-0.22(0.06)	-0.56(0.06)	H γ
H β	0.08(0.11)	0.05(0.07)		0.55(0.05)	-0.27(0.07)	0.03(0.07)	P α
H γ	-1.70(0.05)	-0.45(0.07)	-0.65(0.02)		0.56(0.08)	0.24(0.13)	P β
P	0.71(0.10)	0.02(0.19)	0.12(0.30)	0.15(0.42)		0.1(0.)	G
	H α	H β	H γ	P

Evaluating the energy functions in blind tests

The quality of interface recognition

An important goal of modern genomics is to identify interacting proteins and characterize the structure and function of the corresponding complexes. A promising strategy is to assemble putative complexes from individual protein structures and characterize them by computing a binding free energy. Here, we have developed simple, empirical energy functions for interface recognition. We used a coarse-grained description, at the residue level, along with the simplest possible, stepwise, distance-dependence. More compex models are known to give superior results [10, 18, 19, 21, 26, 33]. Here, however, we pursued two more limited goals. First, we examined whether a large computational effort would lead to better parameterization and performance of this simplest class of models. Second, we examined what size amino acid library is needed for optimal performance.

Thanks to increased computing power, we employed many more structures and many more decoys than previous workers. The decoys were both diverse and structurally-realistic (for example, in terms of interface size), thanks to a realistic, flexible docking method of construction. In particular, van der Waals and electrostatic contacts between partners are reasonable. The decoy set was characterized in some detail; more details will be given elsewhere. Our search of parameter space was also considerably more thorough than some past work. This is important for the 6-, and especially the 20-class amino acid alphabets, which have 21 and 210 adjustable parameters, respectively. For these, we considered a total of over 140,000 starting guesses for the parameters. We also included fold recognition in the optimization criterion (using additional sets of target structures and decoys), which should help to restrain the parameters to physically meaningful values. Finally, by constructing a hierarchy of energy functions, with amino acid alphabets of increasing complexity, we could focus the parameter search effectively for the larger alphabets. For example, our best 20-class function was obtained using a 6-class function as a starting point. The performance of our method is measured by its ability to identify the experimental, native interface from among sets of decoys.

Despite the improvements made here, the performance of our functions is very similar to those (of the same complexity) developed by Bastolla, Lu, and coworkers. For the CAPRI tests, we had six successes out of 13 targets with our 6- and 20-class functions, while our 4-class function actually gave a seventh near-success (the native T04 structure, ranked 11th). The Bastolla function gave seven successes and the Lu function, six. The successful cases are almost the same with all the functions, and behavior for near-native decoys is similar. This is a bit disappointing, given the extensive computational effort in the present work. However, this effort was needed to establish firmly that the performance limit of this class of models has been reached.

A coherent hierarchy of amino acid alphabets and models

Our second goal was to examine a hierarchy of amino acid alphabets and models. This analysis has produced new insights into both the similarities and differences among amino acids, and the diversity required for protein–protein recognition. The amino acid groupings were initially obtained by clustering the Blosum50 similarity matrix [30], which is itself based on the mutational frequencies seen in multiple sequence alignments [34]. The mutational frequencies are mainly determined by factors other than interface recognition. Presumably, the dominant mutational constraint is the fold and stability of the individual protein [35]. Protein–protein recognition is one of several secondary constraints, along with folding kinetics and robustness towards mutations [36]. Nevertheless, the Blosum clustering leads to good quality interface recognition. At the 4-class level, for example, two thirds of the native structures in our data set are ranked first, among over 1300 high-quality decoys. This is consistent with the idea that similar physical-chemical effects operate in protein folding and protein–protein association.

The best energy functions at each level were obtained from random initial guesses, except for the 20-class function, which was obtained from a 6-class function. Nevertheless, the energy functions form a consistent hierarchy, with qualitatively similar parameters at all levels. This is especially striking in the contour representation of Fig. 8. From the complete energy matrix, we can also infer amino acid groupings (left of Fig. 1). The 20-class matrix gives back the Blosum 6-class grouping, as expected (because of its construction from a 6-class function). At the 4-class level, there is an interesting difference between the two groupings. Blosum splits the aromatics {FWY} into their own group, whereas the energy matrix classification splits the polar amino acids into two groups, based partly on electrostatic charge: {DENQ} and {HKR}. At the 6-class level, polar–polar interactions are repulsive, except for the {DENQ}{HKR} interaction, which is reasonable. On the other hand, F, W, and Y have much higher interface propensities than the other hydrophobics except Met, so that an aromatic group is also reasonable. The 6-class function captures both of these effects. The internal coherence of the 4- and 6-class energy groupings is very high (see correlation coefficients of each group, left of Fig. 1). In future work, it would be interesting to also consider 5-class functions.

Monomeric structure selection

We used 810 X-ray structures from the PDBselect collection [37]. Sequence identity between all pairs was less than 25%. The structures were monomeric and did not contain any cofactors or metal ions. Crystallographic resolution was better than 3.3 Å (average of 2.0 Å). They were divided into two sets, having similar average chain lengths (around 190 amino acids). 615 proteins were used as an Optimization set, for the parameterization of the energy functions below. The other 200 were used as a Test set, for evaluating the parameters' performance.

Monomeric decoys

For each native structure, we performed a gapless threading of the sequence onto segments of matching size from all the structures in our 810 protein set. On average, we obtained 349 decoys per native structure. For each amino acid position in a decoy, the most probable rotamer for the given amino acid type was arbitrarily taken from the Dunbrack backbone-independant rotamer library [38] and constructed. Thus, steric constraints could be violated in the decoy structures, because the rotamers were not optimized with respect to the particular backbone and sidechain environment. With the energy functions used below, this is not a difficulty. Furthermore, this simple method for constructing monomeric decoys was deemed sufficient, because our main focus is not fold recognition, but the identification of protein–protein complexes.

Dimeric structure selection

We used 219 dimer structures from the DIMER-1 set of Lu et al [19], which contains 340 structures. We retained only dimers, as opposed to higher-order complexes. This distinction is not always straightforward, and one higher-order complex was identified in the data set at a late stage. It is not expected to affect the results, judging by the similar parameters obtained with OS1 and OS2, or when an optimized parameter set is refined further with the higher-order structure left out. All the structures have a crystallographic resolution better than 2.5 Å, 30 interface contacts or more (using a 4.5 Å distance cutoff), and a pairwise sequence identity with any other member of the set below 35%. All correspond to biological dimers, as opposed to crystal contacts. This was established by running an automated test with the PQS server [39] and by checking the literature.

Our final set included 195 homodimers and just 17 heterodimers. A second set of 23 heterodimers was set aside to perform blind tests after parameterization was finished. The predominance of homodimers was largely imposed by the structures available in the PDB [40]. A typical homodimer is expected to have a larger, somewhat more hydrophobic interface, and a stronger binding constant, compared to a typical heterodimer. This is both a limitation and an advantage for our purposes. We expect that there will be a stronger dimerization signal in the homodimer interfaces and, furthermore, the signal is repeated twice, since the complexes are almost always symmetrical. This should facilitate parameterization, without damaging the heterodimer performance. Indeed, heterodimer association is driven by largely the same physical effects as homodimerization, so that the easier homodimer parameterization should carry over and be applicable to heterodimers. This is confirmed by our blind tests on the separate set of 23 heterodimers, which are not used in any step of the parameterization, and are not homologous to any of the proteins used.

Dimeric decoys

where N_dec,N _nat are the number of residue–residue contacts in the decoy and native structure, respectively. If F was at least 80%, the decoy coordinates were output at this stage. If not (the vast majority of cases), minimization was continued, without the pulling restraint, but allowing now intramolecular deformation of one of the partners. Also, an electrostatic contrast was introduced by enhancing the atomic charges on the two partners (with respect to the Charmm19 values), to increase the interprotein attraction. The charges of one partner were all increased by 0.25e, while the charges of the other partner were decreased by 0.25e. Intraprotein sidechain–sidechain electrostatic interactions were omitted and intraprotein sidechain-backbone electrostatic interactions were reduced by 1/2.

After 50 steps of Powell minimization [42], we compared the decoy structure to the native one. For the "deformable" partner, we performed a best fit of the monomer on its native structure and evaluated the "intramolecular" rms deviation between the two, considering its sidechains only. If the sidechain deviation was above 4.5 Å, the decoy was discarded. Decoys with F < 0.45 were also discarded. Decoys with a van der Waals energy greater than 4000 kcal/mol were also discarded. If the decoy was too similar to the native complex it was discarded. The criterion for similarity was the rms deviation for the "mobile" partner in the decoy and the native complex: if the rms deviation was less than 3 Å, the decoy was discarded. The deviation was computed without superimposing the two structures, so that it measures the overall motion of the mobile partner away from its starting, native position. Finally, to eliminate redundant decoy structures, we computed rms deviations between all pairs of decoy structures, rejecting decoys that were within 3.5 Å of another decoy. In all, we obtained a total of about 290000 decoys, for 219 native complexes, giving an average of 1326 decoys per native complex, with a minimum of 277 and a maximum of 1843 decoys. Typically, several tries were needed per decoy, so that the first step of the protocol (random displacement of one partner) was done over one million times.

To assess the quality of decoy structures, we used a detailed, atomic energy function, corresponding to the Charmm19 molecular mechanics force field [41]. The solvent environment was modelled with the ACE variant [43] of the Generalized Born model, parameterized as in [44]. These calculations were done with the XPLOR program [45].

Empirical, residue-based energy function

where the sum is over all pairs of amino acids i, j in the protein. U(S _i , S _j ) is an interaction energy, which depends on the amino acid types S _i , S _j at positions i, j. C _ij is one if the pair i, j is close together, and zero otherwise. Specifically, i, j are assumed to interact if they have at least one pair of nonhydrogen atoms within a distance of 4.5 Å, and if they are not first or second neighbors along the polypeptide chain; i.e., |i – j| ≥ 3. We refer to U as an energy matrix, and to C as a "contact map".

With a complete amino acid alphabet of twenty types, the 210 elements in the energy matrix U are all different. With a reduced alphabet, amino acid types belonging to the same group share the same coefficients in the energy matrix U. For example, with a binary, hydrophobic/hydrophilic alphabet, the 210 elements in U fall into just three categories, and have just three possible values, corresponding to hydrophobic–hydrophobic, hydrophobic–hydrophilic, or hydrophilic–hydrophilic pairs. In what follows, we consider alphabets with two, three, four, six, and twenty classes.

Reduced amino acid alphabets

M is the Blosum50 matrix [34], and the sums are over the twenty amino acid types. The clustering method assigns the two amino acid types with the highest similarity to a group. Then, the pair with the next-highest similarity value is considered. If one member of this pair belongs to the first group, the other member is added to that group. If not, the pair is assigned to a new group. This process is repeated until all the amino acid types are divided into a predefined, desired number of groups.

Here, x and y are two different amino acids, U_{x, i}is the pairwise energy of interaction between the residues x and i, $U_{\overline{x}}$ is the average over the U_{x, i}, and similarly for y.

Optimizing the energy function: general method

where the sum is over all the decoys, C is the contact map of a decoy, C _nat is the native contact map, k is Boltzmann's constant, and T is the absolute temperature. If the Boltzmann-averaged overlap Q is close to one, then the lowest-energy conformation is either equal to the native structure or very close to it. At physiological temperatures, a high value of Q implies that structures very different from the native one have high energies, while native-like structures have low energies. The energy landscape is said to be well-correlated [20].

where t, t + 1 are successive iterations; the sum is over the proteins S in the Optimization set; q₀ is the overlap between the minimum energy conformation and the native structure of S; δ and γ are numerical parameters, chosen to facilitate convergence to a maximum. Typical values were 0.2 and 0.0075. For more details, see [20].

For monomer structure recognition, we optimized the energy function using the 615 native structures in our Optimization set, along with their associated decoys. For dimer recognition, we optimized the energy function using a combination of monomeric and dimeric structures. We considered two Optimization sets, each containing 109 or 110 dimeric structures and their associated decoys, along with 110 monomeric structures and their decoys. We denote the Optimization sets by OS1 and OS2. The dimeric structures were assigned randomly to OS1 or OS2. The monomeric structures were chosen so that their average chain length is equal to that of the Monomeric Optimization Set, above (189 amino acids). The OS2 dimer structures form the Test Set 1 (TSl); the OS1 dimer structures form Test Set 2 (TS2). The same monomeric structures were used for both OS1 and OS2. Notice that cross-validation, below, is done using the performance of the energy functions on interface recognition only. Parameter bias towards the monomer structures can exist but is not a concern. Notice also that with the normalization of Q in Eq. (6), the monomers and dimers have the same weight in the combined optimization, despite the larger number of decoys for the dimers.

Optimizing the energy function: starting parameters

For the more complex alphabets, the number of energy parameters to be optimized is large (210 in the case of the complete amino acid alphabet). This means that there are many local maxima for $\overline{Q}$, and the optimization must be done repeatedly, with multiple starting points. We used three methods, both for the monomeric and dimeric protein sets. For the simplest alphabets, with two or three amino acid types, we performed a systematic scan of possible parameter choices. The range of values considered was from -10 to + 10 kcal/mol, with a step of 0.25 (binary alphabet) or 0.5 (ternary alphabet). For the larger alphabets, a systematic scan was no longer possible; instead, we used 64000 random starting points, with each parameter drawn randomly from the same range of -10 to +10 kcal/mol. In addition, for alphabets with three amino acid types or more, we used starting points taken from the optimization of the smaller alphabets. Indeed, the alphabets form a hierarchy, where the finer groupings are subdivisions of the larger groupings. For example, the binary alphabet uses the groups {LVIMCAGSTPFYW} (hydrophobic) and {EDNQKRH} (hydrophilic), while the ternary alphabet splits the hydrophobic group into two subgroups.

The ternary optimization can be started from the binary optimum, assigning initially the same, "binary" parameters to the two hydrophobia subgroups. As the optimization proceeds, the parameters associated with two subgroups will shift to different values.

For the alphabets with 6 and 20 classes, to explore the effect of an overall scaling of the energy parameters, we also generated sets of random starting parameters with a fine-grained sampling of possible values, limited to the range -1 to +1 kcal/mol. These parameter sets were optimized in combination with four different values of the temperature T (Eq. 6), corresponding to kT products of 0.6, 1.0, 1.5, and 4.0 kcal/mol. 2000 parameter sets were optimized with each temperature, for each each alphabet, giving 8000 additional parameter sets per alphabet.

Testing the energy function

To test the potentials, we used a cross validation procedure. As described above, the set of native structures were split into Optimization and Test sets. Only the Optimization set was used in Eqs. (7, 8), above. The total $\overline{Q}$ was then evaluated for the Test set. A value close to 100% means that the native structure is almost always ranked as the lowest-energy conformation. We also compute a "discrimination" percentage. We mostly employ a "strong" discrimination, D, which is the fraction of native structures that are ranked as the lowest-energy conformation (compared to their decoy set). We occasionally compute a "weaker" discrimination, D _k , which is the fraction of native structures that are ranked among the k lowest-energy conformations (so that D = D₁).

Another form of cross validation was done by performing several blind tests. In one test, we used our energy functions to rank native and decoy structures available on the web, from Vakser, Sternberg and coworkers. Vakser's decoys sets were produced by a rigid body docking method, using the GRAMM program at "high resolution" (grid step ≤ 2 Å) [27]. They correspond to five biological dimeric complexes. Each decoy series includes 99 decoys. A few of these are within 5 Å (rms deviation) of the native structure; the others are farther away. The Sternberg decoys were produced by a rigid body docking method using the MultiDock program [28]. They correspond to ten biological dimeric complexes, including one that is present in the Vakser set (with different decoys). Each series includes 99 decoys. An electrostatic filter is used during the docking simulations, and all ten decoy series include at least three decoys that are close to the native structure. The complexes are listed in Supplementary Material.

Finally, in another blind test, we applied our energy functions to native and decoy structures of protein dimers submitted to the CAPRI experiment for protein–protein complex structure prediction (rounds 2–5) [8, 9]. We considered 13 native structures, with an average of 173 decoys each and a minimum of 66.