GPCR modeling tips and tricks

After speaking with John and giving a general overview of important considerations in modeling GPCRs, he highlighted the number of "branch points" at which human decisions need to be made before a final set of models is ready for docking and selection. I wrote up some documentation summarizing key practices and insights from experience that help in generating high-quality docking-ready models.

Template selection

For orphan GPCRs, the first stop in choosing a template structure is the Gloriam group's GPCR database. In the current website design, choose "Template Selection" from the "Receptors" menu and enter the name of your favorite receptor, and you will get a table listing all available template structures ranked by sequence similarity in the TM region calculated from the GPCRdb pre-computed, optimized sequence alignment (see below). Templates over 20% sequence identity in the TM region are preferable. It's generally better to separately build model ensembles from each template you like, rather than building multi-template models. (These often produce poorly packed TM helices.) There is surprising functional conservation of the molecules predicted from a model, so it may be preferable to try an active-state template if one is available.

Alignment

Particularly for orphans of low sequence identity to the available templates, it may be useful to build models based on more than one sequence alignment. Several tools exist for generating candidate alignments:

• GPCRdb from the Gloriam lab: GPCRdb alignments can be obtained between any two receptors by clicking Receptors > Structure-based alignments and entering your target GPCR and your favorite template. (This can also be done via the API.) These have been manually reviewed and are generally high quality in the TM helices, but require additional review (see below) - I've observed errors frequently in ECL2 and in the register of TM5. A set of precomputed models based on these alignments, using a library of precedented sidechain rotamers, is available for download and these make useful "baseline" comparators when performing docking later in the process.
• GPCR-HGmod from the Zhang lab: uses the GPCR-I-TASSER method; has a database of precomputed alignments, and precomputed models based on these alignments, for most GPCRs, and a webserver interface for modeling those not in the database.
• HHpred: an integrated system for searching for templates, producing alignments, and producing models (the webserver runs MODELLER in the back end). Not GPCR specific.
• EVfold: a tool for predicting structures based on evolutionary constraints, runs JackHMMer by default. Not GPCR specific.
• PROMALS3D: multiple sequence and structure alignment. Expects the user to provide candidate templates and performs best with multiple template structures. This is not GPCR specific and will often differ from GPCRdb in the predicted register of TM5.

Review of key sequence features in the candidate alignments is advisable before building model ensembles, particularly if the alignment was generated from one of the non-GPCR-specific tools. Common trouble spots include:

• Disulfide bonds: if you expect that the disulfide bond(s) in the extracellular loops should be conserved, double-check that the cysteines are aligned correctly, particularly in ECL2 where this is a common error in GPCRdb alignments.
• Conserved GPCR sequence motifs: some orphan GPCRs lack common features (DRY motif, sodium site, NPXXY motif) and consequently are slightly more difficult to align correctly, but if your target sequence has these features, a candidate alignment should be viewed with skepticism if it locates these out of register.
• The position and register of TM5 is a common source of error; for remote orphans, different tools may identify completely different regions of sequence as corresponding to the TM5 helix, and even when they agree, i+4 variations (that is, alignments out of register to one another by one helical turn) are not uncommon. In my experience it is best to build test models for more than one alignment if TM5's predicted position varies; often one alignment will produce noticeably worse models despite being similar by alignment-scoring methods. (For example, this was a recurring issue with GPR162, a remote orphan with a weak hit reported by the Roth lab).

Preliminary models

It is often convenient to build a few preliminary models using your preferred template using the Chimera GUI interface to give them a quick look before committing to a larger ensemble. There is a tutorial here.

1. Fetch your favorite template by PDB ID (File > Fetch by ID) or load your .pdb file locally.
2. Open your target-template alignment chosen by any of the methods above: Tools > Sequence > MultiAlign Viewer.
3. Associate the template structure with its corresponding sequence: in the MultiAlign Viewer, click Structure > Associations and click through the prompts. Red boxes will highlight sequence variations/mutations/missing sequence.
4. Build your test models: Structure > Modeller (homology) and follow the prompts
   1. You have the option to use a webserver (you will need to enter the Modeller key) or the local install (currently in /nfs/soft/modeller/current/bin/mod9.15 on the new cluster)
   2. The tool will write out copies of your template PDB, the modeller input file (ModellerModelling.py), and various log and config files, load the models back into the Chimera workspace, and show a new window with a score summary.
5. Evaluating models: the tool shows the GA341 and zDOPE (normalized DOPE) scores for each model. GA341 is insensitive and should be high for good models - above 0.8, ideally close to 1. zDOPE is more sensitive but not parameterized for membrane proteins.
6. You will probably not get good models from this preliminary look - depending on whether you pre-processed your alignments, you may see long loops, untrimmed termini, etc.
7. To improve your test models, trim the termini and loops (you can add a short linker of 4 or 5 alanines to avoid creating gaps in the structure) and add disulfide constraints using the disulfide residue patch.
8. You can run a sample model through MolProbity for a separate check on its quality.

Model ensemble options

It is desirable to sample broadly in conformational space when building models. Modeller can generate sets of models from a single alignment and template input, but these are usually relatively similar in conformation, varying mostly in loop conformation and sidechain rotamers.

• 3kENM: Elastic network model sampling around a given conformation. (Credit to Joel Karpiak and Ryan Coleman for the idea of using this for GPCR orphans). The tool is described at Enm explorer. (The .crd file can be split into PDBs in various ways; a VMD TCL script is available in my scripts directory at /mnt/nfs/home/kstafford/scripts/breakpdbs.tcl (for use within VMD) or /mnt/nfs/home/kstafford/scripts/breakpdbs_all.tcl (for use on the command line). This can be applied to the template structure and then the resulting backbone conformations used as template input to Modeller.
• Rosetta: Best used to sample rotamer packing for models with reasonable backbone conformations.
• PLOP: Best used to optimize loops, especially after a round of docking.
• AllosMod: Available as a webserver here, takes a couple of hours to sample around a single conformation (ie, within a single well). I tried using this with multiple conformations (wells) to sample models for which the best templates were in the inactive state with the aim of constructing active-state models, but this was mostly a lot of time for small effects in enrichment.

Docking setup

With a diverse model ensemble, it can be difficult to choose a high-throughput way to set them up for docking in order to find candidates with the best enrichment.

• Simplest route: Make sure each model is aligned to its template and use the crystallographic ligand as the xtal-lig. A script to do this is available at /mnt/nfs/home/kstafford/scripts/model_ensemble_setup.py.
• The FTMap webserver is often used in the lab to sample possible positions of small fragments in a protein in order to choose an xtal-lig. This is usually done by visual inspection of the fragment clusters, but once prepared the xtal-lig file can be used for multiple aligned models with similar binding site shapes.
• Thin spheres are useful especially for large or well-solvated sites (e.g. GPCRs predicted to bind peptides or using a peptide-binding receptor as a template).
• After a round of docking that produces good enrichment, it is often possible to use the known ligands' docked positions as the "ligand" in subsequent rounds of docking/screening.

Ligand building/decoys

When working with orphans, we typically had information from the Roth lab about one or a few molecules that bound to a given protein, and a correspondingly large number of compounds in their screening set which could be considered known inactives. Therefore most models were screened against multiple sets of compounds:

• The full set of experimentally tested compounds in the Roth lab screen, from which the actives are designated "ligands" for purposes of enrichment calculations
• Any known molecules from other sources, e.g. CHEMBL, literature/patent searches, information from collaborators, etc., with a preference for compounds listed as potential ligands by IUPHAR (however, use with caution, as many reported ligands were tested by the Roth lab and their activity was not reproduced).
• A set of computationally generated decoys, usually 50 decoys for each ligand, for both the Roth lab hits and the "known" hits from other sources (if any).

In June 2015, I built a set of the compounds from the Roth lab, available in /nfs/work/kstafford/roth_collection_201506. The ligand pipeline has since been updated; see subdirectory zinc15_update. These compounds are also available as the IDG catalog in ZINC.

Enrichment

Enrichment is the primary metric for judging the predictive quality of candidate models. Models with apparently similar properties in e.g. quality scores or RMSD may have quite different enrichments due to relatively small factors such as choice of rotamers within the binding pocket. Often it is necessary to further optimize models that perform well in order to fine-tune rotamers, and sometimes models with apparent high enrichment are on inspection not very good, for example because they suggest implausible poses.

• Enrichment should generally be calculated against both the Roth lab inactives and the generated decoys, since these may offer different challenges (though in practice they're quite similar).
• Enrichment should generally be calculated using just the Roth lab hit list and separately including the "known" ligands - sometimes these are probable false positives and do not activate the receptor in the Roth lab's hands.
• High rank correlation (that is, producing similar compound rankings) between different models may be a useful validation step, as this indicates the models are likely responding to similar features of the molecules.