Integrating cross-linking experiments with ab initio protein-protein docking

Overall approach

We used ZDOCK [18-21] with input component proteins obtained from X-ray crystallography or through homology modeling using X-ray crystallography template structures. The ZDOCK algorithm was integrated with experimental cross-linking data to generate only predictions that satisfy the cross-links. The following three approaches were tested: (1) Filtering the predictions from a standard ZDOCK calculation using the Euclidean distance between cross-linked sites. Although Euclidean distances are fast to compute and therefore applicable to large sets of predictions, the cross-linking distances could be underestimated because the Euclidean path is allowed to pass through protein-occupied space. (2) Filtering the ZDOCK predictions using the Xwalk algorithm [27,36] to determine the shortest path that is allowed to only pass through protein-unoccupied space (void-volume). Although physically more accurate than Euclidean distances, computationally the grid-based algorithm is orders of magnitude more expensive to evaluate. (3) Restrict ZDOCK to search only the space that satisfies the Euclidean cross-linking constraints. This approach yields more retained predictions than the filtering methods and therefore may improve performance.

We performed cross-linking and mass spectrometry experiments with the lysine-reactive BDP-NHP chemical [30] and then used the ReACT [30] algorithm to identify the cross-linked sites. We used our previously published cross-linking data [25,29-34,37] and unpublished data, and because this data was all generated with the same cross-linking chemical, it allows for a systematic study of the computational algorithms. We then retained only the 2000 heteromeric cross-linked proteins, as most ab initio protein-docking algorithms are designed to predict such complexes. We plan to investigate homomeric complexes in future studies. A search of the Protein Data Bank (PDB) [35] by requiring a minimal sequence identity of 30% resulted in 219 complex structures. However, many of these were not suitable for our first effort of developing a crosslink-guided protein-protein docking algorithm. We excluded 44 ribosome complexes because of interactions with nucleic acids. Seventy-four complexes were excluded because the interaction required three or more components. The sequence identity of the components was so high for 28 pairs that they were effectively homodimers. For 12 pairs we only found bound structures but no unbound structures, and we further removed the complexes that cannot be handled by docking algorithms because they showed peptide-like interactions (9 structures), co-folded chains (6 structures), or had one protein wrapped around the other (9 structures). Twenty complexes were excluded for various other reasons (covalently bound proteins, no direct interaction between proteins, incomplete structure, only C_α coordinates listed, and low-sequence-identity complexes that do not agree with cross-linking data—each of these reasons for exclusion was seen at most five times in our set of 219 structures). The resulting collection of 19 complexes was used to assess the three approaches for integrating cross-linking experiments and protein-protein docking. We realize proteome-wide cross-linked peptide-guided docking will require inclusion of many additional capabilities, such as co-folding, homo-multimers, or other considerations. However, the efforts presented here with 19 complexes demonstrate how the combination of large-scale interactome data with computational advances can correctly predict complex structures and give motivation for developments on this front.

Test set

Ideally, a protein-protein docking test set has the unbound structures available for all component proteins, or unbound templates that can be used for homology modeling the components. However, in order to maximize the number of entries in our set, we allowed ten tests that had one of the components in the bound form. One of these cases had the unbound structure for the other component available (unbound/bound docking), and for the remaining nine the other component needed to be homology modeled (homology/bound docking). For two cases in our set we had unbound structures for both components (unbound/unbound docking), for five cases only unbound templates (homology/homology docking), and for the two final cases one template and one unbound structure (homology/unbound docking). The 19 complexes could be divided into 11 groups based on fold-similarity and are summarized in Table 1 and shown in Figure 1.

The I-RMSD values in Table 1 indicate the conformational differences between the bound and unbound interfaces [6]. Except for 2C, 2D, 3A-C, 6, and 10 that had I-RMSDs ranging from 2.5 to 5 Å, the other complexes had I-RMSDs under 2.5 Å and would be classified as having low to medium difficulty for ab initio docking algorithms such as ZDOCK [6]. Table S1 lists the experimentally detected cross-links. Case 3C had 35 cross-links, and other cases between 1 and 7 with an average of just over 2. Based on the bound structures, the Euclidean distances between the cross-linked sites were under 35 Å and the void-volume distances under 40 Å, except for cases 2C, 2D, and 9, whose cross-linked sites were at slightly greater distances, and case 3C, which showed several much larger distances, which are likely related to higher order oligomers that are not reflected in the crystal structure that we used [38]. The distance distribution is consistent with previously reported data [29].

Three of the groups (1-3) display multiple interfaces in the bound structures. When they involved different binding sites they were considered both separately and combined in the assessment of the prediction algorithms. The first example, group 1, consists of the complex of tubulin α and β chains (1). The bound structure contains a stathmin chain as well, although it was not incorporated in the docking experiments [39]. The structure shows three interfaces between the components, two of which are distinct and involve two different binding sites for each component.

Group 2 consists of six members, human (2A) and non-human (2B-F) Succinyl-CoA ligases, each being a complex of a α-subunit and a β-subunit. The β-subunit consists of two domains that both bind the α-subunit. Although some sequence identities amongst the members of the group were high (2A and 2E over 90% for both components), we considered these complexes as distinct test cases because there is little redundancy among the cross-links. Only the single cross-link of 2C was equivalent to one of the cross-links of 2B, see Table S1. Furthermore, the stoichiometries were different for 2A and 2B, for which crystal structures with nearly exact sequence identities were available. One α-subunit and one β-subunit form the heterodimer 2A [40], whereas 2B is a tetramer with a homodimer of two β-subunits at the center [41]. Indeed, one pair of the cross-linked sites for 2B (P0AGE9 residue 43 with P0A836 residue 172) requires the tetramer structure, whereas both pairs of cross-linked sites of 2A are consistent with a dimeric complex structure. Consequently, we used the monomeric β-subunit as the docking input for 2A and the dimeric β-subunit for 2B. Based on the high sequence identities among 2B, 2C, and 2D as well as the need of the dimeric β-subunit to rationalize some of the cross-links, we used the dimeric β-subunit structure as the input for docking 2C and 2D. Complexes 2E and 2F have high sequence identity with 2A and the observed cross-links did not require the tetramer structure, thus we used the monomeric β-subunit as input structure.

Group 3 consists of F1-ATP synthases. The bound structures show three α-subunits and three β-subunits that are close to the C₃ symmetry but broken by the γ-subunit that binds to the center of the complex. The α and β-subunits have the same fold, but within each complex, the sequence identities are only 26%. We assumed that the synthases were stable without the γ-subunit and ignored the γ-subunit in docking and subsequent analysis. This is supported by the thermophilic Bacillus 1-ATP synthase, which has been crystallized both with and without a symmetry breaking subunit [42,43].

The remaining groups (4-11) involve complexes of two components with a single interface. Group 4 has a single member, the complex of human profilin-1 with β-actin (4), and is one of only two cases for which all bound and unbound structures were available with high sequence identity (over 94%). This case is also an entry of the protein-protein docking benchmark that we maintain [6,44-47].

Group 5 contains a single complex, the heterodimer of human Alu binding proteins SRP9 and SRP14 bound to Alu RNA (5). We ignored the RNA component during docking assuming that the two proteins could form a complex without the RNA.

Group 6 is formed by two subunits of succinate dehydrogenase (6). Although two additional subunits of the enzyme were present in the bound structure [48], they were ignored in our calculations.

Group 7 is the complex of the GTPase Rab14 with a Rab GDP dissociation inhibitor (7). For the bound structure, we used the complex of the inhibitor with the prenylated YPT1 GTPase, having sequence identities of about 50% with the target.

Group 8 consists of the human ElonginB-ElonginC complex (8) [49]. The bound structure includes the VHL tumor suppressor, but it was not included in the docking since VHL only contacts ElonginC and not ElonginB.

Groups 9, 10, and 11 represent different interfaces of the five-protein barrel assembly machinery (BAM) complex, responsible for the proper assembly of β-barrel proteins into the outer membrane of E. coli. The single member of group 9 contains the two interacting components BamC and BamD (9). BamC contains a 73-residue-long unstructured region essential for binding BamD (Figure 1) [50]. The unbound structure of the full-length BamC was not available, but even if it had been available, it might not have been suitable for rigid-body docking due to the unstructured region; consequently, we use the bound structure in our docking. Groups 10 and 11 are the complexes of BamA with BamB (10) and BamD (11A and 11B), respectively. The bound structures of group 11 show that two domains of BamA form separate interfaces with BamD. But since the BamA unbound structure contained only one of these domains, we considered only a single interface in this work.

Docking without cross-linking data

ZDOCK was used to predict complex structures for the 19 cases in the test set, using the component proteins as described above. We used interface RMSD between the predicted and bound structures (iRMSD) to assess the predictions. We applied an iRMSD cutoff of 5.0 Å to denote a prediction as a near-native structure, or a ‘hit’, [51-53] if this resulted in at least one hit in the ZDOCK calculation without constraints, otherwise we used 7.5 Å (Table 2). In case of distinct interfaces between the components (1, 3A-C), we assessed the docking predictions for each interface separately and also combined (claiming a prediction correct if either one of the interfaces observed in the bound structure matched the prediction). For 3A-C, whose bound complexes showed similar interfaces between the components with the differences caused by the symmetry breaking central chain, we based the assessment on the bound interface that yielded the lowest iRMSD.

Tables 2 and 3 show the results of unconstrained docking as well as docking combined with the various approaches of applying cross-linking constraints. Without constraints, more than half of the docking runs had a hit within the top 100 predictions (15 out of 23 interfaces), often within the top 10 (9 interfaces), but a top ranked hit was only found three times. These results are in line with the ZDOCK performance on the protein-protein docking benchmark [6]. Interestingly, three of the interfaces that had no hits at all within the top 100 predictions corresponded to the cases with multiple distinct interfaces (1, 3A, and 3B). We can speculate that the formation of the complex occurs in stages, in which the B:C (1) and A:D (3A, 3B) interfaces are only stable after the B:A and A:E interfaces have formed. Alternatively, the chains that we ignored during docking (indicated in green in Figure 1) may be required for the complete complex formation.

Separating cross-links by interface

In all of our calculations, we assumed that the experimental cross-link data did not include false positives. Consequently, we applied hard cutoffs to the cross-link distances calculated for the predictions and required a prediction to satisfy all cross-links. It was straightforward to apply these requirements for binary complexes 2A, 2E, 2F, and 4-11. Furthermore, 3B has two distinct interfaces but only one cross-link, so we simply focused on the interface associated with the cross-link. For cases 1, 3A, and 3C however, we had two distinct interfaces and multiple cross-links, and the cross-links for one interface may not be satisfied by the other interface of the same complex. Thus, it was essential that we assigned each cross-link to the appropriate interface. A similar issue arose for case 2B, whereby the cross-links could occur between the α-subunit and either chain of the β-subunit homodimer. In these situations, we need to group the cross-links by interface (1, 3A, and 3C) or chain (2B) so that we simultaneously apply only the cross-links that belong to the same interface or chain. Such grouping information, however, is not directly provided by the experimental cross-link data.

To group the cross-links by interface or chain, we assumed that two cross-links that are both satisfied by a docking prediction are more likely to be associated with the same interface (or chain) than with different interfaces. Thus, we defined a binary matrix with rows corresponding to the top 200 ZDOCK predictions, columns corresponding to the cross-links, and the elements set to 1 for the predictions that satisfy the cross-link (using Euclidean distance with a 30 Å cutoff) and 0 otherwise. We then calculated the correlation coefficients r between all pairs of columns to obtain a correlation matrix for each test case with each dimension corresponding to the total number of cross-links. Figure 2 depicts the nine test cases with multiple cross-links, with the r>0.1 elements shaded. For cases 1, 2B, and 3A a block-diagonal pattern arose, with each block corresponding to the cross-links that belong to the same interface or chain. Case 3C, with had many more cross-links, showed two large blocks and several smaller blocks. The large blocks indeed contain cross-links separated for the two interfaces, and the remaining cross-links may reflect higher order F1-ATP synthase oligomers that were not represented in the crystal structure [38]. The latter were ignored in the docking calculations. As expected, the cross-links formed a single group for binary complexes with multiple cross-links. We applied each group of constraints separately during the docking runs. In addition to improving docking performance, the occurrence of multiple groups of cross-links provides information on the stoichiometry and topology of the complex, which we explored further to determine the symmetry of the complexes (see below).

Filtering using Euclidean distances

We filtered out predictions with cross-linked sites further than 40 Å apart in Euclidean distance, which was based on the distance distribution in the bound structures (Table S1). This corresponds roughly to the distance between the C_β atoms of two cross-linked lysine residues when the linker is fully extended (34 bonds, assuming tetrahedral configurations). Although we did not observe distances below 12 Å in the bound structures, we did not apply a lower limit cutoff as shorter cross-links were observed for monomers and homomers [34]. Even with the loose cutoff of 40 Å, the highest ranking hit was eliminated for four docking runs— the second hit was retained for 3C and 9, and only the fourth hit for 2C and 2D. As a result, incorporating cross-link data worsened the rank of the top hit for 2C and 2D. Nevertheless, the Euclidean filter resulted in seven cases with a hit ranked as number one, a substantial improvement from the three without the filter (Table 3).

When we tightened the cutoff to 35 Å, there were more cases for which the top-ranked hit did not pass the filter and several cases for which all hits were filtered out. However, the overall results slightly worsened with the 40 Å Euclidean filter or without filtering, judged by the number of cases with at least one hit in the top 1, top 5, or top 10 predictions (Table 3). Tightening the filter further to 30 Å worsened the overall results further. Thus, a 40 Å cutoff provided the Euclidean filter with the best performance for the present data set.

Filtering using void-volume distances

Since void-volume distances for the bound structures were somewhat larger than Euclidean distances (Table S1), we increased the cutoff by 5 Å. Again, we observed a tradeoff between losing hits and improving the ranks of the retained hits, and we found the optimal cutoff to be 45 Å. The void-volume filter only performed slightly better than the Euclidean filter.

Euclidean constraints within FFT

The third type of constraint also uses Euclidean distances, but it restricts the translational search space of the docking algorithm and is therefore implemented within the fast Fourier transform (FFT) step of the ZDOCK algorithm. We tested the same cutoffs as for the Euclidean filter, and the performance was comparable with the above two post-processing approaches. The FFT-based constraint with a cutoff of 30 Å showed the best performance for the top 5 and top 10 predictions, and with 40 Å for the top 1 prediction. Overall, the FFT-implemented constraint approach with a 30 Å cutoff yielded the best performance among all method and cutoff combinations if we sum the number of hits for the top 1, 5, and 10 predictions. It is reasonable to consider the top 5 predictions in follow-up computational or experimental work, for which this approach had a success rate of 64% (14 out of the 22 interfaces for which we had cross-links). Note that without incorporating cross-link data we had to make 61 predictions for each test case to obtain the same success rate. Of the eight tests cases for which this approach failed to generate a hit among the top 5 predictions, only two were correctly predicted in the top 5 by any of the other method and cutoff combinations.

Symmetry analysis

The occurrence of multiple groups of cross-links indicates (Figure 2) that at least one component of the complex is a homo-multimer (as in the case of 2B) or the components form multiple distinct interfaces. The latter can then lead to symmetric complexes (as in the case of 3A and 3C) or a linear configuration (as in the case of 1). To differentiate these three possibilities, we integrated the experimental cross-link data with docking to predict whether a complex had symmetry and whether this symmetry could be used to improve the docking performance.

Cases 1, 3A, and 3C showed two distinct interfaces each, and we asked whether predictions for these interfaces could lead to symmetric complexes. We started with the top 5 predictions for each interface, obtained from the 30 Å FFT-implemented distance constraints. Combining these top 5 predictions yielded 25 predicted interface pairs. Starting from a monomer, we used each of these 25 interface pairs to sequentially add components, creating tetrameric (C₂), hexameric (C₃), and octameric (C₄) complex structures for each interface pair. The resulting structures were not symmetric because the docking was performed without any symmetry information. We moved the component proteins as rigid bodies to reach a symmetric structure while keeping the deformation of the predicted interfaces to a minimum (see Methods), and retained the symmetrized structure only if the iRMSD between the predicted interfaces and symmetrized interfaces did not exceed 2.5 Å, and all the cross-link constraints were still satisfied using a cutoff of 32.5 Å (the original cutoff of 30 Å increased by 2.5 Å to account for the allowed interface deviation). Figure 3 outlines the procedure for a single interface pair.

As shown in Tables S2, S5, and S8, none of the cases had a predicted interface pair that was consistent with the C₂ symmetry. Case 1 showed two symmetry-consistent interface pairs, which had the C₄ symmetry (Table S4). For 3A, six predicted interface pairs were consistent with the C₃ symmetry (Table S6), and nine with C₄ (Table S7). From the clash count (Tables S4 and S7) we see that a few octameric complexes symmetrized to unrealistic structures, but these examples also had starting interfaces that are far from the bound form. For 3C, we only found a single predicted interface pair consistent with C₃ symmetry (Table S9), and nine with C₄ (Table S10).

Table 4 integrates the results and they roughly agree with the symmetries observed in the bound structures: we found nine symmetry-consistent pairs for both 3A and 3C (C₃ in bound) and only two pairs for 1 (non-symmetric in bound). The bound structure of 3A has C₃ symmetry, and indeed the majority (six of the nine) of the interface pairs symmetrized to a hexamer. For 3C on the other hand, all nine symmetry-consistent interface pairs symmetrized to an octamer, while the bound structure is hexameric. Hexameric and octameric structures are difficult to distinguish by our algorithm because the difference between the C₃ (hexamer) and C₄ (octamer) symmetries is only 15 degrees per interface. In agreement with the bound structures, no predicted interface pairs yielded the C₂ symmetry. Our results based on these two cases suggest that combining docking with cross-linking data can reveal whether a protein-protein complex is symmetric, although the predicted fold of the symmetry is not precise.

Finally, we tested whether such a symmetry analysis could be used to improve the accuracy of complex structure prediction. For the three cases combined we considered 75 pairs of interface predictions, of which six showed both interfaces below the iRMSD cutoff for denoting a correct prediction (5.0 Å and 7.5 Å for cases from group 1 and 3, respectively), representing a success rate of 8%. Twenty of the predicted interface pairs were symmetry-consistent, and five of the corresponding symmetrized structures had interfaces below the iRMSD cutoff (Table 4). Thus the success rate increased from 8% to 25% if we only retained the symmetry-consistent interface predictions.

Data sets

We compiled the cross-linking data sets from previous work and new experiments [25,29-34,37]. In all cases, cross-linking was performed using the BDP-NHP cross-linker described extensively in previous work [30]. Briefly, cross-linker was added to live cells resuspended in phosphate buffer (170 mM KH₂PO₄, pH 8.0), the cells were lysed, and protein lysates were digested with trypsin. Cross-linked peptide pairs were fractionated by strong cation exchange, enriched with monomeric avidin beads (Thermo) and loaded onto an in-house pulled 45cm C8 reverse phase column for injection into an LTQ-Velos FT-ICR instrument. ReACT [30] was run as previously described to identify cross-linked peptides in real time, and peptides were searched using SEQUEST.

Protein-protein docking

We used our ZDOCK3.0 algorithm for the prediction of protein-protein complex structures [3,4,18-21]. ZDOCK inputs the structures of two constituent proteins and performs an exhaustive grid-based rigid-body search to predict their binary complex. The search returns an ensemble of predictions ranked using a scoring function, which includes shape complementarity, electrostatics, and the IFACE statistical potential [54,55]. Optimal results are typically obtained using X-ray crystallography structures as input, but also NMR, homology modeled, or ab initio modeled input structures can be utilized [56].

ZDOCK separates the full six-dimensional rigid-body space into a three-dimensional translational space and a three-dimensional rotational space. For each point in the rotational space, the best scoring translational pose is used as a prediction. In this study, we used 15° rotational sampling. Each docking run resulted in 3,600 predictions, which were ranked according to the docking scores.

Cross-linking constraints

We considered three methods for incorporating cross-linking constraints in ZDOCK. The first two approaches involved post-processing, or filtering the set of predictions after the ZDOCK calculation. Cross-link distances were computed based on the predicted structures, and only when the distance was below a cutoff the constraint was considered satisfied, and the prediction retained. The two filtering approaches differed in the calculation of the cross-link distances. In the first, we used Euclidean (straight-line) distances between the C_β atoms. In the second, we used void-volume distances, computed with the command line version of the Xwalk program by Kahraman, Malmström and Aebersold (v0.6, using the C_β atom as anchor and the -bb flag in addition to the default options) [36]. Xwalk is grid-based and uses a breadth-first algorithm to obtain the shortest residue-residue distance that passes only through solvent accessible space.

The third approach, FFT-based constraints, integrates the cross-linking constraints within the translational search of ZDOCK, following the algorithm presented by Xia [28]. In ZDOCK, Fast Fourier Transform (FFT) is used to generate a three-dimensional matrix that, for a point in the rotational space, contains the scores for all the (grid-based) translational coordinates. In a standard ZDOCK calculation, this matrix is then searched for the best score, and the corresponding complex structure is retained as a prediction. This structure may or may not satisfy the cross-links, hence the need for the filtering steps outlined above. In the FFT-constraint approach, on the other hand, we exploited the one-to-one correspondence between the score matrix elements and complex structures. For each score matrix element we could trivially compute its corresponding Euclidean cross-link distance, and when larger than the cutoff, marked the element as excluded. We repeated this procedure for all the cross-links observed for the complex. The subset of non-excluded elements was then scanned as usual, which yielded the best scoring complex structure that satisfied all the cross-links. The modified version of ZDOCK is available at http://zdock.umassmed.edu/software/download/.

Test set construction

We used BLAST [57], with a threshold of 30% sequence identity, to identify heteromeric complex structures in the PDB [35] for which we had cross-linking data available. We applied the following restrictions: First, the cross-linked sites needed to be part of the aligned regions. Second, the cross-linked sites needed to be resolved in the X-ray crystallography structure. The complexes were then investigated for suitability for protein-protein docking, using similar requirements as used for our protein-protein docking benchmarks [6]. For example, co-folded chains were excluded, as well as three-body (or higher order) interactions and protein-peptide like complexes. For the resulting complex list, we then searched the PDB for unbound structures. When finding bound or unbound structures that were less than 100% sequence identity with the cross-linked proteins, we used Modeller v.9.12 [57] to generate homology models, except for case 5 that had minimal sequence identity with the template, and yielded a more reasonable structure using I-TASSER [58-60]. We required at least one of the components to be available in its unbound form or as homology model based on an unbound template.

Prediction assessment

To measure the quality of a prediction we used the C_α interface root-mean-square-deviation (iRMSD), which results from superposing the predicted interface onto the bound interface [6]. The interface includes all residues that have at least one atom within 10 Å of the binding partner in the bound structure. When the bound structure had multiple interfaces, we calculated the iRMSD for each interface separately and used the lowest value.

Symmetry analysis

Starting with a single component and predictions for the two interfaces, we built multimeric structures by repeatedly adding components according to the predicted interfaces. We used PyMOL [61] for the superposition, and components were allowed to overlap if this followed from the interfaces. The resulting tetrameric, hexameric, and octameric structures were not symmetric because the predicted interfaces resulted from docking runs without any symmetry considerations. Therefore, we optimized each starting structure to the symmetric structure with the smallest deviation from the predicted interfaces, while keeping the components rigid. To achieve this, we used an optimization function that consisted of two components. The first component was the root mean square of the iRMSDs (defined above) between the predicted interfaces and the interfaces at the current optimization step. For the second component we followed the approach by Nilges [62]. We defined six centers for each of the two unique components (maximum and minimum along the three Cartesian axes) and calculated the distances between the paired centers located on different components. In a symmetric structure, the distances between centers across similar interfaces are identical. The optimization function thus contained the root mean square of the distance differences, which goes to zero at convergence. Although in principle the iRMSD component of the composite function is sufficient, we found that adding the specific symmetrizing component according to Nilges significantly improved the convergence behavior. For the optimization, we used steepest descent, with numerically calculated gradients. To speed up the optimizations we performed several thousand steps with only the interface similarity component, followed by additional steps using the full composite optimization function until a symmetric structure was obtained. The resulting function value, which is the root mean square of the iRMSDs (as the symmetrizing component was zero for the optimized structure) was denoted the symmetrizing iRMSD and used to measure how much the predicted interfaces deviated from symmetry.