- Research article
- Open Access
Adverse drug reactions triggered by the common HLA-B*57:01 variant: a molecular docking study
© The Author(s) 2017
- Received: 27 July 2016
- Accepted: 24 February 2017
- Published: 4 March 2017
Human leukocyte antigen (HLA) surface proteins are directly involved in idiosyncratic adverse drug reactions. Herein, we present a structure-based analysis of the common HLA-B*57:01 variant known to be responsible for several HLA-linked adverse effects such as the abacavir hypersensitivity syndrome.
First, we analyzed three X-ray crystal structures involving the HLA-B*57:01 protein variant, the anti-HIV drug abacavir, and different co-binding peptides present in the antigen-binding cleft. We superimposed the three complexes and showed that abacavir had no significant conformational variation whatever the co-binding peptide. Second, we self-docked abacavir in the HLA-B*57:01 antigen binding cleft with and without peptide using Glide. Third, we docked a small test set of 13 drugs with known ADRs and suspected HLA associations.
In the presence of an endogenous co-binding peptide, we found a significant stabilization (~2 kcal/mol) of the docking scores and identified several modified abacavir–peptide interactions indicating that the peptide does play a role in stabilizing the HLA–abacavir complex. Next, our model was used to dock a test set of 13 drugs at HLA-B*57:01 and measured their predicted binding affinities. Drug-specific interactions were observed at the antigen-binding cleft and we were able to discriminate the compounds with known HLA-B*57:01 liability from inactives.
- HLA variant
- Molecular docking
- Virtual screening
Certain patients develop harmful adverse drug reactions (ADRs) after taking a medication . Unfortunately, these undesired reactions to a drug or its metabolite(s) can potentially be serious and even life threatening. According to the FDA Adverse Event Reporting System (FAERS), over one million cases of ADRs were observed in 2014 and this number is steadily increasing as reporting systems become more accessible to physicians . There are two main classifications of ADRs: Predictable ADRs occur due to the pharmacological activity of a drug or its metabolites, whereas idiosyncratic ADRs are primarily observed as an immune system response [3, 4]. The major biological pathway capable of triggering such idiosyncratic ADRs is activated by a drug’s direct binding with human leukocyte antigen (HLA) protein variants [4–6].
Even though the vast majority of HLA variants are rare, there are several well-known cases of HLA-mediated ADRs. For instance, carbamazepine (brand name Tegretol) is used for the treatment of seizures, but patients carrying the HLA-B*15:02 variant [12, 13] have an increased likelihood to suffer from the Steven-Johnson Syndrome (SJS) or toxic epidermal necrolysis (TEN) . Another well-established example is flucloxacillin (Floxapen), a beta-lactam anti-microbial, that can cause typical drug-induced liver injury (DILI) in patients possessing the HLA-B*57:01 variant . Additionally, HLA-B*57:01 is responsible for abacavir hypersensitivity syndrome (AHS) which is a severe, life-threatening ADR occurring in patients prescribed abacavir for treating HIV infection .
Clearly, the ability of a drug to bind with a given HLA-variant plays a significant role in determining whether that compound or its metabolite(s) may potentially trigger ADRs in a given subpopulation carrying that particular HLA variant. Therefore, computational approaches able to forecast such HLA–drug molecular interactions reliably could have serious implications in preventing ADRs and thus potentially contribute to the development of precision medicine. Using in silico techniques to predict drugs’ liabilities has now become a common, reliable enough, and cheap enough screening approach, especially for toxicity and ADRs evaluations [17–30]. Typically, forecasting potential ADR involves analyzing the chemical space with chemical similarity techniques [17–20]. Liu et al. recently used a 2D structural alert—based screening for chemical similarity to forecast ADR , while Vilar et al. employed a 3D pharmacophore-based similarity search in order to predict ADR . Alternatively, quantitative-structure activity relationship (QSAR) models have been developed in order to forecast drug-induced SJS  or DILI . New methods can use a systems chemical biology approach in order to predict drug hepatoxicity through the integration of chemical and biological data [21, 22].
However, there are very few molecular modeling studies in the literature attempting to analyze and predict the molecular interactions between drugs and HLA variants. For instance, Luo et al.  modeled the interactions between various HLA variants and some endogenous peptides using a network analysis approach, but the authors did not examine potential drug binding events. Recently, Paul et al. developed two approaches for predicting HLA-Class I and -Class II epitopes using TepiTool  and the Immune Epitope Database and Analysis Resonance . Another group developed a very useful, online database compiling all of the known HLA–drug interactions resulting in ADRs (HLADR, http://pgx.fudan.edu.cn/hladr/); however, this database is based solely on measured odds ratios (ORs) obtained from existing literature . Recently, Yang et al.  conducted a preliminary molecular docking study on abacavir using AutoDock Vina, but little details were discussed regarding the actual binding mode of the drug. In a proof-of-concept study, we used molecular docking to predict the binding modes of clozapine with several HLA-variants and explore some possible clozapine–HLA interactions . Clozapine (Clozaril) is an efficient antipsychotic that may result in agranulocytosis/granulocytosis when a patient has HLA-DQB1 or specific HLA-B variants [26, 27].
Alfirevic et al.  attempted to establish a HLA-typed DNA archive that could be used to map DILI events between class I and II HLAs using distance trees. Recently, Schotland et al.  attempted to data-mine FAERS reports for SJS/TEN associations with HLA-variants using the Molecular Analysis of Side Effects (MASE) approach. Additionally, Schotland et al. performed a homology docking model of carbamazepine (and several other drugs) at the HLA-B*15:02 variant . This homology model was developed from work previously performed by Wei et al. . Both research groups were able to successfully verify the importance of the ARG62 residue for carbamazepine binding at B*15:02 [29, 30]. Furthermore, using this same homology model Zhou et al. conducted a molecular dynamic simulation exploring the T-cell signaling mechanism of bonded carbamazepine with HLA-B*15:02 . Overall, the literature on the computational modeling of HLA–drug complexes is limited but definitely emerging.
Recently, Metushi et al. conducted a virtual screening of the ZINC database in order to attempt HLA-B*57:01 liable chemicals . The ZINC database contains over 35 million commercially available compounds . Using concatenated FP2 and FP4 structural fingerprints, Metushi et al. conducted a 2D-similarity search of abacavir on 3.5 million compounds from the ZINC database followed by a 3D-similarity search using pharmacophoric features of abacavir . From this initial screening, 54 compounds were identified and selected for molecular docking using the X-ray crystal 3UPR. Next, Metushi et al. identified seven compounds that were tested for HLA-B*57:01 affinity from which acyclovir was identified as a potential candidate . But, when acyclovir was subjected to a CD8+ T-cell response assay , it was determined that acyclovir did not induce a CD8+ T-cell response . This study by Metushi et al.  represents a full in silico to in vitro screening for HLA-B*57:01 liable compounds from ZINC.
Developing molecular docking protocols that effectively identify hits can be a challenging undertaking, especially when it comes to the preparation of complex proteins, such as HLA . Moreover, when molecular docking was conducted in published studies involving HLA proteins, it was not specified if a co-binding peptide was present or absent. As such, we believe that it is of the utmost importance that a thorough analysis of molecular docking targeting the HLA-B*57:01 variant be conducted in order to properly identify the limitations of this molecular modeling technique for forecasting a drug’s likelihood to bind a HLA variant and thus potentially cause a HLA-mediated ADR.
Overall, this study underlines the relevance of employing molecular docking for analyzing and predicting HLA–drug interactions. Due to the large number of HLA-alleles  and their varying frequency of occurrence by ethnicity , the ability to forecast idiosyncratic ADRs is extremely difficult and thus the development of HLA–drug specific virtual screening procedures could become a key component of future precision medicine protocols. Importantly, the goal of this research was to test molecular docking as a computationally inexpensive virtual screening approach for the reliable prediction of critical drug–HLA-B*57:01 interactions triggering ADR events. As such, molecular dynamic simulations were not considered in this study, especially when considering the screening of drugs towards thousands of HLA variants.
Three X-ray crystals including abacavir bound to the antigen-binding cleft of HLA-B*57:01 with an endogenous peptide were downloaded from the Protein Data Bank: 3VRI (resolution 1.6 Å, peptide P1, Fig. 2), 3VRJ (resolution 1.9 Å, peptide P2), and 3UPR (resolution 2.0 Å, peptide P3) [16, 39]. The three crystal structures are highly similar and consist of bound abacavir covered by a differing endogenous peptide. In the case of 3VRI and 3VRJ, there are three distinct chains constructing the protein, chain A (275 residues for both crystals), chain B (100 and 99 residues, respectively), and chain C corresponding to the peptides P1 and P2, respectively (10 residues). Crystal 3UPR is a dimer with matching chains A and C (275 residues each), an interlinking chain B and D connecting chains A and C (99 residues), and chains P and Q corresponding to the peptide P3 (9 residues). The binding pocket is located on chain A for 3VRI, 3VRJ, and 3UPR (as well as chain C for 3UPR because it is a dimer).
In order to identify any significant differences between the binding sites of 3VRI, 3VRJ, and 3UPR, we analyzed the overall protein structure, ligand conformation, and co-binding peptides in several different ways. First, a general all-atom alignment was performed between the protein structures as implemented in the Schrodinger Suite . In general, a measured RMSD value can be used to determine structural similarity between closely related proteins (as is the case with our three crystals of the HLA-B*57:01 protein) . Next, a peptide backbone alignment was performed on the three co-binding peptides present in addition to superimposition of bound abacavir from the three crystals. Additionally, overlay similarity scores were calculated between the three structures for the protein, co-binding peptides, and bound abacavir using Discovery Studio . Finally, a binding site-specific alignment was performed using residues within 5 Å of bound abacavir.
Next, the physical chemical characteristics of the binding pocket were explored using SiteMap [43, 44]. SiteMap characterizes the possible binding sites of a protein by analyzing several physical chemical properties, such as size, volume, exposure to solvent, hydrophobic and hydrophilic space, and H-donor/-acceptor ability [43, 44]. Using these descriptors, two scores are generated: Site Score (Sscore) and Drugability Score (Dscore). A binding pocket that is likely to bind a ligand will have an Sscore greater than 0.8 and a Dscore greater than 0.83 [43, 44]. Using the default SiteMap parameters (6 Å buffer region, a minimum of 15 site points, restrictive hydrophobicity, and a standard grid), the HLA-B*57:01 binding pockets were analyzed under three conditions. First, the binding pocket was analyzed in the presence of the co-binding peptide using abacavir as the reference ligand. Second, the ligand binding environment was analyzed in the absence of co-binding peptide with abacavir as the reference compound. Third, the peptide was used as the reference ligand to map the binding pocket. Analyzing the peptide binding pocket of HLA-B*57:01 under these three conditions afforded a detailed analysis of how the ligand and peptide could impact the binding environment.
Abacavir’s binding mode with HLA-B*57:01 occurs through an altered repertoire mechanism. The altered repertoire binding mechanism occurs when an antigen binds non-covalently to the HLA active site and then, an endogenous or exogenous peptide binds non-covalently across the active site [4, 9, 16, 39]. This prevents the antigen from exiting the binding cleft, while also serving as a signaling trigger to T-cells resulting in an immune system response. As such, the model developed for this study used an altered repertoire binding mechanism. A cartoon schematic of this mechanism is provided in Fig. 1.
Drugs used to construct test set for docking with their proposed HLA-association
Uric acid inhibitor
Hypercholesterolemia cardiac heart disease
Seizures bipolar disorder
Acute Generalized Exanthematous Pustutosis (AGEP)
DILI ALT concentration increase
Agranulocytosis aplastic anemia neutropemia
Prior to the modeling, we conducted a 3D alignment of 3VRI, 3VRJ, and 3UPR in order to evaluate any significant deviations between the protein structures, bound ligand, and peptides. Molecular docking was conducted using the three aforementioned X-ray crystals preprocessed and curated (e.g., removal of water, addition of explicit hydrogens) using the Protein Preparation Wizard from the Schrodinger Suite [37, 57–60]. Missing side chains were generated using Prime [59, 60] while the protonation states of each side chain were generated using EPIK at pH = 7 [57, 58]. Protein minimization was performed using the OPLS3 force field [61–64]. Internal and external receptor grid boxes of 10 × 10 × 10 and 20 × 20 × 20 Å were defined using abacavir bound in the antigen cleft. The optimized structures of the test set were then docked with Schrodinger’s GLIDE software using both SP and XP scoring functions with a rigid protein, flexible ligand, and rigid peptide (when docked with peptide) [65–68]. Due to the presence of a unique peptide for each X-ray crystal, we conducted the docking with and without a peptide covering the solvent-accessible surface of the antigen-binding pocket. Each docking result was analyzed by comparing the docking and eModel scores in addition to the analysis of the drug’s binding mode in the B*57:01 site. The docking score (DS) consists of a sum of the Glide Score, measured from the SP or XP scoring functions, and the measured EPIK state penalty; the eModel score (eM) is a measure of the ‘favorability’ of a docked pose [57, 58, 65–67]. The DS may be used for comparing different ligands, but the eM score is suitable only to rank different conformations of the same ligand and should not be used to compare different ligands. A drug was considered to be B*57:01 liable (active) if the two following empirical thresholds were met: First, the DS had to be at less than or equal to −7 kcal/mol and second, the eM score had to be less than or equal to −50 kcal/mol. These thresholds were previously used in virtual screening protocols for discerning micromolar binders [69–71]. However, these scoring thresholds are specific for our model using GLIDE docking with SP and XP scoring functions and are obviously subject to change depending upon the studied protein, the software and method employed for the virtual screening. Our docking results were also evaluated using accuracy, sensitivity, specificity, positive performance value (PPV), and negative performance value (NPV) ; these values can be found in Additional file 1: Table S1.
Alignment of 3VRI, 3VRJ, and 3UPR
Even though these three peptides contain different residues, they do possess similar physical chemical attributes. For example, with the exception of P2, the ends of each peptide possess a basic residue and a hydrophobic residue (arginine and isoleucine for P1, histidine and valine for P3). Peptide P2 has hydrophobic residues at either extremities (leucine and isoleucine). Additionally, it should be noted that the center of each peptide contains a hydrophobic residue next to a hydrophilic residue. Indeed, peptide P1 has a leucine next to a glutamate at residue positions five and six, P2 has a leucine next to a tyrosine for residues five and six, and P3 has a tyrosine at position six with both a hydrophilic threonine at position four and hydrophobic leucine at position six.
For all three peptides, the AA residues at either end form non-covalent interactions with the binding pocket to anchor the peptide into the pocket. The carbonyl backbone for VAL, TYR, and ILE of P1 serve as H-bond acceptors for the in pocket residues TRP147 and TYR84 in addition to the formation of a salt bridge between ILE carboxylate group and LYS146. On the other end of P1, ARG serves as a H-bond donor with the pocket residues GLU63 and TYR171 while also forming a salt bridge with TYR59. Peptides P2 and P3 are anchored in a similar fashion in crystals 3VRJ and 3UPR. Importantly, the centroid of each peptide does not form interactions with binding pocket due to displacement by abacavir. The 2D-binding modes of peptides P1, P2, and P3 are provided in Additional file 1: Fig. S1.
Furthermore, the physical binding environment of the binding pocket was evaluated using SiteMap [43, 44]. The binding pocket for each crystal was evaluated under three conditions: (1) abacavir was used as the reference ligand with peptide present, (2) abacavir was used as the reference ligand in the absence of a co-binding peptide, and (3) the co-binding peptide was used as the reference ligand with abacavir present. Under conditions 1 and 2, the measured Sscore and Dscore were ranging between 1.1 and 1.3 indicating that the binding environment is extremely favorable (>0.8). Overall, condition 2 afforded slightly lower Sscore and Dscores, which is most likely due to increased solvent exposure due to the missing co-binding peptide. Interestingly, when the co-binding peptide was used as the reference ligand, the Sscore and Dscore were significantly lowered to a range of 0.8–1.1. Additionally, SiteMap identified two binding locations under the third condition (at either end of the peptide), while the center of the peptide was excluded from the binding surface. The observation of two binding pockets occurs as a result of the altered repertoire binding mechanism of abacavir with HLA-B*57:01: Abacavir displaced the center of the peptide from binding with the pocket. All SiteMap-generated binding surfaces are provided in Additional file 1: Figure S2.
Self-docking of abacavir in 3VRI with (+) and without (−) the presence of a co-binding peptide
The next step of our analysis was dedicated to the self-docking of abacavir in both the presence and absence of a co-binding peptide in B*57:01. To do so, we removed the native pose of abacavir from 3VRI, then we re-docked abacavir using both Glide SP and XP scoring functions. This self-docking procedure was conducted in order to test whether molecular docking could accurately reproduce the native binding mode of abacavir and investigate the significance of the co-binding peptide. To do so, we aligned the highest scoring conformation of self-docked abacavir with the native pose of abacavir from 3VRI. Self-docking without P1 was also performed; however, there is limited existing data about the potential binding mode of abacavir without peptide [16, 39]. When attempting to solve the X-ray crystals 3VRI and 3VRJ, Illing et al.  used molecular docking to probe the binding cleft of HLA-B*57:01 to assist their crystallization procedure. Similarly, when Ostrov et al.  solved for the 3UPR crystal, molecular docking was employed to select an optimized co-binding peptide for crystallization. Next, both the measured DS and eM scores of self-docked abacavir (with and without P1) were analyzed followed by the description of the molecular interactions between abacavir and HLA-B*57:01.
DS and eM scores reported as absolute values for abacavir—B*57:01 docking with crystals 3VRI, 3VRJ, and 3UPR in the presence and absence of peptides P1, P2, and P3 using the SP and XP scoring functions
Docking a set of ADR-causing drugs
The test set of ADR-causing drugs with known or putative HLA-binding profiles (Table 1) were docked using both SP and XP scoring functions. Upon inspection, we were able to determine that the test set of drugs used in this study had distinct therapeutic classes and low 2D chemical similarities. All of the drugs from the test set had distinct therapeutic classes with the exception of flucloxacillin and minocycline, which were both antibacterial agents. Furthermore, when the MACCS fingerprints  were calculated, the Tanimoto similarities  were quite distinct with the highest similarity (0.63) between the drugs atorvastatin and ciprofloxacin, and clozapine and ticlopidine, respectively. Abacavir and fenofibrate had the lowest measured structural similarity of 0.15. Surprisingly, the drugs abacavir, flucloxacillin, and pazopanib were all dissimilar with Tanimoto’s ranging from 0.35 to 0.40 even though these drugs are HLA-B*57:01 liable [15, 16, 39, 53]. All pairwise Tanimoto similarities are provided in Additional file 1: Table S2.
Virtual screening was conducted using the three X-ray crystals 3VRI, 3VRJ, and 3UPR in the presence and absence of peptide. Herein, we reported the results of docking using crystal 3VRI with and without peptide P1 (the results pertaining to docking performed with 3VRJ (P2) and 3UPR (P3) are provided in the Additional file 1). Statistical measures of our model’s ability to forecast HLA-B*57:01 liable drugs have been provided in Additional file 1: Table S1. Interestingly, our model performed very well in the absence of peptides P1, P2, and P3 [SEN = 0.67 (SP) and 1.00 (XP)], but when the peptide was included in docking, there was a significant decrease in our model’s ability to predict true positives (SEN = 0.33 for SP and XP). However, our model’s ability to forecast true negatives was extremely good as indicated by the high specificity (0.73–0.82 for SP and XP without peptide and 0.91–1.00 for SP and XP with peptide). The accuracy, positive prediction value (PPV), and negative prediction value (NPV) were also measured and are provided in Additional file 1: Table S1.
Interestingly, our model identified allopurinol and clozapine as two inactive compounds, i.e., failed to meet either DS or eM score criteria at the B*57:01 variant. Allopurinol is a HLA-B*58:01 active compound [48, 49] and was not expected to be active at the B*57:01 variant. Previous work had suggested that clozapine might be a HLA-B*39 or -B*57:01 active drug, but it was unclear which specific B-variant was preferred [26, 27]. In our current model, clozapine afforded a DS above the threshold of −7 kcal/mol and an eM greater than −50 kcal/mol providing further evidence that clozapine is not a B*57:01 binder. In some cases, our model produced conflicting results that would predict a drug as being inactive or active under varying conditions. For example, as shown in Fig. 7a, carbamazepine afforded DS well below the threshold (with the exception of GXP without P1); however, the measured eM score of carbamazepine were greater than our −50 kcal/mol threshold indicating that carbamazepine was not an HLA-B*57:01 liable drug (Fig. 7b).
In other instances, the model identified a drug as a B*57:01 active in the absence of peptide, but when the P1 peptide was present in the cleft, the drug would fail to meet either DS, eM, or both criteria. The drugs atorvastatin (HLA-DRB1*10:10 active), ciprofloxacin (HLA-B*50:02), and minocycline (HLA-B*35:02) were among these cases (Table 1). It was also observed that for some drugs, such as flucloxacillin, the use of SP or XP scoring function influenced the models prediction. The DS and eM scores for 3VRJ and 3UPR with the test set are provided in Additional file 1: Figures S7 and S8, respectively.
Pearson correlation coefficients were calculated for DS and eM scores obtained by the three X-ray crystal docking systems (Additional file 1: Tables S3, S4). The DS results without peptide showed a significant improvement in fit when the XP scoring function was used between 3VRI, 3VRJ, and 3UPR (0.71 ≤ RSP(−)P ≤ 0.88 and 0.92 ≤ RXP(−)P ≤ 0.98); however, when docking was performed in the presence of peptides P1, P2, or P3, there was more fluctuation in the results (0.57 ≤ RSP(+)P ≤ 0.87 and 0.70 ≤ RXP(+)P ≤ 0.74). Interestingly, when we compared eM scores, the similarity between docking grids appeared to be more dependent upon the scoring function employed as opposed to the influence of the peptide. When the SP function was used, the crystals showed excellent agreement in the absence of peptide, but in the presence of peptide there was more variation between the systems (0.73 ≤ RSP(−)P ≤ 0.84 and 0.60 ≤ RSP(+)P ≤ 0.91). Yet, when the XP function was used the overall fit improved in the absence and presence of peptide (0.60 ≤ RXP(−)P ≤ 0.91 and 0.67 ≤ RXP(+)P ≤ 0.83). For both DS and eM scores, the presence of peptide lowered the correlation, because not all of the test set drugs had a favorable conformation in the binding pocket with peptide. The discrepancy observed between scoring functions for eM correlations could result from differences between the SP and XP scoring functions in calculating conformational energy of the binding drug [65–68]. Importantly, the eM score was not used to compare results between different sets of ligands, but to measure how the SP and XP scoring functions performed when different crystals (3VRI, 3VRJ, and 3UPR) under different stresses (the presence or absence of peptides P1, P2, and P3) were used. Essentially, this provides a baseline for measuring the differences between our systems.
There were two top scoring drugs when docking with 3VRJ (see Additional file 1: Figures S9): abacavir docked using the SP function with peptide P2 and pazopanib docked using the XP function without peptide P2. When docking was performed with 3UPR, the top performing drug was again abacavir (the SP and XP scoring functions with peptide P3 both yielded the best DS and eM scores in this case). The plot of eM versus DS for 3UPR is provided in Additional file 1: Figures S10.
We should underline the results obtained in our docking study seem to be dependent on the composition of the co-binding peptides and this can be seen as a clear limitation. Thus, more analysis would be needed to study whether peptide P1 is rather specific to abacavir and therefore does not allow favorable interactions for the binding modes of other drugs (such as flucloxacillin or pazopanib which are both HLA-B*57:01 liable, but were predicted as inactive in the presence of peptide). Thus, HLA-variants may bind drugs with drug-specific (or class-specific) peptides that could significantly enable, boost, and thus impact their binding modes. Another possibility is related to the fact that P1 probably adopts a specific conformation favorable for abacavir binding. It means that other drugs may bind differently in the B*57:01 pocket and these binding modes are not favored in P1’s conformation from 3VRI. Thus molecular dynamic simulations  are needed to test and potentially confirm this hypothesis at the case by case level. Of course, MD simulations cannot be used for screening purposes, especially considering a large collection of drugs (e.g., DrugBank) towards thousands of HLA variants. Therefore, molecular docking should still be considered as the main approach for high throughput HLA–drug screening.
Additional focus on known HLA-B*57:01 actives: flucloxacillin and pazopanib
In addition to abacavir, the test set of compounds contained two other drugs, flucloxacillin and pazopanib, known to bind HLA-B*57:01 and potentially causing drug-induced liver injury (DILI) [15, 53]. Glide was able to successfully forecast pazopanib as a B*57:01 active drug with DS (DSSP = −8.7 kcal/mol, DSXP = −9.2 kcal/mol) and eM scores (eMSP = −78.0 kcal/mol, eMXP = −91.5 kcal/mol) well above the threshold in the absence of P1. Notably, the presence of P1 had a significant impact upon the DS and eM scores for pazopanib. When the SP function was used, pazopanib’s DS surpassed the threshold (DSSP = −9.1 kcal/mol), while the eM score failed to meet it (eMSP = −43.2 kcal/mol). However, when the XP function was used the opposite phenomenon was observed with the DS failing the threshold (DSXP = −6.6 kcal/mol) and the eM score surpassing it (eMXP = −56.3 kcal/mol). As a result, the binding mode of pazopanib was investigated for SP and XP scoring functions with and without peptide P1.
Finally, the binding mode of flucloxacillin was further investigated using the XP scoring function with and without peptide P1 (see Fig. 10b, c). There are very significant differences in the binding mode indicating that the binding location in B*57:01 most likely occurs in a different region of the pocket than where the docking grid was generated using abacavir. More specifically, in the absence of P1, H-bonding is observed between ASH114 and the carboxylate group and between SER116 and the amide linker; π–π stacking also occurs between TYR123 and the isoxazol moiety (Fig. 10b). However, when docking with P1 the ligand-residue interactions change drastically (in addition to the global orientation of the ligand). As shown in Fig. 10c, the H-bonding of SER116 now occurs with the carboxylate group, whereas π–π stacking occurs between TYR99 and the 2-chloro-6-fluoro-phenyl group. Consequently, the position of flucloxacillin has shifted in the protein pocket. These significant variations in flucloxacillin’s binding modes indicate that molecular dynamic studies will be needed in order to further investigate the potential binding mode(s) between B*57:01 and flucloxacillin.
Reliably predicting whether a drug candidate is likely to be an HLA binder at a given variant constitutes a potentially valuable knowledge in evaluating the likelihood of an associated ADR event in a subpopulation of patients. Future work will focus on the B*57:01 screening for the entire DrugBank database  which includes over 7000 drugs and drug candidates to identify unknown drug–HLA interactions. Since we showed that endogenous peptides have a significant impact in the binding mode of drugs with HLA, it is still unclear how the peptide changes conformations in response to different drugs, and this is a clear limitation of molecular docking with rigid peptides. Therefore, we plan to conduct peptide self- and cross-docking  in the HLA binding cleft as well as molecular dynamics simulations with peptides P1, P2, and P3 to monitor and analyze the dynamic variations of the binding interactions between different peptides, different drugs, and B*57:01. However, as demonstrated in this study, molecular docking represents a fast and reliable cheminformatics technique to forecast drug–HLA interactions. This is critical to screen large datasets of chemicals towards thousands of HLA variants and determine whether a particular drug has a ‘best-binding’ peptide partner that is unique given a particular HLA variant.
GV conducted all the calculations and wrote the manuscript. DF conceived the study and wrote the manuscript. Both authors read and approved the final manuscript.
GV thanks the NCSU Department of Chemistry. DF sincerely thanks the NCSU Chancellor’s Faculty Excellence Program, the NCSU Bioinformatics Research Center (BRC), and the NCSU's department of Chemistry.
Both authors declare that they have no competing financial interests.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Edwards IR, Aronson JK (2000) Adverse drug reactions: definitions, diagnosis, and management. Lancet 356:1255–1259View ArticleGoogle Scholar
- U.S. Food and Drug Administration (2016) FDA adverse event reporting system (FAERS) statistics [Internet]. http://www.fda.gov/Drugs/GuidanceComplianceRegulatoryInformation/Surveillance/AdverseDrugEffects/ucm070093.htm. Cited 29 Mar 2016
- Hunziker T, Bruppacher R, Kuenzi UP, Maibach R, Braunschweig S, Halter F et al (2002) Classification of ADRs: a proposal for harmonization and differentiation based on the experience of the Comprehensive Hospital Drug Monitoring Bern/St. Gallen, 1974–1993. Pharmacoepidemiol Drug Saf 11:159–163View ArticleGoogle Scholar
- Bharadwaj M, Illing P, Theodossis A, Purcell AW, Rossjohn J, McCluskey J (2012) Drug hypersensitivity and human leukocyte antigens of the major histocompatibility complex. Annu Rev Pharmacol Toxicol 52:401–431View ArticleGoogle Scholar
- Daly AK (2014) Human leukocyte antigen (HLA) pharmacogenomic tests: potential and pitfalls. Curr Drug Metab 15:196–201View ArticleGoogle Scholar
- Yip VLM, Alfirevic A, Pirmohamed M (2015) Genetics of immune-mediated adverse drug reactions: a comprehensive and clinical review. Clin Rev Allergy Immunol 48:165–75Google Scholar
- Marsh SGE, Parham P, Barber LD (2000) The HLA facts book. Academic Press, London, pp 55–57Google Scholar
- Madden DR (1995) The three-dimensional structure of peptide-MHC complexes. Annu Rev Immunol 13:587–622View ArticleGoogle Scholar
- Pichler WJ, Beeler A, Keller M, Lerch M, Posadas S, Schmid D et al (2006) Pharmacological interaction of drugs with immune receptors: the p–i concept. Allergol Int 55:17–25View ArticleGoogle Scholar
- Robinson J, Hayhurst J, Flicek P, Parham P, Marsch SGE (2015) The IPD and IMGT/HLA database: allele variant databases. Nucleic Acids Res 43:D423–D431View ArticleGoogle Scholar
- Cao K, Hollenbach J, Shi X, Shi W, Chopek M, Fernández-Viña MA (2001) Analysis of the frequencies of HLA-A, B, and C alleles and haplotypes in the five major ethnic groups of the United States reveals high levels of diversity in these loci and contrasting distribution patterns in these populations. Hum Immunol 62:1009–1030View ArticleGoogle Scholar
- Chung WH, Hung S-I, Hong H-S, Hsih M-S, Yang L-C, Ho H-C et al (2004) A marker for Stevens–Johnson syndrome. Nature 428:6–7View ArticleGoogle Scholar
- Wu XT, Hu FY, An DM, Yan B, Jiang X, Kwan P et al (2010) Association between carbamazepine-induced cutaneous adverse drug reactions and the HLA-B*1502 allele among patients in central China. Epilepsy Behav 19:405–408View ArticleGoogle Scholar
- Mockenhaupt M, Messenheimer J, Tennis P, Schlingmann J (2005) Risk of Stevens–Johnson syndrome and toxic epidermal necrolysis in new users of antiepileptics. Neurology 64:1134–1138View ArticleGoogle Scholar
- Daly AK, Donaldson PT, Bhatnagar P, Shen Y, Pe’er I, Floratos A et al (2009) HLA-B*5701 genotype is a major determinant of drug-induced liver injury due to flucloxacillin. Nat Genet 41:816–819View ArticleGoogle Scholar
- Illing PT, Vivian JP, Dudek NL, Kostenko L, Chen Z, Bharadwaj M et al (2012) Immune self-reactivity triggered by drug-modified HLA–peptide repertoire. Nature 486:554–558Google Scholar
- Liu R, Yu X, Wallqvist A (2015) Data-driven identification of structural alerts for mitigating the risk of drug-induced human liver injuries. J Cheminform 7:1–8Google Scholar
- Vilar S, Hripcsak G (2016) Leveraging 3D chemical similarity, target and phenotypic data in the identification of drug-protein and drug-adverse effect associations. J Cheminform 8:35View ArticleGoogle Scholar
- Low YS, Caster O, Bergvall T, Fourches D, Zang X, Norén GN et al (2015) Cheminformatics-aided pharmacovigilance: application to Stevens-Johnson Syndrome. J Am Med Inform Assoc 0:1–11Google Scholar
- Rodgers AD, Zhu H, Fourches D, Rusyn I, Tropsha A, Hill C et al (2010) Modeling liver-related adverse effects of drugs using k nearest neighbor quantitative structure—activity relationship method. Chem Res Toxicol 23:724–732View ArticleGoogle Scholar
- Chen B, Dong X, Jiao D, Wang H, Zhu Q, Ding Y et al (2010) Chem2Bio2RDF: a semantic framework for linking and data mining chemogenomic and systems chemical biology data. BMC Bioinform 11:255View ArticleGoogle Scholar
- Low Y, Sedykh A, Fourches D, Golbraikh A, Whelan M, Rusyn I et al (2013) Integrative chemical–biological read-across approach for chemical hazard classification. Chem Res Toxicol 26:1199–1208View ArticleGoogle Scholar
- Luo H, Ye H, Ng H, Shi L, Tong W, Mattes W et al (2015) Understanding and predicting binding between human leukocyte antigens (HLAs) and peptides by network analysis. BMC Bioinform 16(Suppl 1):S9View ArticleGoogle Scholar
- Tingting Du, Yang L, Luo H, Zhou P, Mei H, Xuan J et al (2015) HLADR: a database system for enhancing the discovery of biomarkers for predicting human leukocyte antigen-mediated idiosyncratic adverse drug reactions. Biomark Med 9:1079–1093View ArticleGoogle Scholar
- Yang C, Wang C, Zhang S, Huang J, Zhou P (2015) Structural and energetic insights into the intermolecular interaction among human leukocyte antigens, clinical hypersensitive drugs and antigenic peptides. Mol Simul 41:741–751View ArticleGoogle Scholar
- Goldstein JI, Jarskog LF, Hilliard C, Alfirevic A, Duncan L, Fourches D et al (2014) Clozapine-induced agranulocytosis is associated with rare HLA-DQB1 and HLA-B alleles. Nat Commun 5:4757View ArticleGoogle Scholar
- Legge SE, Hamshere ML, Ripke S, Pardinas AF, Goldstein JI, Rees E et al (2016) Genome-wide common and rare variant analysis provides novel insights into clozapine-associated neutropenia. Mol Psychiatry 0:1–7Google Scholar
- Alfirevic A, Gonzalez-Galarza F, Bell C, Martinsson K, Platt V, Bretland G et al (2012) In silico analysis of HLA associations with drug-induced liver injury: use of a HLA-genotyped DNA archive from healthy volunteers. Genome Med 4:51View ArticleGoogle Scholar
- Schotland P, Bojunga N, Zien A, Trame MN, Lesko LJ (2016) Improving drug safety with a systems pharmacology approach. Eur J Pharm Sci 94:84–92View ArticleGoogle Scholar
- Wei CY, Chung WH, Huang HW, Chen YT, Hung SI (2012) Direct interaction between HLA-B and carbamazepine activates T cells in patients with Stevens–Johnson syndrome. J Allergy Clin Immunol 129:1562–1569View ArticleGoogle Scholar
- Paul S, Sidney J, Sette A, Peters B (2016) TepiTool: a pipeline for computational prediction of T cell epitope candidates. Curr Protocol Immunol 114:18.19.1–18.19.24Google Scholar
- Paul S, Lindestam Arlehamn CS, Scriba TJ, Dillon MBC, Oseroff C, Hinz D et al (2015) Development and validation of a broad scheme for prediction of HLA class II restricted T cell epitopes. J Immunol Methods 422:28–34View ArticleGoogle Scholar
- Zhou P, Zhang S, Wang Y, Yang C, Huang J (2016) Structural modeling of HLA-B*1502/peptide/carbamazepine/T-cell receptor complex architecture: implication for the molecular mechanism of carbamazepine-induced Stevens–Johnson syndrome/toxic epidermal necrolysis. J Biomol Struct Dyn 34:1806–1817View ArticleGoogle Scholar
- Metushi IG, Wriston A, Banerjee P, Gohlke BO, English AM, Lucas A et al (2015) Acyclovir has low but detectable influence on HLA-B*57:01 specificity without inducing hypersensitivity. PLoS ONE 10:e0124878View ArticleGoogle Scholar
- Irwin JJ, Sterling T, Mysinger MM, Bolstad ES, Coleman RG (2012) ZINC: a free tool to discover chemistry for biology. J Chem Inf Model 52:1757–1768View ArticleGoogle Scholar
- Lucas A, Lucas M, Strhyn A, Keane NM, McKinnon E, Pavlos R et al (2015) Abacavir-Reactive Memory T Cells Are Present in Drug Naïve Individuals. PLoS ONE 10:e0117160View ArticleGoogle Scholar
- Madhavi Sastry G, Adzhigirey M, Day T, Annabhimoju R, Sherman W (2013) Protein and ligand preparation: parameters, protocols, and influence on virtual screening enrichments. J Comput Aided Mol Des 27:221–234View ArticleGoogle Scholar
- Convertino M, Dokholyan NV (2016) Computational modeling of small molecule ligand binding interactions and affinities. Methods Mol Biol 1414:23–32View ArticleGoogle Scholar
- Ostrov D, Grant BJ, Pompeu Y, Sidney J, Harndahl M, Southwood S et al (2012) Drug hypersensitivity caused by alteration of the MHC-presented self-peptide repertoire. Proc Natl Acad Sci 109:9959–9964View ArticleGoogle Scholar
- Yang A-S, Honig B (2000) An integrated approach to the analysis and modeling of protein sequences and structures. I. Protein structural alignment and a quantitative measure for protein structural distance. J Mol Biol 301:665–678View ArticleGoogle Scholar
- Mizuguchi K, Go N (1995) Seeking significance in three-dimensional protein structure comparisons. Curr Opin Struct Biol 5:377–382View ArticleGoogle Scholar
- Dassault Systemes BIOVIA (2016) Discovery studio modeling environment. Dassault Systemes, San DiegoGoogle Scholar
- Halgren T (2007) New method for fast and accurate binding-site identification and analysis. Chem Biol Drug Des 69:146–148View ArticleGoogle Scholar
- Halgren TA (2009) Identifying and characterizing binding sites and assessing druggability. J Chem Inf Model 49:377–389View ArticleGoogle Scholar
- Berka N, Gill JM, Liacini A, O’Bryan T, Khan FM (2012) Human leukocyte antigen (HLA) and pharmacogenetics: screening for HLA-B*57:01 among human immunodeficiency virus-positive patients from southern Alberta. Hum Immunol 73:164–167View ArticleGoogle Scholar
- Martin AM, Nolan D, Gaudieri S, Almeida CA, Nolan R, James I et al (2004) Predisposition to abacavir hypersensitivity conferred by HLA-B*5701 and a haplotypic Hsp70-Hom variant. Proc Natl Acad Sci USA 101:4180–4185View ArticleGoogle Scholar
- Saag M, Balu R, Phillips E, Brachman P, Martorell C, Burman W et al (2008) High sensitivity of human leukocyte antigen-b*5701 as a marker for immunologically confirmed abacavir hypersensitivity in white and black patients. Clin Infect Dis 46:1111–1118View ArticleGoogle Scholar
- Hung S, Chung W, Liou L-B, Chu CC-C, Lin M, Huang H-P et al (2005) HLA-B*5801 allele as a genetic marker for severe cutaneous adverse reactions caused by allopurinol. Proc Natl Acad Sci USA 102:4134–4139View ArticleGoogle Scholar
- Park HJ, Kim YJ, Kim DH, Kim J, Park KH, Park J-W et al (2016) HLA allele frequencies in 5802 Koreans: varied allele types associated with SJS/TEN according to culprit drugs. Yonsei Med J 57:118View ArticleGoogle Scholar
- Genin E, Chen D-P, Hung S-I, Sekula P, Schumacher M, Chang P-Y et al (2014) HLA-A*31:01 and different types of carbamazepine-induced severe cutaneous adverse reactions: an international study and meta-analysis. Pharmacogenomics J. 14:281–288View ArticleGoogle Scholar
- Athanasiou MC, Dettling M, Cascorbi I, Mosyagin I, Salisbury BA, Pierz KA et al (2011) Candidate gene analysis identifies a polymorphism in HLA-DQB1 associated with clozapine-induced agranulocytosis. J Clin Psychiatry 72:458–463View ArticleGoogle Scholar
- Dettling M, Cascorbi I, Opgen-Rhein C, Schaub R (2007) Clozapine-induced agranulocytosis in schizophrenic Caucasians: confirming clues for associations with human leukocyte class I and II antigens. Pharmacogenomics J 7:325–332View ArticleGoogle Scholar
- Xu C-F, Johnson T, Wang X, Carpenter C, Graves AP, Warren L et al (2016) HLA-B∗57:01 confers susceptibility to pazopanib-associated liver injury in patients with cancer. Clin Cancer Res 22:1371–1377View ArticleGoogle Scholar
- Hirata K, Takagi H, Yamamoto M, Matsumoto T, Nishiya T, Mori K et al (2008) Ticlopidine-induced hepatotoxicity is associated with specific human leukocyte antigen genomic subtypes in Japanese patients: a preliminary case-control study. Pharmacogenomics J 8:29–33View ArticleGoogle Scholar
- Durant JL, Leland BA, Henry DR, Nourse JG (2002) Reoptimization of MDL keys for use in drug discovery. J Chem Inf Comput Sci 42:1273–1280View ArticleGoogle Scholar
- Willett P, And JMB, Downs GM (1998) Chemical similarity searching. J Chem Inf Comput Sci 38:983–996View ArticleGoogle Scholar
- Greenwood JR, Calkins D, Sullivan AP, Shelley JC (2010) Towards the comprehensive, rapid, and accurate prediction of the favorable tautomeric states of drug-like molecules in aqueous solution. J Comput Aided Mol Des 24:591–604View ArticleGoogle Scholar
- Shelley JC, Cholleti A, Frye LL, Greenwood JR, Timlin MR, Uchimaya M (2007) Epik: a software program for pKa prediction and protonation state generation for drug-like molecules. J Comput Aided Mol Des 21:681–691View ArticleGoogle Scholar
- Jacobson MP, Friesner RA, Xiang Z, Honig B (2002) On the role of the crystal environment in determining protein side-chain conformations. J Mol Biol 320:597–608View ArticleGoogle Scholar
- Jacobson MP, Pincus DL, Rapp CS, Day TJF, Honig B, Shaw DE et al (2004) A hierarchical approach to all-atom protein loop prediction. Proteins Struct Funct Genet 55:351–367View ArticleGoogle Scholar
- Harder E, Damm W, Maple J, Wu C, Reboul M, Xiang JY et al (2016) OPLS3: a force field providing broad coverage of drug-like small molecules and proteins. J Chem Theory Comput 12:281–296View ArticleGoogle Scholar
- Jorgensen WL, Maxwell DS, Tirado-Rives J (1996) Development and testing of the OLPS all-atom force field on conformational energetics and properties of organic liquids. J Am Chem Soc 118:11225–11236View ArticleGoogle Scholar
- Jorgensen WL, Tirado-Rives J (1988) The OPLS [optimized potentials for liquid simulations] potential functions for proteins, energy minimizations for crystals of cyclic peptides and crambin. J Am Chem Soc 110:1657–1666View ArticleGoogle Scholar
- Shivakumar D, Williams J, Wu Y, Damm W, Shelley J, Sherman W (2010) Prediction of absolute solvation free energies using molecular dynamics free energy perturbation and the opls force field. J Chem Theory Comput 6:1509–1519View ArticleGoogle Scholar
- Friesner RA, Banks JL, Murphy RB, Halgren TA, Klicic JJ, Mainz DT et al (2004) Glide: a new approach for rapid, accurate docking and scoring. 1. Method and assessment of docking accuracy. J Med Chem 47:1739–1749View ArticleGoogle Scholar
- Friesner RA, Murphy RB, Repasky MP, Frye LL, Greenwood JR, Halgren TA et al (2006) Extra precision glide: docking and scoring incorporating a model of hydrophobic enclosure for protein-ligand complexes. J Med Chem 49:6177–6196View ArticleGoogle Scholar
- Halgren TA, Murphy RB, Friesner RA, Beard HS, Frye LL, Pollard WT et al (2004) Glide: a new approach for rapid, accurate docking and scoring. 2. Enrichment factors in database screening. J Med Chem 47:1750–1759View ArticleGoogle Scholar
- Tubert-Brohman I, Sherman W, Repasky M, Beuming T (2013) Improved docking of polypeptides with glide. J Chem Inf Model 53:1689–1699View ArticleGoogle Scholar
- Shityakov S, Förster C (2014) In silico structure-based screening of versatile P-glycoprotein inhibitors using polynomial empirical scoring functions. Adv Appl Bioinform Chem 7:1–9Google Scholar
- Shityakov S, Förster C (2014) In silico predictive model to determine vector-mediated transport properties for the blood-brain barrier choline transporter. Adv Appl Bioinform Chem 7:23–36Google Scholar
- Fourches D, Muratov E, Ding F, Dokholyan NV, Tropsha A (2013) Predicting binding affinity of CSAR ligands using both structure-based and ligand-based approaches. J Chem Inf Model 53:1915–1922View ArticleGoogle Scholar
- Fawcett T (2006) An introduction to ROC analysis. Pattern Recognit 27:861–874View ArticleGoogle Scholar
- Proctor EA, Dokholyan NV (2016) Applications of discrete molecular dynamics in biology and medicine. Curr Opin Struct Biol 37:9–13View ArticleGoogle Scholar
- Wishart DS, Knox C, Guo AC, Shrivastava S, Hassanali M, Stothard P et al (2006) DrugBank: a comprehensive resource for in silico drug discovery and exploration. Nucleic Acids Res 34:D668–D672View ArticleGoogle Scholar