12.2.1 HLA Proteins

In humans, the major histocompatibility complex (MHC) is encoded by the HLA, which is responsible for determining cell similarity and signaling to T-cells for immune system activation against pathogens. HLA proteins can be categorized into two classes: Class I and Class II. To date, the IMGT/HLA database (accessible at www.ebi.ac.uk/ipd/imgt/hla/) has over 15,000 HLA allele sequences with 11,000 Class I alleles and about 4,000 Class II alleles [7]. In general, Class I and Class II HLA perform the same function by presenting peptides to T-cells for immune system activation [8]. However, beyond their functions, the two classes of proteins become quite distinct. One major distinction is that Class I HLA presents a peptide to CD8⁺ T-cells; whereas Class II HLA presents a peptide to CD4⁺ T-cells [8]. Another major difference between the two classes occurs due to the source of peptide. The normal peptide source for Class I HLA occurs from an endogenous antigen, while the normal peptide source for Class II HLA comes from an exogenous antigen [9]. However, under particular conditions these peptide sources can be reversed. In the case of cross-presentation, Class I HLA will present an exogenous peptide [10] and when autophagy protein degradation occurs, Class II HLA can present an endogenous peptide [11]. There are also notable differences in the structures of Class I and Class II HLA which are presented in Figure 12.1.

Even though T-cell signaling occurs at the cell surface, the formation of Class I HLA occurs inside the cell in the endoplasmic reticulum (ER) [8]. The Class I heterodimer is assembled in the ER consisting of a 270 amino acid residue peptide binding domain and a 90–100 amino acid residue β₂-microglobulin light chain. Interestingly, it is not stable until after an 8–10 amino acid residue endogenous peptide is translocated into the ER to stabilize the HLA protein [8, 12]. The C- and N-termini of the peptide form H-bonds that anchor the peptide in the Class I pocket [8, 12]. After the HLA-peptide complex is formed, it is released from the ER and presented at the cell surface [9]. A more detailed exploration of the pathway can be found in the excellent review of MHC Class I and II pathways by Neefjes et al. [8].

The peptide binding domain of Class I HLA can be subdivided into three 90-residue α-domains. Two of these α-domains (α₁ and α₂) are directly involved in stabilizing the co-binding peptide between two α-helical walls connected by a β-pleated floor [12]. Once Class I HLA is presented at the cell surface, the third α-domain (α₃) and the β₂-microglobulin chain provide anchoring sites to the cell surface [12]. Interestingly, even though there are over 11,000 Class I alleles, the overall protein structure is highly conserved [7]. A schematic of the structure for Class I HLA variants is provided in Figure 12.1.

Class II HLAs share a lot of similarities with Class I HLAs such as the polymorphism, three-dimensional structures, genetic location, and function [8, 12]. However, there are several notable differences between the two classes beginning with the fact that there are significantly fewer known sequences of Class II HLA compared to Class I (4,000 known Class II alleles) [7]. Two other major differences concern the tissue distribution and the antigenic peptides between the Classes. Class I HLA is ubiquitously expressed in most cells, whereas Class II HLA is expressed in specific cells (APCs, DCs, macrophages, and B-cells) [8]. As discussed earlier, Class I antigenic peptides are 8–10 residues in length, while Class II antigenic peptides can be 12–24 residues in length [8, 12].

Antigenic peptide presentation in Class II HLA begins in the ER by the formation of an α- and a β-heterodimer that are then stabilized by an invariant chain (Ii) [8]. Then the Ii-Class II HLA is transported from the ER to the MHC Class II Compartment (MIIC) where the Ii chain is digested to form a residual Class II-associated Ii peptide (CLIP) [8]. This CLIP peptide is typically 83–107 residues in length [13]. Finally, the CLIP peptide is replaced by a 12–24 mer antigenic peptide degraded from proteins in the endocytic pathway [14], before the Class II HLA protein is translocated from the MIIC to the cell surface. Again, a more thorough analysis can be found in the very detailed review by Neefjes et al. [8].

Structurally, Class II HLA consists of an α and a β heterodimer, with each section containing 90–100 residues [12]. Each of these regions can then be subdivided into two subdomains that are membrane-distal (α₁ and β₁) and membrane-proximal (α₂ and β₂). The membrane-distal region is where peptide binding occurs, while the membrane-proximal region holds Class II HLA to the cell surface. Comparing the structure of Class II with Class I, it becomes apparent that the α₁ and β₂ subdomains of Class II correspond to the peptide binding pocket α₁ and α₂ subdomains of Class I. As such, the membrane-distal region of Class II is composed of two alpha helical walls with β-pleated strands forming the peptide binding pocket floor. However, instead of antigenic peptides non-covalently binding in anchor fashion (as in Class I), Class II peptides undergo H-bonds between the pocket and center of the peptide which allows the C- and N-termini to hang past the edges of the peptide pocket (allowing for longer peptides to bind Class II HLA) [8, 12]. Continuing the comparison between Class I and Class II, it is known that the α₃ and β₂-m subdomains of Class I correspond to the membrane-proximal subdomains of Class II (α₂ and β₂). An overview of the structure of Class I and II HLA alleles is provided in Figure 12.1.

12.2.2 ADR–HLA Associations

Idiosyncratic ADRs can occur anytime a drug binds with an HLA variant and activates either CD4⁺ or CD8⁺ T-cells. Upon T-cell activation a patient can suffer from a multitude of ADRs such as drug-induced liver injury (DILI), Stevens-Johnson syndrome/toxic epidermal necrolysis (SJS/TEN), hypersensitivity, drug reaction with eosinophilia and systemic symptoms (DRESS), malignant pleural effusion (MPE), or agranulocytosis/granulocytosis to name the most well-known examples. Table 12.1 provides examples of drugs and their reported OR for a particular HLA-variant; interestingly, some drugs, such as carbamazepine, are listed as being active at multiple variants [18–20]. The OR is a population-based statistical measure that assesses the likelihood of an exposure (drug-HLA variant interaction) to cause a particular outcome (ADR event occurs) [31]. This measurement is obtained using the following equation:

12.1

^a Putative data, unpublished.

There are three possible outcomes in Equation 12.1; first, OR = 1 in which case it can be concluded that exposure does not affect the odds of the outcome [31]. Second, if OR > 1, then the exposure is associated with a higher odds of the outcome, and finally, if OR < 1, then the exposure is associated with a lower odds of the outcome [31]. The recently built HLADR database compiles all the known drug-HLA associations with their respective ORs [32].

Let us take a real world example using ORs to determine a potential drug-HLA relationship. Clindamycin is an antibiotic that causes either delayed maculopapular exantherma or SJS/TEN in patients expressing the HLA-B*51:01 variant [21]. However, Yang et al. found that there were three other variants that had OR greater than one, indicating that the HLA variants B*15:27 (OR: 55), B*51:01 (OR: 24), C*14:02 (OR: 9), and C*15:02 (OR:19) could all be binding clindamycin to trigger an ADR event [21]. These results raise several questions: Do all of these alleles result in ADR events? Are drug-HLA interactions more favorable for a particular variant compared to the others? Will drug-HLA binding be strong enough to induce the ADR? Clearly the limitations of using OR to diagnose drug-mediated ADR events have been reached and alternative techniques are needed to analyze and elucidate drug-HLA complexes. In this particular case, Yang et al. attempted to perform a molecular docking study using a HLA-B*51:01 crystal (PDB:1E27) and homology models for HLA-B*15:27 (PDB Template:1XR8), HLA-C*14:02, and -C*15:02 (PDB Template: 1EFX) to determine which drug-HLA combination was most likely [21]. From their docking results, Yang et al. concluded that HLA-B*51:01 was most likely responsible for inducing ADR (GLIDE docking score DS_XP = −10.1 kcal/mol) [21]. Unfortunately, the docking results reported by Yang et al. indicated that docking at all four variants was actually possible and predicted as favorable (all docking scores < −7 kcal/mol) [21]. Additionally, the authors did not mention how the co-binding peptide was analyzed or taken into account in their docking procedure despite the fact that the peptide has a significant role in HLA stability [8, 12]. The role of the co-binding peptide with HLA-B*57:01 and abacavir will be discussed in greater detail in Section 12.3.2.

Another well-studied example is the drug carbamazepine. Carbamazepine is used to treat epilepsy and bipolar disorder, but in patients who have one of the following HLA-variants, SJS/TEN may occur: HLA-B*15:02 (OR: 895), HLA-B*15:11 (OR: 31), or HLA-A*31:01 (OR:33) [18–20]. In 2012, Wei et al. confirmed carbamazepine's affinity for HLA-B*15:02 through the use of T-cell proliferation and cytotoxic assays in addition to a HLA-B*15:02 peptide-specific assay [33]. Interestingly, the assay results also indicated that HLA-B*15:08, -B*15:11, and B*15:21 would have some binding affinity with carbamazepine [33]. Furthermore, based on the basis of their site-directed mutagenesis work and homology model (using HLA-B*15:01 as a template, PDB: 1XP8), Wei et al. were able to identify ARG62 and GLU63 as two critical amino acids for carbamazepine binding [33].

Using this previous work, two interesting studies have recently emerged. First, Schotland et al. used molecular analysis of side effects (MASE) to data-mine the FAERS database for drug-HLA associations for SJS/TEN and found that patients taking acetominophen (APAP) and carbamazepine showed a reduced chance of SJS/TEN [34]. Using the homology model developed by Wei et al. [33], Schotland et al. tested the binding affinities of carbamazepine and APAP and found that APAP had a higher binding affinity with HLA-B*15:02 [34]. The second study used the same homology model, but now involved performing molecular dynamic (MD) simulations to further understand the T-cell signaling pathway [35]. Zhou et al. [35] have proposed an interesting signaling pathway that is highlighted in Section 12.2.3.

Thus far, only interactions between Class I HLA and drugs have been discussed. It is important to remember that drug-HLA interactions are also observed with Class II HLA. Clozapine is used for the treatment of schizophrenia, but in certain cases may result in agranulocytosis or granuloctopenia. In 2007, Dettling et al. reported significant OR values for clozapine with HLA-B*57 (OR:22) and HLA-DRB5*02:01 (OR: 22) [36] Then in 2014, Goldstein et al. reported a possible connection for clozapine binding with HLA-B*57:01 or HLA-DQB1*05:02 variants [22]. Most recently, Legge et al. confirmed clozapine's HLA-dependence with HLA-DQB1*05:02 with a reported OR of 15 [23]. Remarkably, these reports indicate that clozapine may bind with either Class I or Class II HLA [22, 23, 36].

The use of ORs provides evidence for a correlation between an ADR event and specific HLA variant(s), but it does not imply a significant binding affinity between a particular drug and HLA-variant. In some instances, (as with carbamazepine) binding assays are available to confirm drug-HLA associations, but even in vitro assays are limited in their ability to elucidate information about specific drug and HLA binding pocket interactions. The best method to visualize HLA-drug binding interactions is through the use of NMR or X-ray crystallography.

The association between abacavir and HLA-B*57:01 is the best understood one, notably owing to the availability of three X-ray crystals identified in 2012 by Illing et al. (PDB: 3VRI and 3VRJ) [37] and Ostrov et al. (PDB: 3UPR) [38]. These three crystals not only confirmed the HLA-drug complex previously reported has OR values (greater than 900) [15, 39] but also provided accurate starting structures for molecular docking and dynamics [40–42]. Recently, Metushi et al. conducted a virtual screening of the ZINC database to identify potential HLA-B*57:01-active compounds with a high similarity to abacavir [41]. After performing a 2D/3D similarity search, they next used a combination of pharmacophore and molecular docking (using PDB: 3UPR from their previous work) to identify several chemicals for an HLA-binding assay and a peptide-elution assay [38, 41]. Using this protocol, Metushi et al. identified acyclovir as a potential HLA-B*57:01-active compound and subjected it to a CD8⁺ T-cell response assay [41, 43]. Using the T-cell response assay, they were able to confirm that acyclovir does not mediate a T-cell response (as expected, this is due to historical lack of ADRs associated with acyclovir) [44–46].

Interestingly, in each of the HLA-drug interaction examples discussed above, molecular docking was used as the main cheminformatics tool to explore the binding modes of drugs in the HLA pocket [21, 22, 33–35, 40–42]. Best practices and limitations when performing molecular docking are discussed in Section 12.3.1. This is followed by an in-depth analysis of our own work in developing a virtual screening protocol at the HLA-B*57:01 variant using molecular docking, in Section 12.3.2 [42].

12.2.3 HLA-Drug-Peptide Proposed T-Cell Signaling Mechanisms

In order for a drug to cause an adverse event, it must signal to T-cells (CD8⁺ for Class I and CD4⁺ for Class II) after forming an HLA-drug-peptide complex. Currently, there are three proposed mechanisms for HLA-T-cell signaling: altered repertoire, pharmacological interaction (p. i.) complex, and the hapten complex [4, 47]. An adaptation of these signaling modes from Illing et al. [48] is provided in Figure 12.2. The altered repertoire signaling mechanism occurs when a drug non-covalently binds in the peptide pocket beneath the peptide as shown in Figure 12.2a. A p. i. complex is based on non-covalent interactions, but forms with the peptide surface outside of the binding pocket as shown in Figure 12.2b. The hapten complex is the only occurrence where the drug is covalently bonded to the peptide and can occur either inside or outside the binding pocket as shown in Figure 12.2c. Interestingly, the altered repertoire and hapten complex can both form either at the cell surface or in the ER during HLA formation, whereas the p. i. complex forms solely on the cell surface [4, 47, 48].

These three signaling mechanisms have been suggested (and/or observed) for abacavir, carbamazepine, and flucloxacillin [35, 37, 38, 48, 49]. To date, the only structurally verified signaling mode is for HLA-B*57:01-abacavir-peptide in an altered repertoire complex (PDB: 3VRI, 3VRJ, 3UPR) [37, 38]. Carbamazepine is believed to form a p. i. complex with HLA-B*15:02-peptide-carbamazepine [48], but recent MD simulations by Zhou et al. [35] suggests that carbamazepine may non-covalently bind to a CD8⁺ cell before non-covalently binding to the peptide surface as shown in Figure 12.3. Without an NMR or crystal structure, it will be difficult to confirm this supposed mechanism. Flucloxacillin is believed to form a hapten complex due to the nature of β-lactams and their ability to form covalent bonds with lysine residues [48–50]. Without NMR structure or crystals of HLA-drug complexes, the use of molecular docking and molecular dynamics offers an opportunity to identify drugs that favor altered repertoire or p. i. complexes with HLA. Unfortunately, molecular docking and dynamics will be unable to identify hapten complexes owing to their use of classical force fields; a quantum mechanical technique would be needed to identify these types of complexes due to the formation of a covalent bond.

12.3 Structure-Based Molecular Docking to Study HLA-Mediated ADRs

Owing to the serious health implications of ADRs, the ability to predict whether a drug can bind with HLA and induce an immune system response is extremely important. However, solely relying on experimental assays to identify all possible drug-HLA relationships is perhaps not the best strategy. Currently, there are over 15,000 HLA-variants listed in the IMGT/HLA database [7]. If an assay were developed for every single variant and an experimental screening of DrugBank were performed (currently experiments performed. The screening of larger libraries such as ZINC database (over 15 million compounds) would simply become infeasible. That number does not even account for experimental errors and replicates in the different assays. Clearly, brute-force experimental screening for HLA-drug interactions is hardly possible if many variants and a large library of compounds are considered.

Meanwhile, in silico virtual screening offers a comparatively quick alternative for identifying potentially HLA–active drug compounds. The use of similarity searching, pharmacophore, quantitative structure-activity relationships (QSAR), and molecular docking could all provide valuable insight into the identification of HLA-active compounds. However, owing to very limited structural information about HLA-drug binding affinities (K_i), statistical techniques such as QSAR modeling cannot be conducted at this time. In contrast, molecular docking represents the most appropriate cheminformatics technique to perform a virtual screening tool for HLA owing to the availability of crystal structures (e.g., three structures for HLA-B*57:01-abacavir crystal structures 3VRI, 3VRJ, 3UPR) and the possibility to build additional homology models for the other variants.

12.3.1 Structure-Based Docking

Molecular docking consists of virtually inserting a ligand into the binding pocket of a target (e.g., protein) and then scoring the actual molecular interactions between the ligand and protein according to the binding mode of the ligand and a scoring function. Docking scores can be used to rank a set of ligands on the basis of their predicted binding affinities toward the target, and thus help identifying and prioritizing those ligand(s) with the most interesting interactions with the target. A classical docking procedure can be performed in three ways: rigid ligand-rigid protein, flexible ligand-rigid protein, and flexible ligand-flexible protein. Rigid-rigid docking is considered to be the least accurate and is only used for screening extremely large databases as a first and rapid initial step. Flexible-rigid docking is the most commonly used technique and is generally considered to produce reliable binding modes and binding affinities between ligands and proteins. In theory, the best approach would be to use fully flexible-flexible docking as this would account for induced conformational changes in proteins upon ligand binding, but this technique is still rarely usedowing to the large computational cost. The rapid development of ensemble docking (i.e., docking flexible ligands in many conformations of the protein) and GPU-accelerated molecular dynamics will accelerate the transition to flexible-flexible docking.

When performing flexible-rigid molecular docking, the starting structure for the protein comes from an experimentally validated 3D structure (X-ray or NMR). Using this known structure, two approaches can be used: (i) docking is performed directly on that structure or (ii) docking is performed by developing a homology model built from this initial 3D structure. In the case of homology modeling, the residue sequence for a protein is known (for example HLA-B*15:02) but the 3D structure is unknown. However, there may be an NMR or crystal structure available of a closely related allele (such as HLA-B*15:01). Using the 3D structure of this closely related variant as a template, a theoretical homology structure can be built. This technique was used in the study of HLA-B*15:02 with carbamazepine [33–35]. Homology modeling is usually followed by several cycles of full-atom molecular dynamics to optimize the overall geometry of the protein before conducting molecular docking.

In some instances, the actual binding mode of a drug with a given protein is directly known from NMR or X-ray crystallography. If this is the case, as with HLA-B*57:01 and abacavir, then molecular docking may be performed directly on that structure after in-depth preprocessing steps have been completed. First, unessential counter ions and waters are removed (the ones in the binding site can be preserved if necessary) and explicit hydrogens are added. Then, missing side chains are identified and corrected. Protonation states and H-bonding are optimized. The overall structure is minimized [51–55]. In our group, we are notably relying on the Protein Preparation Wizard available in the Schrödinger Suite (www.schrodinger.com). Protein Preparation Wizard affords the easy removal of waters, salts, and H-bond optimization [51], while utilizing PRIME to generate missing side chains [52, 53] and protonation states are generated using EPIK at physiological pH of 7 [54, 55]. Finally, protein minimization is performed using the OPLS3 force field [56–60].

After the protein structure has properly been curated, it is important to perform self-docking (sometimes referred to as “native docking”) of the parent ligand from the original crystal. Self-docking is a critical step to verify whether the docking software can accurately retrieve the native binding mode of the ligand. After successful completion of self-docking, the curated and prepared compounds of interest can be docked. Chemical curation is an important step to ensure the correctness of chemical structures to be docked [61, 62]. The curated set of ligands should then be optimized (either with classical- or quantum mechanics) in addition to generating active tautomerized states at physiological pH. Our group performs ligand preparation using LigPrep available in the Schrödinger suite which generates tautomeric states with EPIK [54, 55] and minimizes conformers using the OPLS3 force field [56–60].

With the properly curated and preprocessed ligand database, molecular docking can be performed. There are numerous docking software available: Autodock Vina [63], GOLD [64], FLEXX [65], and GLIDE [66–69]. Our group relies on Schrodinger's GLIDE for docking owing to the availability of multiple scoring functions, HVSP, SP, XP, and SP-PEP [66–69], as well as the overall prediction performances of the software [67]. Briefly, GLIDE operates by docking in three main steps: (i) performs initial screening of ligand conformers in the active site of the receptor, (ii) minimizes the energy according to a force field (e.g., OPLS3), and (iii) scores the molecular interactions to assess the estimated free energy of binding of the ligand [66]. The GLIDE DS estimates the free energy of binding of a ligand in the binding pocket and is formally the sum of the Glide score and the EPIK state penalty occurred from generating ligand tautomers [51, 66]. The eModel (eM) score is the sum of the Glide score, ligand-receptor interaction energy, and the ligand strain energy [66]. When studying HLA, we used two empirical thresholds to identify HLA-active compounds. A drug was considered HLA binder or “active” when it was found to afford a DS below −7 kcal/mol and an eM score below −50 kcal/mol.

Recently, we developed a virtual screening molecular docking procedure for the HLA-B*57:01 variant using three crystals of HLA-bound abacavir (PDB: 3VRI, 3VRJ, 3UPR) [42]. Interestingly, each structure had a different co-binding peptide which allowed us to explore the importance of that peptide in HLA-drug binding. In this work, we did not use ensemble docking as we wanted to develop a fast virtual screening procedure targeting the HLA-B*57:01 variant and the fact that abacavir is found in an altered repertoire complex led us to hypothesize that there is little variation in the protein structure.

12.3.2 Case Study: Abacavir with B*57:01

Recently, we developed and tested a virtual screening protocol for identifying HLA-B*57:01-active drugs using flexible-ligand and rigid-protein molecular docking [42]. The HLA-B*57:01 variant was selected owing to the availability of three X-ray crystals consisting of HLA-B*57:01-abacavir and a co-binding peptide in an altered repertoire complex (PDB: 3VRI, 3VRJ, 3UPR) [37, 38]. Owing to the availability of these three crystals with a confirmed drug-binding mechanism and known HLA-B*57:01-drug associations resulting in ADR events (AHS and DILI), we decided to use these structures as a starting point for developing a virtual screening protocol. A schematic of our protocol is provided in Figure 12.4.

As described in Section 12.3.1, we started by curating the three protein structures using the Protein Preparation Wizard available in the Schrodinger Suite [51–55]. Once the three proteins were preprocessed, we conducted a protein alignment to identify the differences between the three binding pockets, the three units of abacavir, and the three co-binding peptides. The protein pocket was found to be in the same conformation for all three crystals with measured root mean square Deviations (RMSDs) less than 0.6 Å and overlay similarities greater than 75% [42]. Furthermore, we found that aside from ring puckering in abacavir's cyclopropyl group or orientation of the hydroxyl functional group, the structure of abacavir in the pocket was universal among all three crystals [42]. The largest variation between crystals was in the co-binding peptide as each peptide had a unique amino acid sequence for 3VRI (P1: RVAQLEQVYI), 3VRJ (P2: LTTKLTNTNI), and 3UPR (P3: HSITYLLPV) [42]. Interestingly, when a peptide backbone alignment was performed, it was found that the backbone of each peptide was in a similar orientation with a measured RMSD of less than 1.8 Å [42]. On the basis of alignment results, we moved forward with the self-docking abacavir.

In a molecular docking experiment with a known binding mode between a protein and a drug, self-docking is an essential step in validating the ability of molecular docking to retrieve that particular binding mode with very high accuracy. In the case of HLA, self-docking is even more important owing to the presence of a co-binding peptide. When a flexible-ligand/rigid-receptor protocol is used, GLIDE fixes the entire protein in place. Importantly, in the case of HLA complexes, this includes the co-binding peptide. As such, we conducted docking in the absence and presence of the co-binding peptide using the SP and XP scoring functions of GLIDE [66–68] in order to understand the specific role of the peptide in docking [42]. These results are provided in Figure 12.5. The SP and XP scoring functions were both able to reliably reproduce the native binding pose of abacavir in the absence and presence of peptide P1 with measured RMSD of 1.21-1.28 Å provided in Figure 12.5a [42]. Interestingly, the differences between structures occurred as a result of ring strain in the cyclopropyl moiety and the hydroxyl functional group.

The observed protein-drug interactions between HLA-B*57:01 and abacavir were identical between the native crystal and the docking-generated binding modes except with respect to the hydroxyl group [42]. Strangely, the hydroxyl group was found to participate in the H-bond with three different amino acid residues as shown in Figure 12.5b [42]. Native abacavir from 3VRI has H-bonding between the hydroxyl group and the TYR74 residue, while SP docking without peptide found the hydroxyl group to H-bond with the TYR99 residue, and SP docking with peptide found that the hydroxyl group H-bonds with carbonyl backbone of ALA3 from P1 (Figure 12.5b) [42]. Owing to the proximity of TYR74, TYR99, and ALA3, it is possible that in a dynamic system, the hydroxyl group actually rotates and H-bonds with all three residues. In order to validate this hypothesis, molecular dynamic simulations still need to be performed. We also explored the role of the co-binding peptide in stabilizing the HLA-B*57:01-drug complex. It was observed that all three peptides (P1, P2, and P3) lowered the DS by 1–2 kcal/mol (Figure 12.5c) and that the eM score was lowered by 12-17 kcal/mol (Figure 12.5d) [42].

Next, a small test set of 14 drugs and abacavir were docked at 3VRI, 3VRJ, and 3UPR in the absence and presence of peptides P1, P2, and P3, respectively [42]. The test set consisted of abacavir, allopurinol, atorvastatin, carbamazepine, ciprofloxacin, clozapine, fenofibrate, flucloxacillin, methyldopa, minocycline, pazopanib, sertraline, simvastatin, and ticlopidine [42]. Associated ADRs and indications of these drugs can be found in Table 12.1, while the structures are provided in Figure 12.6. The test set of compounds were docked using both SP and XP scoring functions using protein structures from 3VRI, 3VRJ, and 3UPR in both the absence and presence of peptides P1, P2, and P3 [42]. Heat maps of the DS and eM scores are provided in Figures 12.7 and 12.8, respectively, where green boxes indicate the lowest measured DS or eM score (most favorable), red boxes indicate the highest measured DS or eM score (least favorable), and white boxes indicate that no favorable binding pose between drug and HLA-B*57:01 were obtained.

Interestingly, the role of peptide becomes quite apparent when viewing either Figure 12.7 or 12.8: there are three drugs that have favorable DS scores (less than or equal to −7 kcal/mol) and favorable eM scores (less than or equal to −50 kcal/mol) in the absence of P1, P2, or P3, but do not pass the thresholds when peptide is present [42]. Interestingly, fenofibrate in the presence of P1 passes both thresholds, but fails both thresholds when docked in the presence of P2 or P3. For the test set of drugs DS values range from −1.8 (flucloxacillin, XP + P2) to −10.3 (fenofibrate, XP + P1); however, there is significantly more variance in the eM scores. The lowest measured eM score is −91.6 (pazopanib, XP – P1), while the highest eM score is positive 112.2 (simvastatin, XP + P1) [42]. Intriguingly, there are four drugs that exhibit positive eM scores when the XP scoring function is employed in the presence of either P1, P2, or P3 (clozapine, fenofibrate, flucloxacillin, and simvastatin) [42]. When these positive eM scores are observed, it is expected that the corresponding DS would also be high; however, in some instances this is not the case and in fact the corresponding DS for fenofibrate, flucloxacillin, and simvastatin are rather negative. For example, when simvastatin was docked using XP + P1, the resulting DS was −10.3 kcal/mol, while the eM score was 112.2 kcal/mol. Simvastatin's DS indicates that the drug forms extremely favorable interactions with the binding pocket, but the eM scores indicate that this pose will be unfavorable. It seems unlikely that both situations can be true to such an extreme extent. These ambiguous results indicate a clear limitation of applying molecular docking to a tripartite complex such as HLA-drug-peptide, where the dynamic interactions between drug and peptide are extremely important. Therefore, molecular dynamics simulations will be needed to further evaluate, analyze, and understand these interactions.

Overall, our model was able to reproduce the binding mode of crystallized abacavir and successfully docked a test set of drugs with the HLA-B*57:01 variant. A full summary of these results are provided in Figure 12.9, where green cells represent drugs that passed both threshold criteria using SP and XP functions, yellow cells represent drugs that only passed using the XP function, orange cells represent drugs that only passed using the SP function, and red cells represent drugs that failed both SP and XP scoring functions. Furthermore, our model was able to successfully distinguish active pazopanib as an HLA-B*57:01-active compound. Interestingly, our model failed to identify flucloxacillin as an HLA-B*57:01-active compound. This result may occur because flucloxacillin is suspected to form a hapten complex with lysine containing peptides [49, 50] and molecular mechanics cannot create covalent bonds. Although there are some clear limitations with applying molecular docking as a virtual screening procedure with HLA variants, we feel confident that our model will be able to identify previously unknown HLA-B*57:01-active compounds when we conduct a screening of the DrugBank and Tox21 databases.