The ClusPro web server for protein-protein docking

The PIPER docking algorithm

As stated, the ClusPro server is based on PIPER, which performs the sampling. The center of mass of the receptor is fixed at the origin of the coordinate system, and the possible rotational and translational positions of the ligand are evaluated at the given level of discretization. The rotational space is sampled on a sphere-based grid that defines a subdivision of a spherical surface in which each pixel covers the same surface area as every other pixel²⁸. The 70,000 rotations we consider correspond to about 5 degrees in terms of the Euler angles. The step size of the translational grid is 1 Å. It is easy to see that for an average size protein this amounts to sampling 10⁹–10¹⁰ conformations.

PIPER represents the interaction energy between two proteins using an expression of the form E = w₁E_rep + w₂E_attr + w₃E_elec + w₄E_DARS, where E_rep and E_attr denote the repulsive and attractive contributions to the van der Waals interaction energy, and E_elec is an electrostatic energy term. E_DARS is a pairwise structure-based potential constructed by the Decoys as the Reference State (DARS)²⁹ approach, and it primarily represents desolvation contributions, i.e., the free energy change due to the removal of the water molecules from the interface¹⁶. The coefficients w_1, w_2, w₃, and w₄ define the weights of the corresponding terms, and are optimally selected for different types of docking problems. Unless specified otherwise in Advanced Options, the server generates four sets of models using the scoring schemes called (1) balanced, (2) electrostatic-favored, (3) hydrophobic-favored, and (4) van der Waals + electrostatics (see Table 1). As will be described, for complexes that do not fit into these categories and are classified as “others” in the protein-protein docking benchmarks^22–25, ClusPro generates structures using three different coefficient sets (Table 1). To understand the magnitude of these coefficients we note that E_rep and E_attr, while defined on the grid, are scaled to the van der Waals energy¹⁶, and hence w₁ < 1.0 and w₂ <1.0 yield “softening” of both repulsive and attractive van der Waals terms. Such softening is necessary, since the bound and unbound conformations of the proteins to be docked generally differ, in some cases substantially. The DARS potential has been parameterized on a set of complexes that included a large number of enzyme-inhibitor pairs and multi-subunit proteins, and hence the resulting potential assumes good shape and electrostatic complementarity²⁹. Since E_DARS is scaled to the magnitudes of protein-protein binding free energies, w₄ =1.0 is the “neutral” choice. E_elec is represented by a truncated and smoothed Coulomb expression, also defined on the grid¹⁶. Since the charges are given as multiples of the electronic charge and the interatomic distance is measured in Å, obtaining the electrostatic interaction energy in kcal/mol using the Coulomb expression and these units³⁰ would require the coefficient w₃ =332. However, the truncation and smoothing reduce the magnitude of E_elec relative to the value based on the Coulomb expression, and hence we use w₃=600 in the balanced option of the parameter set. This set was shown to generally provide very good results for enzyme-inhibitor complexes. If it is known or assumed that the association of two proteins is mainly driven by their electrostatic interactions, then we select results obtained by the electrostatic-favored weights, in which the weight of the electrostatics is doubled relative to the balanced energy expression. In contrast, for complexes primarily stabilized by hydrophobic interactions we use the hydrophobic-favored potential, which thus doubles the weight of the desolvation term. In the fourth option, van der Waals + electrostatics, the pairwise potential E_pair is not used. The need for this option occurs for proteins that are very different from the ones used for the parameterization of E_DARS. The selection of the parameters in the Others Mode will be discussed further in this protocol. We recognize that the weights in the PIPER energy function are somewhat arbitrary, and were selected to optimize the success rates for different classes of protein complexes. However, as described below, the most likely models of the complex are selected on the basis of cluster population rather than energy value. In fact, PIPER does not aim to estimate the true interaction energy, and the score provided by ClusPro should not be considered as a measure of binding affinity.

Model selection by cluster population

The second step of ClusPro is clustering the lowest energy 1000 docked structures using pairwise IRMSD as the distance measure.^12,21 We calculate IRMSD values for each pair among the 1000 structures, and find the structure that has the highest number of neighbors within 9Å IRMSD radius. The selected structure will be defined as the center of the first cluster, and the structures within the 9Å IRMSD neighborhood of the center will constitute the first cluster. The members of this cluster are then removed, and we select the structure with the highest number of neighbors within the 9Å IRMSD radius among the remaining structures as the next cluster center. These neighbors will form the next cluster. Up to 30 clusters are generated in this manner. After clustering we minimize the energy of the retained structures for 300 step with fixed backbone using only the van der Waals term of the Charmm potential³¹. The minimization removes the steric overlaps but generally yields only very small conformational changes. In its basic operation ClusPro outputs the structures at the centers of the 10 most populated clusters.

Considering the centers of the large clusters of low energy structures rather than simply low energy structures appears to be unique to ClusPro,³² and may implicitly account for some entropic effects.³³ We have recently shown that under fairly natural assumptions the cluster populations are proportional to cluster probability. This argument is based on considering the approximate partition function Q = Σ_j exp(-E_j /RT), where E_j is the energy of the jth pose, and the summation includes all structures generated by the docking. Similarly, for the kth low energy cluster the partition function can be approximated by Q_k = Σ_j exp(E_j /RT), where now we consider only the structures in the cluster. In terms of these approximate partition functions the probability of the kth cluster is given by P_k = Q_k/Q. Furthermore, within the low energy cluster we introduce the approximation E_j=E, i.e., we assume that that the energy values are the same for all structures in the cluster, because the energies are from a narrow range that is comparable to the uncertainty of the energy function. Due to this assumption the probability of the cluster is given by P_k=exp(-E/RT)×N_k/Q, and thus it is proportional to the number N_k of the structures in the kth cluster. Therefore we suggest finding the most populated clusters rather than the lowest energy structures and using the centers of the clusters as putative models. Although this approach to model selection is somewhat unusual and has not been widely adopted, we believe that the consistent success of ClusPro in the CAPRI experiment^34–38 supports its use (see the section Comparison to other methods for details on the performance of ClusPro). Since it was requested by a number of users, ClusPro also outputs the PIPER energies of cluster centers, as well as the lowest PIPER energy within each cluster. However, since these values do not include entropic contributions and the weights of the energy components are selected to yield near-native structures rather than correct thermodynamics of binding, the PIPER energy does not provide valid information on binding free energy, and should not be used for ranking the models. We repeatedly emphasize that the latter is based on cluster population rather than cluster energy score.

Applications of the method

The usage of the ClusPro server has been growing beyond our expectations. By June 2016 the server had over 17,000 users (among them over 7,000 registered, which is not required), and completed more than 172,000 docking calculations, currently adding about 3,000 per month. Models built by the server have been reported in over 400 publications. These statistics demonstrate that there is substantial need for protein-protein docking, and that implementing our algorithms as a server expanded the availability of the method to parts of the research community without experience in computational structural biology. The large number of publication also implies that we can follow how the server is employed. We describe here some of the typical applications, each illustrated by examples.

Experimental Design

We describe here a number of advanced options for use in ClusPro as follows.

Comparison to existing methods

Table 2 shows a classification of direct protein-protein docking methods based on the amount of the information that is required in addition to the structures of the component proteins⁸. Each class of methods has its own strengths and limitations. Global methods are the most useful when no a priori information on the complex is available, and hence the entire 6D conformational space must be sampled. Since such methods use rigid body approximation, they allow only for moderate conformational change upon binding, which is a major limitation. The medium-range methods such as RosettaDock⁶⁷ and ATTRACT⁶⁸ can sample only selected regions of the conformational space and hence require some information on the putative complex, but the Monte Carlo or similar search algorithms can account for some level of flexibility, primarily in side chains. Therefore these methods are particularly useful when side chain conformations are very uncertain, for example when one of the component protein structures is a homology model, but some information is available on the structure of the complex to identify the region of interest in the conformational space. Finally, restraint based docking, exemplified by the program HADDOCK (High Ambiguity Driven biomolecular DOCKing)^69,70, incorporates interaction restraints into the scoring function to guide the search toward regions of the conformational space in which the restraints are satisfied. The method can work very well if a sufficient number of correct restraints are available, but performance may deteriorate without them.

The quality of docking methods is continuously monitored by CAPRI (Critical Assessment of Predicted Interactions), the ongoing communitywide experiment devoted to protein docking⁷¹. In the CAPRI challenge, initiated in 2003, participating research groups and automated servers are given prediction targets, each being an unpublished experimentally determined structure of a protein-protein complex. Given the atomic coordinates of the component proteins or of their homologues, the participants have to model the complexes. Servers must submit results within two days, whereas human predictor groups have several weeks and can use any available information. Each group can submit ten predictions for each target. The submitted models are evaluated by independent assessors. For each group or server, the assessors report the number of targets with acceptable or better quality predictions, and also note the number of targets for which highly accurate (***) or medium accuracy (**) models were submitted. These categories have been defined on the basis of the fraction of native contacts, the backbone root mean square deviation of the ligand (LRMSD) from the reference ligand structure after superimposing the receptor structures, and the backbone RMSD of the interface residues (IRMSD). The calculation of these measures and the exact definitions of categories were given in the first CAPRI evaluation paper⁷². Although the CAPRI quality score is based on a number of characteristics including several RMSD measures and the predicted fraction of native contacts, for simplicity we focus on LRMSD, and note that for the highly accurate, medium accuracy, acceptable, and incorrect models the ligand LRMSD is given by LRMSD < 1Å, 1Å < LRMSD < 5Å, 5Å < LRMSD < 10Å, and LRMSD > 10Å, respectively.

Since its release in 2004, ClusPro has consistently been the best automated server in the CAPRI challenge. Table 3 shows the performance of the four best servers based on the results of the last three CAPRI evaluation meetings in 2009³⁶, 2013³⁷, and 2016. These results confirm that the global systematic sampling of the entire conformational space as performed by ClusPro is useful when essentially no a priori information on the structure of the target complex can be used, which is generally the case for servers that must submit results within two days, directly produced by the server. The only freedom we had using ClusPro was the choice of the energy coefficient set. The comparison of different methods is more complex if the performance of both “human” predictor groups and servers is considered (Table 4). Because predictors can use additional information, medium-range methods such as ATTRACT by Zaccharias^68,73, SWARMDOCK by Bates⁷⁴, and FRODOCK by the Chacon group⁷⁵, recently combined with a new scoring scheme by Guerois⁷⁶, performed well. The restraint based HADDOCK method by Bonvin and co-workers also provided good predictions for many targets^69,70. However, good results were also obtained by several methods that, similarly to ClusPro, perform global search assuming rigid proteins, including ZDOCK^53,77, GRAMM⁷⁸, and PatchDock⁵⁰. Our own group’s “human” predictions are based on running ClusPro, followed by refinement and selection of most likely clusters”⁷⁹, in some cases involving Monte Carlo minimization runs either using RosettaDock⁶⁷ or our own implementation of the Monte Carlo algorithm⁸⁰. For a few targets the refinement improved medium accuracy structures into highly accurate ones, but generally the impact was moderate, and in 2016 we actually lost an acceptable prediction. According to Table 4, the relative performance of ClusPro has been improving over the years, and based on the 2016 results the server appears to be competitive with the best “human” predictor groups. We note that the PIPER scoring function has been recently used in an algorithm that expands FFT-based sampling to five rotational coordinates⁸¹. Working in a space with spherical coordinates defined as a manifold, the new 5D algorithm retains the accuracy of PIPER for globular proteins while providing at least 10-fold speedup. However, moderate loss of accuracy may occur for proteins that are very elongated along some direction. An additional advantage of the method is that increasing the number of correlation function terms is computationally inexpensive, which enables using even more complex energy functions as well as accounting for any number of pairwise distance constraints⁸¹.

Limitations

PIPER performs rigid body docking in the 6D space of rotations and translations. The rigid body assumption is a good approximation for several classes of protein complexes, e.g., for most enzyme-inhibitor complexes (Supplementary Table 1). In fact, PIPER uses a “soft” potential that allows for certain steric overlaps. However, the protein-protein docking benchmark set^22–25 also includes a number of “difficult” targets, i.e., complexes in which at least one of the proteins substantially changes backbone conformation upon binding. Neither PIPER nor any other rigid body method can produce acceptable docked structures for such complexes. However, we note that extending ClusPro to docking short flexible peptides to proteins is in preparation.

As discussed, unless requesting the “other” mode, ClusPro yields four sets of docked structures using the scoring schemes called (1) balanced, (2) electrostatic-favored, (3) hydrophobic-favored, and (4) van der Waals + electrostatics. In the “other” mode ClusPro generates structures using three different scoring schemes. The differences in the weighting of energy coefficients represent real differences in the biophysical forces that dominate interactions between two proteins, as the association in some complexes is driven primarily by electrostatic interactions, while in others desolvation of hydrophobic interfaces may be the main driving force. Unfortunately at present we are unable to perform automated selection of the best scoring scheme. Thus, it is left to the user to make such selection based on the assumed properties of the particular complex. If no such information is available, we suggest using the balanced option, and the “other” mode for signal transduction complexes. It is generally a good sign if the selected parameters lead to a cluster of low energy structures that is substantially more populated than the others. It is also useful if different parameter sets yield the same model. Once a model is selected, the user can also explore whether the assumption used was correct by analyzing the properties of the interface. Nevertheless, we recognize that the problem of parameter selection is not resolved, and further research is needed.

Even assuming that a parameter set is selected, ClusPro returns ten clusters of low energy structures, the center of each cluster representing a putative model of the complex. Again, the existence of a large cluster provides some level of certainty, but generally additional information is needed for reliable model selection. As mentioned, we found over 400 publications that reported models constructed by ClusPro. In many applications the models were validated by experimental techniques, including site-directed mutagenesis with NMR, calorimetry, FRET, or surface plasmon resonance^44,82–101, cross-linking^102–106, various spectroscopic methods^107–113, X-ray scattering¹¹⁴, electron self-exchange reaction¹¹⁵, radiolytic protein footprinting with mass spectrometry^39,40, hydrogen/deuterium exchange¹¹⁶, or intermolecular Nuclear Overhauser Effect restraints⁴⁷. ClusPro results can also be used to design low resolution and hence cost effective validation experiments. We believe that such combination of computational and experimental tools is the most meaningful use of docking.

In addition to the above fundamental limitations we have several technical shortcomings that we hope remove in the future. Some of these shortcomings are as follows: (1) At this point the molecules to be docked can include only standard amino acids, nucleic acids to define an RNA molecule as the receptor, and heparin as the ligand. Thus, the docking cannot account for the presence of co-factors or other small ligands that can be important in the modulation of protein-protein interactions. (2) Similarly, we are unable to account for non-standard or modified amino acids, e.g. phosphorylation, another important factor that can affect the interactions. (3) Currently ClusPro can create only dimers and trimers in the Multimer Docking mode, and does not provide an option for constructing more complex multimers. (4) There are two limitations on size of the proteins to be docked. First, uploaded files cannot be more than 10MB each. Second, the total grid size cannot be more than roughly 40×10⁶ ų, which corresponds to a cube of ~350 Å in each dimension. This grid has to contain both the receptor and ligand in each potential relative orientation of the two molecules.

Input data for docking TIMING ~1 to 2 min

• 1| Locate the server at https://cluspro.org. ClusPro can be used without a user account or with a user account if you have an educational or governmental email address. If you decide to use the server without a user account, click ‘Use the server without the benefits of your own account’. To create an account, register on the ClusPro server website. A password will be sent to the email address and can be later changed by clicking on the label Preferences. You can also request an e-mail to be sent to you when your job is completed. If you already have a username and password, fill in the boxes and click Login. Either with or without a user account you will get to the server home screen (Fig. 2). From this page you will be able to submit a new docking job and select a number of Advanced Options.
  

  Open in a separate window

  FIGURE 2

  Home screen of the ClusPro server. The example shows submission of the job to dock soybean trypsin inhibitor (ligand, PDB ID 1BA7) to the X-ray structure of porcine trypsin (receptor, PDB ID 1QQU).

  CRITICAL STEP Users who run ClusPro without an account will have their results publicly accessible.

• 2| (Optional) Provide a job name for this submission. If you choose to leave this blank, a unique ID will be created for this field.
• 3| Select the type of the computer the job will use. Selecting cpu will lead it to use computer clusters at the Massachusetts Green High Performance Computing Center (MGHPC). Selecting gpu will lead it to use the graphic computing units, also located at MGHPC. Since we generally cannot predict the actual usage on these systems, selecting cpu versus gpu has no predictable consequence for the user. Currently we use the cpu version, because ample computer time is available and usually the jobs start immediately after submission.
• 4| Input the coordinates of the receptor structure using PDB format. Only atoms of 20 standard amino acid residues and nucleotides will be retained. All HETATM (hetero atom) records, including waters, ligands, and cofactors, will be automatically removed. There are two options for inputting a structure: use Option A to import coordinates from the PDB, or option B to upload a structure.
• (A) Import coordinates from PDB
• i) Import coordinates directly from the PDB by typing the four digit PDB ID into the field PDB ID.
• (B) Upload structure
• i) Upload a structure from your computer by clicking on the Upload PDB option under the PDB ID field. Select Browse to upload a file containing a structure in PDB format.

  CRITICAL STEP At this point only structures containing standard amino acid residues and nucleotides will be docked. ATOM records of nonstandard amino or nucleic acids will cause an error.

  CRITICAL STEP Since the protein considered the receptor is placed on a fixed grid whereas the protein considered the ligand is placed on a rotating and translating grid, it is computationally advantageous to consider the larger protein as the receptor.

  CRITICAL STEP If the target is a protein-RNA complex, the RNA must be defined as the receptor rather than the ligand.

  (Optional) To see the supported RNA residue and atom names click on the label RNA.

  ? TROUBLESHOOTING

• 5| In the Chains field, enter the protein chains of the receptor that you wish to include in the docking. List chains using their chain id and separate multiple chains with a whitespace. If no chains are specified, then all chains in the PDB file will be used for docking.
• 6| Input the coordinates of the ligand structure using PDB format. Only atoms of 20 standard amino acid residues or a standard heparin structure can be used. If the ligand is a protein, you can use Option A to import coordinates from the PDB, or option B to upload a structure directly as in Step 4.

  CRITICAL STEP If you will use the Advanced Option Multimer Docking, no ligand structure needs to be specified.

  CRITICAL STEP If you will use the Advanced Option Heparin Docking, no ligand structure needs to be specified.

  ? TROUBLESHOOTING

• 7| In the Chains field, enter the chains of the ligand protein that you wish to include in the docking as in Step 5.

Selection of Advanced Options TIMING ~5 to 45 min (see below)

• 8| (Optional) At this point, it is possible to define additional advanced options depending on your requirements. See the following table for a summary of these options. These options are accessible by clicking the “Advanced Options” label.

  Step	Option	Description
  8A	Structure Modification	Removal of unstructured terminal residues
  8B	Attraction and Repulsion	Setting attraction or repulsion on selected residues
  8C	Restraints	Selection of pairwise distance restraints
  8D	Others Mode	Selecting a special scoring scheme for “other” type complexes
  8E	Antibody Mode	Selecting a special scoring scheme for docking antibody-antigen pairs
  8F	Multimer Docking	Constructing homodimers or homotrimers
  8G	SAXS Profile	Accounting for experimental SAXS data
  8H	Heparin Ligand	Global docking of a generic heparin molecule

  1. Selecting Structure Modification TIMING ~5 min

     1. Structure modification is suggested if any of the proteins has an unstructured terminal region. Click on Advanced Options and select Structure Modification. You can then request removal of unstructured terminal residues from the receptor, the ligand, or both.

        CRITICAL STEP This mode can be combined with other advanced options. However, editing the PDB files manually provides greater flexibility, e.g., the potential for removing flexible loops with uncertain structure that would interfere with docking.

  2. Selecting Attraction and Repulsion TIMING ~30 min

     1. If a priori experimental information indicates that certain residues are in the binding interface, you can influence the results of the docking by setting attraction on these residues. Alternatively, if some information indicates that certain residues remain surface accessible upon binding, you can influence the results of the docking by setting repulsion on these residues. Click on Advanced Options and select Attraction and Repulsion. You will then see the entry fields for attraction and repulsion.
     2. Attracting and repulsing residues can be selected by typing whitespace separated “chain-residue” entries, e.g., a-23 a-25 a-26 a-27, in the appropriate boxes.
     3. (Optional) Repulsing residues can also be selected by uploading a masking file that includes residues you do not want to be in the interface. Click on the Use PDB Masking File option, and select Browse to upload a pre-constructed masking file. To generate a PDB masking file, we recommend using PyMOL. Open in PyMOL the PDB file of the protein to be docked, and click on the button marked with “S” in the lower right corner to see the sequence of the protein. You can then select the residues to avoid in the interface by clicking on their one-letter amino acid code in the sequence. This will create a selection object called “sele” on the right side of your screen that holds these residues. You can also see the selection by what residues have dots placed on them in the viewer. You can then choose File->Save Molecule from the PyMOL menu. This will give you a window where you choose to save “sele”. You can then upload this PDB file as your masking file.

        CRITICAL STEP This mode can be combined with other advanced options.

  3. Selecting Restraints TIMING ~30 to 45 min

     1. Pairwise distance restraints may be available, e.g., from NMR Nuclear Overhauser Effect (NOE) experiments or from chemical crosslinking. To account for such restraints in docking, click on Advanced Options and select Restraints.
     2. To use this option you need to prepare a restraint file. There are two options for restraint file format: either the AIR (Ambiguous Interaction Restraints) format, which is also used by the HADDOCK program, or our JSON format. To create a restraint file in AIR format, you can download a sample AIR type restraint file from https://cluspro.org/examples/ja026939x_s2.txt. To provide restraints as a JSON file, you can download a sample JSON type restraint file from https://cluspro.org/examples/ja026939x_s2.json. We have developed a web based restraint set generator to help the preparation of the restraint file in JSON format (see Fig. 3), available at https://cluspro.org/generate_restraints.html.
        

        Open in a separate window

        FIGURE 3

        Screen image of the restraint file generator web tool to prepare a JSON file. The generator is available at the URL https://cluspro.org/generate_restraints.html.

        CRITICAL STEP Restraints can either be a JSON file with a .json suffix, or an AIR restraint file from HADDOCK with any other suffix.

     3. Choose the restraint file.
  4. Selecting Others Mode TIMING ~5 min

     1. The use of Others Mode is suggested for complexes that may have less perfect shape and electrostatic complementarity than the enzyme-inhibitor complexes. Examples are the “other” type complexes in the protein-protein docking benchmark set^22–25.Click on Advanced Options and select Others Mode. As described, this mode uses three different sets of energy weight coefficients and selects 500 low energy structures from each. The 1500 conformations are clustered together to select the most populated clusters.
  5. Selecting Antibody Mode TIMING ~10 to 45 min

     1. Using the Antibody Mode is suggested for antibody-antigen docking. Click on Advanced Options, and select Antibody Mode. Click on the label Use Antibody Mode. As described, this mode will use a scoring function specifically developed for antibody-antigen complexes.
     2. (Optional) Avoiding non-CDR (non-complementarity determining region) residues in the interface may improve the docking result. There are two options for specifying CDR-residues: either automatic or manual selection. For automatic masking of non-CDR residues, click on the label “Automatically Mask Non-CDR Regions”. For manual masking of non-CDR residues using the Attraction and Repulsion option, you need to click on the label Attraction and Repulsion. You will then see the entry fields for attraction and repulsion. Repulsing residues can be selected by typing them as whitespace separated “chain-residue” entries as shown in Step 8B(ii). Alternatively, repulsing residues can be selected by uploading a masking file generated as described in Step 8B(iii). In both cases it is necessary to identify the CDRs of the specific antibody for selecting non-CDR residues as repulsive or listed in a masking file. To find CDRs, you can consult the website http://www.bioinf.org.uk/abs/#cdrid, or use the annotation tools in the abYsis server at http://www.abysis.org. On the home page of the abYsis server, select Sequence Input Key Annotation, enter the sequence of either the heavy or the light chain of the antibody in FASTA format, and click Submit. Click on the label Numbering and Regions. The resulting page provides three different annotations of CDR regions (Chothia, ABM, and Kabat). Select the amino acid residues that are reliably not part of any CDR, and generate the masking file as described in Step 8B(ii).
  6. Selecting Multimer Docking TIMING ~5 min

     1. Multimer docking is special in that the receptor and ligand are the same molecule. Therefore you only need to fill in the information for your molecule in the receptor field. Similarly, if you are specifying information such as attraction and repulsion, you only need to specify this information for the receptor. Click on Advanced Options and select Multimer Docking. You can then enter either 2 or 3 subunits for either a dimer or a trimer.

        CRITICAL STEP The current multimer mode in ClusPro supports only dimers and trimers.

  7. Selecting SAXS Profile TIMING ~10 to 30 min

     1. If you have SAXS data, you can submit it to ClusPro for use during the docking process. The server will fit the theoretical profiles generated from the docked structures against the uploaded SAXS profile using the chi scores to filter the structures before clustering. To use SAXS data click on Advanced Options and select SAXS Profile.
     2. To use this option you need to prepare a file containing the SAXS profile. The file should contain the SAXS results in a 3-column format. The first column should be the angle (q) in units of 1/Å. The second column should be the scattering intensity (I). The third column is optional, and should be the experimental error if it is available. You can download a zip file with sample inputs from https://cluspro.org/examples/saxs.zip. The zip file contains two pdb files, rec.pdb and lig.pdb, as well as a sample SAXS profile file, SalGEc_new_autoRg.dat.
     3. Choose the file containing the SAXS profile.
  8. Selecting Heparin Ligand TIMING ~5 min

     1. Click on the labels Advanced Option and Heparin Ligand, and then select Use Heparin as Ligand.

Running ClusPro TIMING ~ 1 to 8 hours

• 9| Click on the DOCK button to begin the docking calculation. The status of the job can be immediately checked on the Queue page, which also shows the jobs that are ahead of yours. The jobs will run in order of their submission. Each job is listed with ID number, job name, user name, and a status update. To see a detailed Status page, click on the ID of your job. The Status page shows the job ID number, job name, job status, and submission time. See Box 1 for the list of status abbreviations. The Status page also shows small pictorial representations of the uploaded and processed input structures (Fig. 4).
  

  Open in a separate window

  FIGURE 4

  Screen image of the ClusPro Status page. The screen was obtained by clicking on the job Id on the queue page, for docking of soybean trypsin inhibitor (ligand) to the porcine trypsin (receptor). The Status page shows the job ID number, job name, job status, and submission time stamp. The page also shows pictorial representations of the uploaded and processed input structures The PDB file for the ligand includes two chains, but only chain B is used for docking. Clicking on the labels rec and lig will download the submitted coordinate files as read by the server.

  BOX 1

  STATUS UPDATES FOR CLUSPRO RUNS

  The progress of the ClusPro run can be monitored in the “Queue” tab. Note that since the timings for each step are highly dependent on the input structure, they are not provided here.

  Processing pdb files

  Downloading PDB file from the www.pdb.org web site, processing chain information, extracting the chains the user specified

  Pre-docking minimization

  Running CHARMM to add missing atoms and polar hydrogens, minimizing the added atoms in the presence of the protein

  Copying to supercomputer

  Copying the PDB files to the cluster where ClusPro will run

  Held on supercomputer

  Files are on cluster, but job is not yet submitted

  In queue on supercomputer

  Job has been submitted on the cluster, but has not started running

  Running on supercomputer

  Job has begun running on the cluster

  Clustering and minimization

  Clustering top structures, selecting representative from top ranked clusters for minimization. Cluster are ranked by cluster size

  Finalizing job

  Compresses files on server for transfer

  Copying to local computer

  Results are being copied from the cluster back to the ClusPro server

  Finished

  Everything is complete

  Computing Saxs profiles

  Preliminary to docking (Saxs Profile mode only)

  Error on local system

  Processing PDB files fails. Check error messages

  Error on supercomputer

  Job fails to run on the cluster. Check error messages

• 10| If requested, an email will be sent when the job has completed or if an error occurred (see Box 2 for a listing of possible errors and their meanings). The email provides a link to the Results page or to an error message. Click the link to get to the Results page. Alternatively, locate the results under the Results tab on the server. All user results will be listed in order of ID number.

  BOX 2

  ERROR MESSAGES AND THEIR MEANINGS

  When the calculation encounters an error, the job will be terminated and the user will receive an email with an error reason. The error will also be listed next to the corresponding job id in the “Results” tab. The error codes and explanations are listed below.

  xxxx not found in PDB

  xxxx is the four letter PDB ID. This error occurs when the computer is unable to download the entered PDB ID from the website, www.pdb.org. Usually this error occurs when the PDB ID does not exist, but it can also occur when the PDB website is down.

  Unknown residue xxx in receptor. Please remove

  xxx is the three letter amino acid code. This error occurs when a residue in an atom record is not recognizable by ClusPro and thus ClusPro does not have parameters for it. Check to make sure that the proper three letter code is used and the amino acid is one of the 20 naturally occurring amino acids.

  Unknown residue xxx in ligand. Please remove

  xxx is the three letter amino acid code. This error occurs when a residue in an atom record is not recognizable by ClusPro and thus ClusPro does not have parameters for it. Check to make sure that the proper three letter code is used and the amino acid is one of the 20 naturally occurring amino acids.

  Receptor chains must be fewer than 20 characters

  Chain specification is incorrect.

  Receptor chains must be white space separated alphanumeric characters

  Incorrect or missing receptor chain specification.

  Receptor PDB id must be 4 alphanumeric characters

  Invalid PDB id.

  Receptor file too large

  PDB file exceeds maximum allowed size. May happen if the protein has multiple domains. Some domains can be docked separately.

  Receptor file only partially uploaded

  PDB file is either too large or a network error has occurred during upload.

  Ligand chains must be fewer than 20 characters

  Chain specification is incorrect.

  Ligand chains must be white space separated alphanumeric characters

  Incorrect or missing ligand chain specification.

  Ligand PDB id must be 4 alphanumeric characters

  Invalid PDB id.

  Ligand file too large

  PDB file exceeds maximum allowed size. May happen if the protein has multiple domains. Some domains can be docked separately.

  Ligand file only partially uploaded

  PDB file is either too large or a network error has occurred during upload.

  Copy of receptor failed

  File did not transfer to computing cluster. Check input for consistency.

  Copy of ligand failed

  File did not transfer to computing cluster. Check input for consistency.

  Processing failed on receptor

  This error occurs during the initial steps when Charmm is used to add missing hydrogen atoms and minimize the structure. Usually, this occurs when the protein structure has sterically clashing atoms or the structure generally does not make physical sense in terms of bonds.

  Processing failed on ligand

  This error occurs during the initial steps when Charmm is used to add missing hydrogen atoms and minimize the structure. Usually, this occurs when the protein structure has sterically clashing atoms or the structure generally does not make physical sense in terms of bonds.

  Not enough lines in output for xxx

  Restraints are wrong or cannot be satisfied.

  Job ran out of memory on server

  Receptor and/or ligand are too large.

  Ligand repulsion file too large

  If the protein is large, the repulsion file is even larger (Attraction and Repulsion mode only).

  Ligand attraction/repulsion file only partially uploaded. Please try again

  A network error may have occurred during upload. (Attraction and Repulsion mode only).

  Ligand attraction/repulsion must be whitespace separated chain-residue

  Incorrect format (Attraction and Repulsion mode only).

  Receptor repulsion file too large

  (Attraction and Repulsion mode only).

  Receptor attraction/repulsion file only partially uploaded. Please try again

  A network error may have occurred during upload. (Attraction and Repulsion mode only).

  Receptor attraction/repulsion must be whitespace separated chain-residue

  Incorrect format (Attraction and Repulsion mode only).

  Number of multimers must be expressed as a number

  Format error (Multimer mode only).

  Number of multimers must be an integer

  Format error (Multimer mode only).

  Number of multimers must be 2 or greater

  Format error (Multimer mode only).

  Copy of heparin failed

  (Heparin Ligand mode only). Restart job.

  ‘SAXS’ profile only partially uploaded. Please try again

  File is either too large or a network error has occurred during upload (Saxs Profile mode only).

  Restraints file too large

  If the protein is large, the restraints file is even larger (Restraints mode only).

  Restraints file only partially uploaded. Please try again

  File is either too large or a network error has occurred during upload (Restraints mode only).

  Error converting AIR restraints to JSON. Please consult the help page

  Submitted restraints may not be formatted properly. A sample AIR type restraint file can be downloaded from the help page (Restraints mode only).

  Restraints file has incorrect format. Please consult the help page

  Restraints files must be in either AIR format or in JSON format. See the help page to use tool for generating restraints in JSON format and for sample AIR and JSON formatting (Restraints mode only).

  Restraints file refers to unknown residues

  Check residue names (Restraints mode only).

  CRITICAL STEP We will store the result files on the server for at least 2 months. After this time, the results may be deleted.

  ? TROUBLESHOOTING

Predicting an “others” type complex

Docking the component proteins of “other” type of complexes is the same as the basic ClusPro run, apart from selecting Others Mode in Advanced Options. As an example, we dock the X-ray structure of the cytoplasmic domain of the unphosphorylated type 1 TGF-beta receptor (PDB ID 1IAS, Chain A, 342 residues), which actually was defined here as the ligand, to the X-ray structure of the FK506 binding protein (FKBP, PDB ID 1D6O, Chain A, 107 residues). Fig. 8 shows the PyMOL visualization of the results. As in the previous example, the receptor (FKBP in this case) is shown as grey surface, the ligand lig.003.00.pdb at center of the largest cluster represented by the cyan cartoon, and the ligand extracted and superimposed from the complex (PDB ID 1B6C) is shown in magenta. We note that the IRMSD between free and bound structures of the TGF-beta receptor is 1.9 Å, and the structural differences between free and bound states are clearly seen in Fig. 8. Nevertheless, docking of the unbound structures yields a large cluster with 267 members, whereas the second largest cluster has only 100 members, again suggesting a stable native complex, occupying a well-defined energy well on the energy landscape. Indeed, the interface is predicted well, resulting in the IRMSD value of 2.96 Å, in spite of the overall conformational change in the ligand.

Open in a separate window

FIGURE 8

Screen image of exploring the results of running ClusPro in Others Mode. Results are shown using PyMOL from docking the X-ray structure of the ligand, to the X-ray structure of the FK506 binding protein (FKBP). The receptor, FKBP, is shown as grey surface, and the ligand (lig.003.00.pdb) at center of the largest cluster is shown as cyan cartoon. For comparison we superimposed the native complex (PDB ID 1B6C) on the receptor in the docked structure, and the corresponding ligand pose is shown as magenta cartoon.

Antibody-antigen docking

We docked the X-ray structure of the extracellular domain of the human tissue factor (PDB ID 1TFH) to the unbound X-ray structure of the FAB domain of the inhibitory antibody 5G9 (PDB ID 1FGN). Both the heavy and light chains were used to represent the receptor. We first used the Antibody Mode with automatic masking of non-CDR regions, selected in Advanced Options. As shown in Fig. 9, we have an acceptable model, lig.000.05.pdb, with IRMSD of 4.7 Å from the ligand position in the native complex (PDB ID 1AHW). However, this model is the center of the 6^th largest cluster with only 48 members, whereas the largest cluster has 118 members, but much larger IRMSD from the native. Thus, while docking produces a good model, without a priori information we would not have been able to identify it. We also prepared a sequence-specific masking file to place repulsion on the non-CDR residues as suggested in Step 8E(ii). The abYsis server provides three different annotations of CDR regions (Chothia, ABM, and Kabat). We put repulsion only on residues that were not considered to be in a CDR by any of these annotations. The manual masking substantially improved the results, and a model (lig.000.01.pdb) with the IRMSD of 4.8 Å moved to the 2^nd most populated cluster with 73 members. The largest cluster, with 190 members, was still a false positive, emphasizing the uncertainty of docking results for antibody-antigen complexes. Nevertheless, the example also shows that manual preparation of a sequence-specific masking file can be useful.

Open in a separate window

FIGURE 9

Screen image of the PyMOL visualization of the results of running ClusPro in Antibody Mode. We docked the X-ray structure of the extracellular domain of the human tissue factor (PDB ID 1TFH) to the unbound X-ray structure of the FAB domain of the inhibitory antibody 5G9 (PDB ID 1FGN). Both the heavy and light chains were used to represent the receptor. The center of the 6^th most populated cluster, lig.000.05.pdb, is shown as cyan cartoon, whereas the antibody is shown in surface representation. The antigen in the native complex (PDB ID 1AHW) is shown as magenta cartoon. The IRMSD between native and predicted ligand poses is 4.7 Å.

Accounting for pairwise distance restraint

As an example we present constructing the complex from the X-ray structures of the Escherichia coli glucose-specific phosphocarrier protein IIAglc (E2A, PDB ID 1F3G, considered the receptor) and the signal transducing protein HPr (PDB ID 1POH, considered the ligand)¹¹⁷. ClusPro works very well for this problem even without any restraint, as the center of the second ranked cluster, has the IRMSD of 3.78 Å from the native state of the complex (PDB ID 1GGR). As shown in in Fig. 10A, this model (cyan cartoon) is oriented similarly to the native ligand (magenta cartoon), but its position is slightly shifted. However, this second cluster has only 221 members, whereas the top cluster has 424, but the ligand structures in this cluster are rotated relative to the native structure, resulting in the IRMSD of more than 20 Å. NOE (Nuclear Overhauser Effect) measurement were available for this complex¹¹⁷, resulting in 20 intermolecular distance restraints. Docking with these restraints was also performed by HADDOCK using a set of AIRs⁶⁹ (https://cluspro.org/examples/ja026939x_s2.txt). We have prepared a JSON file using https://cluspro.org/generate_restraints.html, the on-line restrain set generator tool (see Fig. 3). The JSON file is shown at https://cluspro.org/examples/ja026939x_s2.json. Accounting for the restraints further improves the result, and the center of the largest cluster, with population 268, is shifted to 2.88 Å IRMSD. Indeed, as shown in Fig. 10B, the model (shown as blue cartoon) is now turned by a few degrees around an axis perpendicular to the middle of the receptor binding site, and this yields only small IRMSD. While accounting for the restraints improved the result, the identification of the best structure without additional information still would be difficult, as the second largest cluster has 235 members, thus it is almost as large as the top cluster.

Open in a separate window

FIGURE 10

Best models from docking the signal transducing protein HPr (PDB ID 1POH) to the glucose-specific phosphocarrier protein E2A (PDB ID 1F3G) without and with restraints. The receptor protein, E2A, is shown as grey surface. A. Ligand position at the center of the second largest cluster (shown as cyan cartoon) from docking without restraints. The position is slightly shifted relative to the native ligand binding position, shown as magenta cartoon. B. Ligand position at the center of the largest cluster (shown as blue cartoon) from docking with restraints. The ligand is now turned a few degrees around an axis perpendicular to the center of the receptor.

Accounting for experimental SAXS data

Target 58 of the CAPRI docking challenge was determining the complex between the salmon cold-active goose-type lysozyme and the e. coli PliG lysozyme inhibitor using SAXS data as restraints¹¹⁸. A model of salmon lysozyme was built using the MODELLER v9.0¹¹⁹ program using the template of black swan goose lysozyme (PDB ID 1GBS, 57.8% sequence identity). Aromatic residues (Tyr, Phe, and Trp) that were not present in the template were placed in the most probable non-clashing rotamer positions. Other side chains that were not present in the template were not modeled since they have uncertain localization. The input files for this target are given at https://cluspro.org/examples/saxs.zip. The zip file contains two PDB files, rec.pdb, which provides coordinates of the model for the receptor, the salmon goose-type lysozyme, and lig.pdb, the e. coli PliG lysozyme inhibitor (PDB ID 4DY3). The zip file also includes the relevant SAXS profile file, SalGEc_new_autoRg.dat. Details of the SAXS profile were described in our report on extending ClusPro to accounting for SAXS data⁵⁸.

ClusPro provided a near-native model even without the use of the SAXS restraints, but it was the center of the 6^th largest cluster. Essentially the same near-native model was obtained when accounting for the SAXS data (see cyan cartoon in Fig. 11), but now as the center of the 3^rd largest cluster. The conformation of the ligand, the e. coli PliG lysozyme inhibitor, is shown as magenta cartoon for reference (Fig. 11). The Others Mode can be combined with the SAXS Profile Mode, and this was used for the calculation.

Open in a separate window

FIGURE 11

Docking the E. coli PliG lysozyme inhibitor to the salmon goose-type lysozyme using SAXS data as restraints. A near-native model (shown cyan cartoon) was obtained as the center of the 3^rd largest cluster. The X-ray conformation of the inhibitor is shown as magenta cartoon.

Constructing a dimer in Multimer Mode

As an example of multimer modeling, we constructed the dimer of the sugar aminotransferase AtmS13 from Actinomadura melliaura. Since no monomeric structure of the protein was available, we have first constructed a homology model. The sequence of AtmS13 from Actinomadura melliaura (Uniprot entry {"type":"entrez-protein","attrs":{"text":"Q0H2X1","term_id":"123029125","term_text":"Q0H2X1"}}Q0H2X1) was downloaded from the Uniprot website at http://www.uniprot.org in FASTA format. To build a homology model using the HHpred server http://toolkit.tuebingen.mpg.de/hhpred, the amino acid sequence was pasted into the box provided on the HHpred home page, and submitted for template selection and alignment using the pdb70 database for hidden Markov model (HMM) construction. All other parameters were left at their default values. HHPred returned a large number of potential template sequences aligned to the AtmS13 sequence, with the sugar aminotransferase CalS13 from Micromonospora echinospora (PDB ID 4ZAS) at the top of the list with the highest similarity score (63% sequence identity). The next step was clicking on the label Create model, which provided the options of either creating a model from the manually selected templates or automatically selecting the best templates. We have used the first option. Notice that building a homology model by HHpred requires the user to have a license to the MODELLER program¹¹⁹. The license is freely available for academic users and easily obtainable at http://salilab.org/modeller/registration.shtml. Entering the coordinates of the homology model constructed as the receptor, we used the Multimer Docking option of ClusPro with the number of subunits set to 2 and left all the other options at their default values. Once the docking was completed, the resulting models were compared to the available crystal structure of AtmS13 homodimer (PDB ID 4XAU, chains A and B). In agreement with our prior observation for dimers, the balanced coefficient provided the best ranking of near-native models. In particular, a near native docking model with IRMSD of 2.62 Å was ranked second (Fig. 12).

Open in a separate window

FIGURE 12

Constructing the dimer of the sugar aminotransferase AtmS13 from Actinomadura melliaura by homology modeling and multimer docking. Chains A and B of the target dimer are shown in green and light green, respectively. Chain B, predicted from the homology model of chain A with the IRMSD of 2.62 Å, is shown as orange cartoon.

This investigation was supported by grants R35 GM118078 and R01 GM061867 from the National Institute of General Medical Sciences, and grants DBI 1147082, DBI 1458509, and AF 1527292 from the National Science Foundation.