Understanding Docking Complexes of Macromolecules Using HADDOCK: The Synergy between Experimental Data and Computations

Equipment

1. A computer with internet access and a web browser (preferably Chrome, Safari or FireFox)
2. A visualization program such as UCSF Chimera, or PyMOL
   

Software

1. HADDOCK 2.2
   HADDOCK 2.2 (High Ambiguity Driven protein-protein DOCKing) is freely available to non-profit users via its web interface here (van Zundert et al., 2016).
   Access does however require registration, which is free for non-profit users.
2. Program for protein solvent accessibility calculation
   The NACCESS program is freely available here for researchers at academic and non-profit-making institutions. NACCESS is used to identify solvent accessible residues. It is worth noting that Naccess program requires Linux/Unix systems. A free alternative to NACCESS is freeSASA, freely available here.
   
   

Data requirements

1. 3D coordinates of the two components of the protein complex, preferably in the bound conformation, provided in PDB or mmCIF format.
   Note: The 3D coordinates in PDB format of HCN2 CNBD and TRIP8b_nano are available from the following link: (LINK TO THE PDB). Those structures were solved by solution NMR as described in Saponaro et al., 2018. For docking purposes, the unfolded regions of both CNBD and TRIP8b_nano structures, which are not involved in the binding, were removed. This allows to simplify the docking calculation and avoids spurious interactions between the disordered regions of CNBD and TRIP8b_nano. In the case of TRIP8b_nano the unfolded C-terminal tail was removed (I272-N275), whereas in the case of CNBD the unfolded N-terminal region was removed (E521-F533).
2. Definition of Active and Passive residues
   Active residues can be identified by using a large variety of experimental data, including solution state NMR and biochemical protein-protein interaction data, as described in this protocol. Passive residues can be manually assigned by the operator, as performed in this protocol, or automatically defined by HADDOCK. Indeed, the program uses a 6.5 Å distance cutoff from the heavy atoms of the active residues to define the passive ones. Active and Passive residues are incorporated in the input page as comma-separated list of residue numbers (see Procedure).
   
   Note: In the representative case described in this protocol, we define the following active and passive residues for CNBD and TRIP8b_nano respectively:
   
   CNBD:
   Active residues: A546, N547, R650, E653, T654, I657, D658, R659, R662, I663, K665.
   Passive residues: P543, D549, P646, R649, G664 and K666. Manually defined by the operator.
   
   The above residues are defined based on the NMR Chemical Shift Perturbations (CSP) analysis. CSP highlights the residues of the interacting protein whose HN NMR signal features a significant chemical shift perturbation upon interaction with the peptide. A necessary prerequisite to conduct CSP is the availability of the assignment of the NMR backbone signals of the protein of interest.
   
      The CSP (ΔHN) is given by the following equation:
   
   
   ΔHN={((HNfree-HNbound)²+(((Nfree-Nbound))⁄5)²)⁄2)}^1⁄2
   
   
   where HNfree, Nfree and HNbound, Nbound, are the chemical shift of amide protons (HN) and amide nitrogen (N) in the peptide-free and bound state respectively. In order to select a meaningful set of residues involved in the binding to the peptide, a cutoff which corresponds to the average CSP value plus one standard deviation should be selected.
     While more residues might be affected by the binding, we assume that only solvent-expose residues will form contacts. Accordingly, the residues identified by CSP need to be filtered based on their solvent accessibility (either main chain or side chain relative accessibility should be typically > 40-50%). Solvent accessibility can be calculated with the program NACCESS. Passive residues are defined as the residues close in space to active residues and with at least 50% of solvent accessibility. The same procedure here reported is also valid for protein-protein docking calculations and can be applied for both interacting proteins.
   
   TRIP8b_nano:
   Active residues: we initially define as active residues a solvent exposed (≥ 50% of solvent accessibility) stretch of four residues (E239, E240, E241, F242, E243) previously identified as crucial for the binding to the CNBD (Santoro et al., 2011).
     Passive residues: Since TRIP8b_nano is a small peptide (40 amino acids), all other solvent accessible (≥ 50% of relative solvent accessibility) residues of the peptide are defined as passive.
3. Experimental data (rational mutagenesis and protein-protein interaction assay). Experimental data can be used to validate the HADDOCK output.
   

Procedure

Submitting a HADDOCK calculation:

1. HADDOCK calculation can be launched from the web-sever available at: https://haddock.science.uu.nl/services/HADDOCK2.2.
   Note: In the representative case described in this protocol, the “Guru interface” is used to run HADDOCK since it allows to specify/modify more parameters to fine tune the docking settings. Access to this interface does require “guru” level, which can be requested by users in their own registration page.
2. Fill in the HADDOCK input page (Figure 2):
   
   1. Provide a name for the job.
      
   2. Provide the input data, i.e., the 3D coordinates and the lists of the active and passive residues as a comma-separated list of residue numbers for each components of the complex.
      
   3. It is possible also to personalize the docking run by modifying several default parameters, if needed:
      
      1. Histidine protonation state for each component.
         
      2. Semi-Flexible and Fully Flexible segments.
         
      3. Additional restraints such as, for example, distance restraints, dihedral and hydrogen bonds restraints, noncrystallographic symmetry restraints, residual dipolar coupling.
      4. Sampling parameters.
         
      5. Parameters for clustering.
         
      6. Energy and interaction parameters.
         
      7. Scoring parameters.
         
      8. Advanced sampling parameters.
         
      9. Solvated docking parameters.
         
      10. Analysis parameters.
          
      
      
      Figure 2. Example view of an input page of the HADDOCK2.2 web-server
      
      Note: A detailed description of the various parameters of the HADDOCK scoring function and of the details about the HADDOCK molecular docking process can be found in the online manual here.
   4. By default, HADDOCK generates 1,000 models in the initial rigid body docking. This parameter can be changed by unfolding the tab named “Sampling parameters” (Figure 3).
      Note: In the specific case of CNBD-TRIP8b_nano complex the number of structures is increased by a factor of five.
   5. The default number of structures further subjected to the semi-flexible simulated annealing (it1) and final refinement in water (itw) is 200.
      Note: In the specific case of CNBD-TRIP8b_nano complex this parameter is increased by a factor of two (Figure 3, “Sampling parameters” tab).
   6. The final water-refined models are clustered and ranked based on HADDOCK score value. Clustering is performed based on pairwise Root Mean Square Deviations (RMSD, the default value is set to a 7.5Å cutoff) and a minimum number of members to define a cluster (default value is set to four) (Figure 3, “Parameter for clustering” tab). The clusters are then ranked based on the average HADDOCK score of their top four members. The RMSD cut-off value is crucial for a proper clustering of the structures and different clustering cutoffs might be tried.
      Note: In the specific case of the CNBD-TRIP8b_nano complex the RMSD cut-off was kept to the default value.
      
      
      Figure 3. The sampling parameters and the parameters for clustering can be changed by the user in the Guru Interface level
      
      
   7. It is recommended to set the “Advanced sampling parameters” (Figure 4) according to the specific protein-peptide docking calculation following the settings reported in Trellet et al., 2013.
      Note: In the representative case described in this protocol, we have quadrupled the default number of molecular dynamics integrations steps for the rigid body dynamics at high temperature and the subsequent cooling stages with flexibility at the interface (from 500/500/1,000/1,000 to 2,000/2,000/4,000/4,000) (Figure 4, “Advanced sampling parameters” tab). The decision to increase the values is based on the fact that for one of the two component of the complex, TRIP8b_nano, the number of the active residues was significantly limited (4 residues out of 40 forming TRIP8b_nano). This can affect the sampling of the conformational space during the rigid body docking calculation phase. Furthermore, because CNBD contains a disordered C-terminal tail, which is necessary for the binding (Saponaro et al., 2014), the number of MD steps is also increased to allow the protein to adapt its conformation during the flexible refinement stage.
      
      
      Figure 4. Advanced sampling parameters used for CNBD-TRIP8b_nano complex
      
      
      
3. Provide both username and password to submit a job.
   

Data analysis

1. Submission data
   Once all data have been uploaded and the run submitted, the web-server offers the option to download a parameter file and provides a link to the results page. This link is also e-mailed to the user.
2. Analysis of the docking results
   
   1. The result page of a successful docking run presents the statistics for the top10 clusters ranked by HADDOCK. Those are ranked based on the average HADDOCK score of the top4 member of each cluster, which is a weighted sum of electrostatics, van der Waals, restraints energy, buried surface area and empirical desolvation energy terms (Fernandez-Recio et al., 2004). The page indicates how many water-refined models could be clustered and displays at the bottom a graphical overview of the results.
      Note: For the specific case of CNBD-TRIP8b_nano complex calculation, seventeen clusters are obtained and ranked according to their HADDOCK score.
      
   2. For each cluster, detailed statistics are displayed, representing the average values calculated over the top four best-scoring structures within each cluster. The clusters require to be also manually screened for the relative orientation of the two components of the complex on the basis of available experimental data.
      Note: In the specific case of CNBD-TRIP8b_nano, Double Electron-Electron Resonance (DEER) data are available (DeBerg et al., 2015). They provide nanometer-scale distances between unpaired electrons in paramagnetic systems or tagged with paramagnetic probes. DEER data were used to screen the relative orientation of TRIP8b_nano and CNBD and select two out of the seventeen clusters obtained (Figure 5). The two selected clusters lay within the 30% top ranked clusters. Table 1 summarizes the HADDOCK statistics of the two selected clusters.
   
   
   Figure 5. Representative poses of cluster families obtained from the first docking calculation. CNBD is colored in gray; TRIP8b_nano is colored in light blue in cluster 1 and in green in cluster 2. Residues forming cluster-specific contacts are shown and labeled (modified from Saponaro et al., 2018 with permission).
   
   Table 1. HADDOCK Score and energetic values for the selected clusters
   
   
   
3. Validation of the generated models by mutagenesis
   Rational mutagenesis studies, coupled with an appropriate protein-peptide assay, can be used to validate a cluster-specific contact and eventually perform a second round of modelling introducing as active those residues whose mutations abolish the interaction.
   
   Note: This approach was used for determining the structural model of CNBD-TRIP8b_nano complex.
   In order to select residues for mutagenesis we examine the intermolecular contacts formed in each of the two selected clusters. This reveals that:
   
   1. Residues E264 or E265 of TRIP8b_nano are found to interact with residues K665 or K666 of CNBD in both clusters (see Figure 5).
      
   2. In cluster 1, TRIP8b_nano is oriented in such a way that allows E239-E243 residues to interact with R650 of CNBD; whereas, in cluster 2, TRIP8b_nano is oriented, in almost all the conformers of the cluster, in such a way that CNBD R650 cannot developed contacts with TRIP8b_nano.
   Rational mutagenesis is performed to validate the above mentioned cluster-specific protein-peptide contacts a posteriori. The effect of the mutations on the complex formation was tested by means of Isothermal Titration Calorimetry (ITC) assay (Saponaro et al., 2017 and 2018). Reverse charge mutations K665E or K666E of CNBD, and E264K or E265K of TRIP8b_nano each significantly reduced the binding affinity between the two proteins, thus indicating that this set of residues are involved in the CNBD - TRIP8b_nano complex formation. Moreover, the reversion of the charge of R650 (R650E) causes a significant reduction in binding affinity for TRIP8b_nano (for details see Saponaro et al., 2018).
4. Refining the model by a second round of docking
   On the basis of the experimental observations described above, we performed a second molecular docking calculation including E264 and E265 as additional active residues for TRIP8b_nano. The docking setup was similar to the first docking calculation described above.
      From this second HADDOCK calculation fourteen clusters are obtained. Four of these are consistent with the DEER and mutagenesis data (Figure 6). In order to validate these models, we performed a second round of mutagenesis based again on identified unique contacts. In particular, in cluster 1 D252 of TRIP8b_nano, located at the junction between helix N and helix C, contacts residue N547 of the CNBD, while this contact is absent in the other clusters.
   
   
   Figure 6. Representative poses of cluster families obtained from the second docking calculation. CNBD is colored in gray in all the four poses; TRIP8b_nano is colored in blue in cluster 1, in magenta in cluster 2, in green in cluster 3 and in yellow in cluster 4. D252 of TRIP8b_nano is shown in all the four poses; N547 of CNBD is shown in the pose representing cluster 1 as in this cluster N547 contacts D252. Helices N and C of TRIP8b_nano are labeled.
   
   N547 in CNBD was mutated into aspartate (N547D) to generate an electrostatic repulsion with the supposed partner residue D252, while the carboxyl group in D252 of TRIP8b_nano was removed by mutating it into an asparagine (D252N). Both mutations reduced binding to TRIP8b (for details see Saponaro et al., 2018). These results suggest that Cluster 1 can be chosen as the bona fide model for the CNBD-TRIP8b_nano complex.
      It is worth noting that cluster 1 represents also the top-ranking cluster for energetic, scoring function and number of structures populating the cluster (Table 2). This is explained by the fact that the second docking calculation was further implemented with more ambiguous interaction restraints, which significantly improved HADDOCK results.
   
   Table 2. HADDOCK Score and energetic values for the selected clusters
   
   

Notes

Tutorials and video describing HADDOCK program and its use for various scenarios can be found at http://www.bonvinlab.org/education.

Acknowledgments

This protocol was adapted from Saponaro et al. (2018). A.M.J.J.B. acknowledges financial support from the BioExcel CoE (www.bioexcel.eu), a project funded by the European Union Horizon 2020 program under grant agreements 675728 and 823820.

Competing interests

There are no conflicts of interest or competing interest.

References

1. DeBerg, H. A., Bankston, J. R., Rosenbaum, J. C., Brzovic, P. S., Zagotta, W. N. and Stoll, S. (2015). Structural mechanism for the regulation of HCN ion channels by the accessory protein TRIP8b. Structure 23(4): 734-744.
2. Dominguez, C., Boelens, R. and Bonvin, A. M. (2003). HADDOCK: a protein-protein docking approach based on biochemical or biophysical information. J Am Chem Soc 125(7): 1731-1737.
3. Fernandez-Recio, J., Totrov, M. and Abagyan, R. (2004). Identification of protein-protein interaction sites from docking energy landscapes. J Mol Biol 335(3): 843-865.
4. Groß, C., Saponaro, A., Santoro, B., Moroni, A., Thiel, G. and Hamacher, K. (2018). Mechanical transduction of cytoplasmic-to-transmembrane-domain movements in a hyperpolarization-activated cyclic nucleotide–gated cation channel. J Biol Chem 293(33): 12908-12918.
5. Porro, A., Saponaro, A., Gasparri, F., Bauer, D., Gross, C., Pisoni, M., Abbandonato, G., Hamacher, K., Santoro, B., Thiel, G. and Moroni, A. (2019). The HCN domain couples voltage gating and cAMP response in hyperpolarization-activated cyclic nucleotide-gated channels. Elife 8: e49672.
6. Santoro, B., Hu, L., Liu, H., Saponaro, A., Pian, P., Piskorowski, R. A., Moroni, A. and Siegelbaum, S. A. (2011). TRIP8b regulates HCN1 channel trafficking and gating through two distinct C-terminal interaction sites. J Neurosci 31(11): 4074-4086.
7. Saponaro, A. (2018). Isothermal titration calorimetry: a biophysical method to characterize the interaction between label-free biomolecules in solution. Bio-protocol 8(15): e2957.
8. Saponaro, A., Cantini, F., Porro, A., Bucchi, A., DiFrancesco, D., Maione, V., Donadoni, C., Introini, B., Mesirca, P., Mangoni, M. E., Thiel, G., Banci, L., Santoro, B. and Moroni, A. (2018). A synthetic peptide that prevents cAMP regulation in mammalian hyperpolarization-activated cyclic nucleotide-gated (HCN) channels. Elife 7: 35753. 
9. Saponaro, A., Pauleta, S. R., Cantini, F., Matzapetakis, M., Hammann, C., Donadoni, C., Hu, L., Thiel, G., Banci, L., Santoro, B. and Moroni, A. (2014). Structural basis for the mutual antagonism of cAMP and TRIP8b in regulating HCN channel function. Proc Natl Acad Sci U S A 111(40): 14577-14582.
10. Saponaro, A., Porro, A., Chaves-Sanjuan, A., Nardini, M., Rauh, O., Thiel, G. and Moroni, A. (2017). Fusicoccin activates KAT1 channels by stabilizing their interaction with 14-3-3 proteins. Plant Cell 29(10): 2570-2580.
11. Trellet, M., Melquiond, A. S. and Bonvin, A. M. (2013). A unified conformational selection and induced fit approach to protein-peptide docking. PLoS One 8(3): e58769.
12. van Zundert, G. C. P., Rodrigues, J. P. G. L. M., Trellet, M., Schmitz, C., Kastritis, P. L., Karaca, E., Melquiond, A. S. J., van Dijk, M., de Vries, S. J. and Bonvin, A. M. (2016). The HADDOCK2.2 Web Server: User-Friendly Integrative Modeling of Biomolecular Complexes. J Mol Biol 428(4):720-725.
13. Weißgraeber, S, Saponaro, A., Thiel, G. and Hamacher, K. (2017). A reduced mechanical model for cAMP-modulated gating in HCN channels. Sci Rep 7:40168.