Identification and characterisation of putative drug binding sites in human ATP-binding cassette B5 (ABCB5) transporter.

The human ATP-binding cassette B5 (ABCB5) transporter, a member of the ABC transporter superfamily, is linked to chemoresistance in tumour cells by drug effluxion. However, little is known about its structure and drug-binding sites. In this study, we generated an atomistic model of the full-length human ABCB5 transporter with the highest quality using the X-ray crystal structure of mouse ABCB1 (Pgp1), a close homologue of ABCB5 and a well-studied member of the ABC family. Molecular dynamics simulations were used to validate the atomistic model of ABCB5 and characterise its structural properties in model cell membranes. Molecular docking simulations of known ABCB5 substrates such as taxanes, anthracyclines, camptothecin and etoposide were then used to identify at least three putative binding sites for chemotherapeutic drugs transported by ABCB5. The location of these three binding sites is predicted to overlap with the corresponding binding sites in Pgp1. These findings will serve as the basis for future in vitro studies to validate the nature of the identified substrate-binding sites in the full-length ABCB5 transporter.

Introduction

Human ATP-binding cassette B5 (ABCB5) transporter is an integral membrane protein and a member of the ATP-binding cassette (ABC) transporter superfamily. ABCB5 transporters are predominantly expressed in pigmented cells such as melanocytes and retinal epithelial cells. However, the physiological relevance of these transporters in these cells remains unknown. ABCB5 exists in four isoforms, formed by alternative splicing: isoform 1 (ABCB5α; ABCB5.f), isoform 2 (ABCB5β; ABCB5.a), isoform 3 (ABCB5.e) and isoform 4 (ABCB5 full-length; ABCB5.ts) [9], [27]. Of these four isoforms, the ABCB5β isoform (817 amino acids) and full-length (1257 amino acids) have been identified to function as drug efflux pumps by in vitro studies [29], [44]. The ABCB5α and ABCB5.e isoforms are considered non-functional due to the absence of an intact domain that facilitates drug binding [9]; the significance of these two isoforms is unknown.

Expression of ABCB5β and its full length form have been reported to occur in numerous cancer types such as colon, liver and breast cancer, leukaemia, head and neck, oral squamous cell and merkel cell carcinoma, [27], [28], [32], [33], [48]. In tumour cells, the expression of ABCB5 is linked to multidrug resistance (MDR) [10]. ABCB5 transporter extrudes anti-cancer drugs accumulated within tumour cells, conferring resistance, in a similar fashion to ATP-binding cassette B1 (ABCB1; P-glycoprotein 1; Pgp1) transporter [10], [28], [44].

The drug efflux activity of ABCB5 was first reported by Frank et al. showing that MCF-7 cells (breast carcinoma) expressing ABCB5β isoform had an increased efflux of rhodamine 123 [29]. In a follow up study, decreased accumulation of doxorubicin and 5-fluorouracil was observed in ABCB5β-enriched melanoma G3361 cells compared to the ABCB5β low population [28]. Later, Kawanobe et al. demonstrated the efflux of several chemotherapeutic drugs in HEK293 cells overexpressing the full-length ABCB5 isoform [44]. Furthermore, another study compared the efflux action of the ABCB5β and full-length ABCB5 isoforms in yeast cells (S. cerevisiae)overexpressing the two isoforms. The study reported that full-length ABCB5 was functional and conferred resistance for rhodamine 123, daunorubicin and clorgyline; however, ABCB5β-overexpressing cells did have resistance [46]. As chemoresistance rendered by full-length ABCB5 is now established by multiple studies, characterising its drug binding sites will shed more light on the substrate specificity of the transporter.

A typical functional human ABC full transporter such as Pgp1 comprises of two transmembrane domains (TMDs), each comprised of six α-helices, and two cytosolic nucleotide-binding domains (NBDs) interconnected by a linker region (Fig. 1). The structural topology of the full-length ABCB5 isoform (1257 amino acids in a single polypeptide chain) comprises two TMDs and two NBDs interconnected by a linker region (Fig. 1), similar to Pgp1 [59].

Fig. 1

Predicted structural topology of full-length human ABCB5 isoform. The topology of the full-length isoform is comprised of two TMDs and two cytosolic NBDs interconnected by a linker region. Each TMD consists of six α-helices spanning in and out of the membrane connected by coupling helices (CH) and extracellular loops. TMD = transmembrane domain; NBD = nucleotide-binding domain.

Pgp1 and ABCB5 are close homologues and, in fact, the two transporters belong to the same subfamily ABCB. Pgp1 is the most studied transporter that acts as a MDR mediator in cancer. Numerous site-directed mutagenesis and structure-based studies have described the molecular mechanisms involved in substrate translocation, as well as the conformational changes that occur during this process. While ABC transporters are known to transition between open-inward and closed-outward conformations during the transport cycle, their open-inward facing conformation is the resting state conformation which allows molecules to interact with the transporter [2], [31]. Over the last decade, numerous structures of open-inward and closed-outward facing Pgp1 transporters were solved by X-ray diffraction crystallography and, more recently, by cryoelectron microscopy [1], [2], [26], [47], [49], [83]. These X-ray structures have revealed that crystallographic conditions such as the presence of different detergents, nucleotides or ligands can influence the conformation of the transporter [4], [67]. Extensive mutagenesis studies of Pgp1, through cysteine, alanine or arginine replacement approaches, led to the identification of residues in the TMD, NBD and loops implicated in the binding and translocation of its substrates [57], [53], [54], [55]. Biochemical (mutagenesis) and pharmacological studies, such as drug accumulation, ATPase and transport assays, have revealed that multiple, distinct binding sites exist within Pgp1 and are specific to the class of substrates being transported. Three pharmacological binding sites have been proposed based on competitive binding of different substrates, namely the R-site (rhodamine 123), the H-site (Hoechst 33342) and the vinblastine/prazosin binding site [56], [57], [78]. The R- and H-sites of Pgp1 interact in a positively cooperative manner, thus impacting substrate efflux [78]. In the presence of other substrates, the vinblastine site binds in a non-competitive manner, indicating that it has an additional binding site [56]. In the last two decades, numerous structural models of Pgp1 have been generated for experimental and computational studies [7], [63], [71], [72], [73] to characterise protein-drug interactions by molecular docking simulations [19], [63], [64], [73], [90]. Molecular docking simulations indeed confirmed the existence of R- and H- sites as substrate binding sites, while the vinblastine site was predicted to be modulator binding site [22], [19].

In this work, a high-quality, atomistic model of the human full-length ABCB5 transporter was generated using the crystal structure of mouse Pgp1 in its open-inward facing conformation. To assess the stability and conformations of this structure in a model cell membrane environment, molecular dynamics (MD) simulations were conducted. This was followed by molecular docking simulations with known substrates of ABCB5 to identify their putative drug binding sites. The predicted binding interactions of rhodamine 123 and paclitaxel were contrasted with site-directed mutagenesis data for rhodamine 123 and crystallographic data for paclitaxel in human Pgp1.

Materials and methods

Multiple sequence alignment

The FASTA sequence of mouse Pgp1 (UniProt ID: P21447-1), human Pgp1 (UniProt ID: PO8183-1) and full-length human ABCB5 (UniProt ID: Q2M3G0-4) were downloaded from the UniProt database (www.uniprot.org). Multiple sequence alignment was generated using Clustal Omega (ClustalO), with transition matrix Gonnet, using the default setting with a gap opening penalty of 6, and a gap extension of 1. Multiple sequence alignment of full-length human ABCB5 with the human Pgp1 and mouse Pgp1 FASTA sequences revealed a high sequence identity (~54%) and similarity (~74%) with both human and mouse Pgp1 sequences (Fig. S1 in the Supplementary Information), which is excellent for homology modelling [75]. As expected, a higher sequence identity was observed in the residues corresponding to the NBD regions of these three sequences, relative to the TMDs, as both ABCB5 and Pgp1 belong to the same ABCB subfamily.

Protein modelling

Initially, seventeen different structural models of human full-length ABCB5 transporter were generated by homology modelling using Phyre2 (Protein Homology/analogy Recognition Engine V2.0) [45], SWISS MODEL [6], [88] and MODELLER [23], as well as by threading using I-TASSER [91]. The FASTA sequence of full-length human ABCB5 (Q2M3G0-4) was submitted to Phyre2, and I-TASSER servers, whilst sequence alignment of ABCB5 and the template (PDB ID 4Q9H) was submitted to SWISS MODEL (alignment mode) [6], [88]. Template selection in Phyre2 is automated, and structural models were generated using the structures of Pgp1 from Caenorhabditis elegans (PDB ID 4F4C, with a resolution of 3.4 Å) [42] and mouse Pgp1 (PDB ID 3G60, with a resolution of 4.4 Å) [2] as target proteins. Both the normal and intensive modes in Phyre2 were used. The normal and intensive modes in Phyre2 work in a similar manner up to the point of sequence alignment. In the case of normal mode, the HH search alignment algorithm is used to generate backbone models followed by loop and side chains modelling. In the case of intensive mode, models are generated for each selected template using Poing, a multi-template modelling tool. Poing is used only in intensive mode to create complete models of input protein sequences from separate templates. If templates are not available then Poing uses ab initio modelling for the missing regions. The intensive mode in Phyre2 has been reported to perform a more refined model selection compared to the normal mode based on sequence identity, confidence and coverage [45]. Both modes were thus used for the generation of ABCB5 models.

The X-ray crystal structure of mouse Pgp1 (apo form; PDB ID 4Q9H, with a resolution of 3.4 Å) was chosen as a template for sequence alignment due to its relatively high structural resolution. Sequence alignment of human ABCB5 and mouse Pgp1 was generated using Clustal Omega and visualised using BioEdit (V7.2.5) [34]. All templates used in the modelling of the structure of ABCB5 were in the open-inward facing conformation, which has a high affinity for substrate binding.

The stereochemical quality of all models generated was assessed using Ramachandran plots (Rampage) [69], [70], QMEAN [5], ERRAT [11] and Verify-3D [16]. These programs estimate the quality of the protein structural model upon analyses of the torsional angles of amino acid residues, distance-dependent all-atom interaction potentials, non-bonded interactions between different atoms types, and provide a score by comparing these values with highly refined structures, or by comparing the 3D profile of the protein with its own amino acid sequence [5], [11], [69]. The best model was selected based on the above functions. The linker region of the best ABCB5 model was refined using the loop optimisation function in MODELLER [40].

Molecular dynamics simulations

Construction of lipid bilayer membrane systems

Model cell membrane systems consisting of two different phospholipid compositions were constructed: POPC-cholesterol (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine and cholesterol) and POPC-POPE-cholesterol (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine, 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoethanolamine and cholesterol) using MemBuilder II [30]. The POPC-cholesterol membrane was built as a symmetric lipid bilayer, composed of 160 POPC and 20 cholesterol molecules in an 80:20 wt% ratio (89:11 mol% ratio) in each leaflet. The POPC-POPE-cholesterol membrane was modelled as an asymmetric lipid bilayer composed of POPC, POPE and cholesterol in a 67:13:20 wt% ratio (72:15:13 mol% ratio). The POPC-POPE-cholesterol membrane consisted of 150 POPC and 30 cholesterol molecules in the outer leaflet and 110 POPC, 54 POPE and 16 cholesterol molecules in the inner leaflet. The POPC-cholesterol membrane system has been reported to be more suitable for molecular simulations with Pgp1 than the traditional DMPC-cholesterol (dimiristoylphoshatidylcholine) system [21], [82] because unsaturated lipids better mimic eukaryotic membranes. Since ABCB5 has been reported to be predominantly expressed in pigmented cells [9], a POPC-POPE-cholesterol membrane system was built to mimic the melanoma cell membrane [68]. The POPE levels in melanoma cells have been reported to be relatively low and located predominantly in the inner (cytoplasmic) leaflet [17]. Hence, all POPE molecules in the POPC-POPE-cholesterol membrane system were placed in the inner leaflet. Both membranes were solvated with 10,780 water molecules. In each system, 94Na+ and 94Cl− ions were added to give a final ionic strength of 150 mM NaCl [62]. Each membrane system was energy minimised, after which steric clashes from water and lipid molecules were removed with two short, sequential NVT (1 ns) and NPT (1 ns) simulations. This was followed by a production simulation of 400 ns. All simulation parameters were set as outlined in the simulation parameters subsection. The area per lipid (APL) and bilayer thickness were used to monitor the equilibration of the membranes. The final frames in each simulation trajectory were used to set up the protein-membrane systems, as described below.

Assembly of protein-membrane systems

Histidine (His) residues 145, 581 and 1212 in the ABCB5 model were protonated to give a positive charge to each as these residues were reported to be crucial for the stability of the protein in a lipid bilayer system [62]. Other His residues were modelled in their neutral form. After energy minimisation using the steepest descents algorithm, the protein was combined with each one of the pre-equilibrated membranes to construct two protein-membrane systems, ABCB5-POPC-cholesterol and ABCB5-POPC-POPE-cholesterol. The OPM database [51] was used to predict the position of the transmembrane helices (TMHs) in the lipid bilayer. The ABCB5-POPC-cholesterol protein-membrane system was solvated with 35,267 water molecules, and the g_membed tool in GROMACS was used to remove overlapping lipids, which resulted in a membrane composed of 155 POPC and 21 cholesterol molecules in the outer leaflet and 156 POPC and 17 cholesterol molecules in the inner leaflet. The net positive charge of the protein (6+) was neutralised with 6 Cl− ions and a further 157 Na+, 8 Mg+2 and 179 Cl− ions were added to give a final physiological ionic strength of 150 mM NaCl and a concentration of 5 mM MgCl2 to the ABCB5-POPC-cholesterol system.

Similarly, the ABCB5-POPC-POPE-cholesterol protein-membrane system was solvated with 33,802 water molecules, and overlapping lipids were removed as above, resulting in a membrane composed of 144 POPC and 32 cholesterol molecules in the outer leaflet and 107 POPC, 50 POPE and 13 cholesterol molecules in the inner leaflet. The net positive charge of the protein (6 + ) was neutralised with 6Cl− ions and a further 154 Na+, 8 Mg+2 and 176 Cl− ions were added to achieve the same ionic strength and concentration of MgCl2 as the ABCB5-POPC-cholesterol system. The systems were energy minimised using the steepest descents algorithm, after which the system was relaxed with two short, sequential NVT (1 ns) and NPT (1 ns) molecular dynamics (MD) simulations. Each system was then simulated under NPT conditions for at least 200 ns. Equilibration of the protein was monitored by measuring the root-mean-square deviation (RMSD) of the backbone atoms in the whole protein, the NBD and TMHs, using the starting structure as the reference. Standard deviations of RMSD values and two-sided confidence intervals (CI) were used to assess whether differences in RMSD between different protein domains were statistically significant. For each protein-membrane system, the dominant conformations of the protein were identified by RMSD clustering analysis using the trajectory frames from the final 50 ns of the corresponding simulation, using all atoms in the TM domains with a RMSD cut-off of 0.15 nm, as described previously [13]. Root-mean-square fluctuations (RMSF) for each amino acid residue were computed using the last 50 ns of each simulation in order to assess whether the different membrane environments change the mobility of the various protein domains.

MD simulation parameters

Gromacs 4.6.7 [86] was used for all MD simulations with the Slipids force field [41] for lipids and the Amber99SB force field [50] for protein and ions. Water molecules were represented using the TIP3P water model [43]. Initial velocities were randomly assigned from Maxwellian distributions. Periodic boundary conditions were applied throughout all simulations. All NVT (isochoric-isothermal) ensemble simulations used the velocity-rescale thermostat [8]. All other simulations were carried out in the NPT (isobaric-isothermal) ensemble at 1 atm pressure and 298 K using the Nosé-Hoover [37], [61] thermostat and the Parrinello–Rahman barostat [66] using a semi-isotropic pressure coupling scheme and the isothermal compressibility of water (4.5 × 10−5 bar−1). Long-range electrostatic interactions were treated using the Particle-Mesh Ewald (PME) scheme [12] with a real-space cut-off of 1.4 nm and a Fourier spacing grid of 0.12 nm. Van der Waals interactions were described using a Lennard-Jones potential with a cut-off of 1.4 nm. Long-range dispersion corrections were applied to the pressure and potential energy. The LINCS algorithm [36] was used to constrain all covalent bonds. A time step of 2 fs and 1 fs was used for the membrane simulations and the protein-membrane simulations, respectively.

Molecular docking simulations

Molecular docking simulations were performed using the dominant conformation of ABCB5 calculated using the final 50 ns of the production runs of the ABCB5-POPC-POPE-cholesterol protein-membrane system. This dominant conformation was identified using clustering analysis on 500 structures from the last 50 ns of the corresponding simulation, using all atoms in the TM domains with a RMSD cut-off of 0.15 nm, as described previously [13]. The ligand structures for rhodamine 123, doxorubicin, daunorubicin, paclitaxel, docetaxel, mitoxantrone, etoposide, camptothecin, vincristine and 5-flourouracil (5-FU) were downloaded from the ZINC15 database [81] (Fig. S5 in the Supplementary Information), energy-minimised using the steepest descents algorithm in Avogadro v.1.2.0 [35] using the MMFF9s force field, and exported as PDB files. PDBQT files with appropriate atomic charges were generated for all ligands using AutoDock Tools V.1.5.6 [60], which were then used for molecular docking simulations. The default protonation states at physiological pH may be incorrect for some ligands, such as doxorubicin and rhodamine-123. There is evidence that a different protonation state may be prevalent at physiological pH for these compounds [25]. As the relevant protonated groups in those molecules already form H-bonds, addition of the charge would only strengthen that interaction. As the binding site location within the TMD region is unknown, the grid box was set for the entire transmembrane domain. A 40 × 40 × 40 Å3 grid, centred at the centroid of the internal cavity with a grid spacing of 1.0 Å was used. The exhaustiveness parameter was manually set to 50 due to the large search space volume required, consistent with a prior docking simulation study of Pgp1 [22], [19]. Docking simulations were performed using AutoDock Vina v1.1.2 [84]. For each compound, the ligand–protein contacts were evaluated through the analysis of all top-ranked docking modes using BIOVIA Discovery Studio Visualiser V2.5.5 (Dassault Systèmes) and PyMOL [77]. The number of residues and the nature of interactions (non-bonded and hydrogen bonds) were assessed.

Results and discussion

Model of full-length ABCB5 transporter

As expected, all models comprised of two TMDs and two NBDs with an open-inward facing conformation (Fig. 2a). The two halves of the transporter, TMD1-NBD1 of the N-terminus (coloured in yellow) and TMD2-NBD2 of the C-terminus (coloured in blue) are interconnected by a linker region (coloured in green) (Fig. 2a). The secondary structure content of the ABCB5 models includes eight β-sheets in the NBDs, and twelve α-helices (six in each TMD), similar to the topology predictions reported previously using HMMTOP [59]. The core of each NBD consists of four parallel and four anti-parallel β-sheets (Fig. 2a, close up view in the box), comprising the ABC signature (LSGGQ; highlighted in black), Walker A (GSSGCGKS; highlighted in blue) and Walker B (highlighted in purple) motifs. The ABC signature motif is a small, α-helical fragment of 5 amino acids. The Walker A motifs in the ABCB5 model exhibited a coiled loop and a short α-helix (P-loop) rich in glycine, while the Walker B motifs (ILLLD) have a β-sheet conformation. These three motifs were conserved between Pgp1 and ABCB5 in NBD2 but showed slight variation in NBD1 of ABCB5, where the Walker A motif sequence (GLNGSGKS) varies by two amino acids (Leu422 and Asn423). The ABC signature and Walker B motif sequences in NBD1 of ABCB5 are found to be MSGGQ and ILLID, respectively (Fig. S1), and vary by one amino acid, namely Met525 and Ile555, respectively. Cross-linking of the ABC signature motif and Walker A motif site is essential for coupling drug binding with ATP hydrolysis [52]. The impact of these sequence variations in these crucial motifs on the transport cycle of ABCB5 is unknown.

Fig. 2

Structure of the human full-length ABCB5 model. (a) Representative ABCB5 model generated by SWISS-MODEL. The TMD1 and NBD1 are coloured in yellow, TMD2 and NBD2 are coloured in blue. The linker region is coloured in green. The α-helices are shown as cartoons. A close-up view of NBD1 in the box shows the secondary structure of NBD1, comprising of α-helices (golden yellow), β-strands and sheets (red) interconnected by coils (green). The Walker A, Walker B and ABC signature motifs are coloured in blue, purple and black, respectively. (b) Portal gates 1 and 2 in the ABCB5 model. The model (shown with a pale blue surface view) consists of two potential portal gates: portal gate 1 formed by TMHs 4 and 6 (left), and portal gate 2 formed by TMHs 10 and 12 (right). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The TMDs consists of twelve α-helices that span in and out of the membrane whilst remaining interconnected by coupling helices (CH) or extracellular loops (ECL). Each TMD1 consists of six TMHs, 1–6 (shown in yellow), whilst TMD2 is made of TMHs 7–12 (shown in blue) (Fig. 2a). A central cavity containing the putative drug binding sites is enclosed within the twelve TMHs. Two portal gates that allow the entry of the drug from the cell membrane are present in the ABCB5 model, with similar architecture as in Pgp1 [49]. Portal gate 1 is formed by TMHs 4 and 6, whilst portal gate 2 is formed by TMHs 10 and 12 and is only partially open (Fig. 2b).

The linker region that connects the two halves of the ABC transporter was previously reported to play an important role in the function and stability of human Pgp1 transporters [22], [24], [38], [76]. The models generated by PHYRE2 and MODELLER either have a missing linker or lack an intact linker. All templates used for modelling of ABCB5 (including PDB structures 3G5U, 3G60, 4F4C and 4Q9H) lack an intact linker. The missing linker region in the template was predicted by SWISS MODEL and I-TASSER to have a random coil conformation by ab initio modelling. Deletion of the 34 residues from the linker region of human Pgp1 did not change the membrane trafficking of the protein but led to the loss of catalytic activity [76]. Shortening of linker region in murine Pgp1 by deleting the 34 residues corresponding to the 34 residues removed in human Pgp1 showed loss of ATPase activity due to the binding of molecules only at NBD1 [18]. Therefore, the linker region is considered to have a prominent role in the transporter activity than just connecting the two halves [24].

The stereochemical quality of all models was estimated by Ramachandran plots, QMEAN, ERRAT and Verify 3D. All models had an acceptable global quality score, except models generated by Phyre2 intensive mode and I-TASSER. For models generated by Phyre2 intensive mode and I-TASSER, the percentage of residues in the favourable regions of the Ramachandran plot was <90% (Fig. S1 in the Supplementary Information). The relatively low quality of the models generated by Phyre2 and I-TASSER can be attributed to the template chosen for modelling (PDB ID 3G5U, with a resolution of 3.8 Å). This template corresponds to the crystal structure of mouse Pgp1 previously reported to have structural distortions in the transmembrane helices (TMHs) [49]. A revised crystal structure of mouse Pgp1 (PDB ID 4MIM, with a resolution of 3.8 Å) was later deposited with significant corrections in TMHs 3–5, 8, 9 and 12, coupling helices and intracellular helices (IH) IH1, IH3 and IH4, both elbow helices, both TMD-NBD connectors and minor corrections to all extracellular loops [49]. As the template selection in Phyre2 was automated [45], this was unavoidable. Based on this, five additional models were generated using I-TASSER after ensuring the inclusion of this revised structure (PDB ID 4M1M) and excluding distorted structures (PDB structures 3G5U, 3G60 and 3G61) from template selection for threading. However, this did not noticeably improve the quality of the predicted models (Fig. S1, models 18–22), suggesting that the low quality of these models was not wholly the result of the choice of template. Overall, the model generated by SWISS-MODEL based on the crystal structure of mouse Pgp1 (PDB ID 4Q9H, with a resolution of 3.4 Å) was found to have the highest global score. Hence, this model was chosen for further studies and is hereafter referred to as the ABCB5 model. Loop refinement of the linker region of the ABCB5 model was performed using the loop optimisation function in MODELLER [40].

Upon superimposition of the entire structure of the template (PDB ID 4Q9H) with the ABCB5 model, the TMDs and NBDs aligned with minimal deviations (Fig. S2 in the Supplementary Information). Larger deviations were observed in the extracellular loops and intracellular coils. Large structural changes are known to occur in ABC transporters when transitioning from the resting to active state and vice versa, as evidenced by crystal structures and molecular simulation studies [2], [20], [31], [39], [47], [87]. This could result in high plasticity in the putative binding sites of the ABCB5 transporter, affecting the accessibility of the amino acids residues that may potentially interact with the ligand and thereby, making it challenging to predict the location of binding sites by molecular docking. Hence, MD simulations of the ABCB5 model in the resting-state conformation were conducted to characterise its conformational stability in a model cell membrane environment, as well as to determine if a single structure or an ensemble of conformations is required for subsequent molecular docking simulations aimed at predicting the putative binding sites in the transporter.

Conformational dynamics of the resting-state conformation of full-length ABCB5 in a model cell membrane

The conformational behaviour of the ABCB5 model in a model cell membrane environment was investigated using unrestrained MD simulations. Several studies have reported the influence of the lipid bilayer and the linker region on the conformation of ABC transporters [3], [14], [21], [39], [62], [65]. These studies revealed structural deformations when mouse Pgp1 was placed in a POPC-cholesterol membrane in the presence of 150 mM NaCl. However, the structure was found to be stabilised by the addition of Mg+2, which match physiological conditions and the protonation of the catalytic His residues [62]. Therefore, these same conditions were used in our MD simulations. The protein was simulated in two different model cell membranes to characterise the influence of lipid composition and leaflet lipid content asymmetry on the structural dynamics of ABCB5 in its resting conformation and the conformation of the putative drug binding sites.

The lipid composition of the first membrane was POPC-cholesterol in a 89:11 mol % ratio, which is similar to several previous simulations of Pgp1 [21], [39], [62]. The lipid composition of the second membrane was POPC-POPE-cholesterol in a 72:15:13 mol % ratio, which is based on the lipid composition of melanoma cells [68]. This approach is supported by a recent study that investigated the impact of different membrane systems (POPE, DOPE, DPPE, POPC, DOPC and DPPC bilayers) on Sav1866, a bacterial ABC transporter, and which revealed that the conformational dynamics of Sav1866 was lipid-dependent, with large-scale conformational changes observed in POPE bilayers and, to a lesser extent, in DOPE and DPPE bilayers [39].

The PPM server was used to predict the relative position and rotation of the twelve TMHs of ABCB5 in the membrane systems (Table S2 in the Supplementary Information). As expected, the resulting protein arrangement in the membranes revealed that the N-terminus of the ABCB5 model has an intracellular or cytosolic location, as in Pgp1 [42], [49], [83].

The backbone root-mean-square deviation (RMSD) of the entire protein, NBDs and TMHs were used as a measure of the average deviation with respect to the initial structure of the ABCB5 model to assess the stability of the transporter (Fig. 3a). A stable conformation was reached after ~100 ns in the POPC-cholesterol membrane system and after ~50 ns in the POPC-POPE-cholesterol membrane system. Root-mean square fluctuations (RMSF) (Fig. S4) show that there is no difference in the conformational flexibility of the different protein domains in the two membrane environments. The RMSD data further shows that in both membrane systems, the TMHs appeared to be the more stable domain, with lower structural deviations compared to the NBDs. The RMSD of the NBDs was found to be significantly higher (confidence level 99%) than that of TMHs. Superimposition of the initial structure of the ABCB5 model (shown in grey in Fig. 3b) with the structures obtained at the end of 200 ns of simulation in POPC-cholesterol (left) and POPC-POPE-cholesterol (right) shows a small shift in the entire TMHs towards the internal cavity occurred in both membrane systems, which decreases the volume of the internal cavity. The NBDs also shifted towards each other, decreasing the volume of the ATP binding cavity (see also Movies 1 and 2 in the Supplementary Information). Previous studies also reported a decrease in the volume of the central cavity of Pgp1 (PDB ID 3G5U), with more water molecules being replaced from the centre of the TM pore in the simulated conformation [20], [62]. In particular, the motion of the NBDs to form a more closed internal cavity was observed in simulations of Pgp1 in a membrane environment [62].

Fig. 3

Conformational behaviour of ABCB5 during the molecular dynamics simulations. (a) Evolution of the RMSD of the entire protein, transmembrane helices (TMHs) and nucleotide-binding domains (NBDs) in the POPC–cholesterol membrane (left), and the POPC–POPE-cholesterol membrane (right). (b) The final conformations of ABCB5 at the end of the simulations in the POPC-cholesterol (left, pale pink) and POPC-POPE-cholesterol (right, orange) membrane systems are overlayed on to the initial ABCB5 model (shown in grey). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Based on the time evolution of the RMSD of the protein, we considered the first 150 ns of the simulation as equilibration and only the last 50 ns were used for subsequent analysis. For both protein-membrane systems, 500 frames from these last 50 ns were used for clustering analysis of the putative rhodamine 123 binding residues in human full-length ABCB5. The rhodamine 123 binding amino acid residues in human Pgp1, identified by previous site-directed mutagenesis studies [15], [53], [54], [55], were mapped onto the human ABCB5 sequence (Table S3 in the Supplementary Information) and used for clustering analysis. The mapped residues in human ABCB5 include Phe330, Tyr337, Phe709, Tyr750, Tyr933, Tyr963, Met966, Ala967 and Glu970. RMSD clustering was conducted for the main chain of the TMDs containing these nine putative rhodamine 123 binding residues with a RMSD cut-off of 0.15 nm, which resulted in a single cluster. This indicates that there is limited motion of the main chain of the putative binding region during the simulations. Fig. 4 shows the superimposition of the rhodamine binding residues in the ABCB5 model with their corresponding conformations during the final 50 ns of simulation in the POPC-POPE-cholesterol model. Relatively minor shifts were observed in the orientation of the sidechains of most residues, except Phe709, Tyr963 and Glu970. Based on this finding, we used a central structure of the single cluster as the conformation for docking simulations.

Fig. 4

The registered shift in the sidechain atoms of the putative rhodamine 123 binding residues in the ABCB5 model simulated in a POPC-POPE-cholesterol bilayer. Overlay of the sidechain atoms (shown as sticks; coloured in aqua blue) of the mapped putative rhodamine 123 binding sites in the ABCB5-POPC-POPE-cholesterol system in the final 50 ns. The protein structure is shown in cartoon representation. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Comparison of the Ramachandran plots of the initial ABCB5 model and its structure after 200 ns of simulation in the POPC-cholesterol and POPC-POPE-cholesterol membrane environments revealed a slight reduction in the percentage of residues in outlier (unfavoured) regions from 1.2% to 0.6 and 0.7%, respectively (Table 1). Some residues in the outliner regions in the structures at the end of the simulations were now found in the favourable or allowed regions of the Ramachandran plot. This could be due to the structural reorganisation of the NBDs and may be behind the observed increase in RMSD (Fig. 3). The remaining outlier residues corresponded largely to the linker region (Fig. S3 in the Supplementary Information), except for Leu1036, Ser1051 and Arg1244 of the NBD2 region. The quality of the ABCB5 models in the two membrane environments was very similar as judged by the Ramachandran plots, QMEAN scores and ERRAT plots. However, the VERIFY-3D score was slightly higher for ABCB5 in the POPE-POPC-cholesterol system. A slightly distorted and inward pivot (towards the internal cavity) was observed in the NBD2 of the ABCB5 model at the end of the simulation in the POPC-cholesterol membrane. Previous MD simulations of mouse Pgp1 in a POPC-cholesterol membrane also reported inward pivots in the NBD regions [21], [58], [62]. Interestingly, no such pivots were observed in the NBDs in the structure of ABCB5 at the end of the simulation in the POPC-POPE-cholesterol membrane. Small changes in the NBDs are known to correspond to large deviations in the TMHs and vice versa [52]. Consequently, the ABCB5-POPC-cholesterol membrane system was not considered any further for the molecular docking studies.

Table 1

Quality of the ABCB5 model before and at the end of the MD simulations.

System	Residues in Ramachandran plot (%)			QMEAN	ERRAT	VERIFY-3D (%)
	Favoured	Allowed	Unfavoured
ABCB5 model	94.4	4.4	1.2	−2.96	95.70	48.05
ABCB5-POPC-cholesterol	94.4	5.0	0.6	−4.29	96.69	53.91
ABCB5-POPC-POPE-cholesterol	93.9	5.5	0.7	−4.15	96.26	56.60

Putative drug binding sites of full-length human ABCB5 identified by molecular docking

The entire internal drug-binding pocket of the ABCB5 model was used as the docking simulation cell. Three putative binding sites, sites 1, 2 and 3, were identified in full-length human ABCB5 from the top-ranked docking modes of ten chemotherapeutic drugs, previously identified as substrates of the human full-length ABCB5 transporter (Fig. 5a). The location of these three putative binding sites in human ABCB5 overlapped with the location of three binding pockets in human Pgp1 referred to as the prazosin-, H- and R-sites, respectively [2], [22], [19], [49], [57], [78]. Furthermore, we used site-directed mutagenesis data for rhodamine 123 and the structure of the complex of paclitaxel with human-mouse chimeric Pgp1 (PDB structure 6QEX) as a reference for our molecular docking predictions. Both rhodamine 123 (a fluorescent dye) and paclitaxel (a chemotherapeutic drug) are common substrates of human Pgp1 and ABCB5 [44], [46], [57], [85].

Fig. 5

Putative binding sites in the ABCB5 model predicted by molecular docking. (a) Location of the three putative binding sites: site-1 (blue), site-2 (red) and site-3 (green) in human ABCB5 model (shown in grey cartoon). (b) Top-ranked docking pose for rhodamine 123 (shown in cyan with hetero atoms such as oxygen and nitrogen coloured in red and blue, respectively) in site-1 of the ABCB5 model (shown in golden yellow in cartoon representation). (c) Ligand-protein interactions of rhodamine 123 in human ABCB5. Rhodamine 123 is shown as sticks and coloured cyan (hetero atoms such as oxygen and nitrogen are coloured in red and blue, respectively. Hydrogen atoms are not shown). The interaction residues in ABCB5 are shown in golden yellow, and Pgp1 is shown in grey and (d) Two-dimensional representation of the ligand–protein interaction of the top-docked pose of rhodamine 123. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

All nine, top-ranked binding modes of rhodamine 123 in human ABCB5 were predicted to be located in a cavity closer to the outer leaflet (Fig. 5b) with a free energy of binding (ΔGb) ranging from −8.4 to −7.2 kcal/mol. This suggests that binding of rhodamine 123 in ABCB5 could be restricted to a single binding pocket, as is the case in Pgp1. A total of 17 residues outlining this cavity in human ABCB5 were identified, and the binding site is hereafter referred to as site-1 (Table S3). Site-1 in human ABCB5 involved TMHs 1, 5, 6, 7 and 12, which is consistent with the site-directed mutagenesis data for the rhodamine 123 binding site in human Pgp1 (TMHs 1, 6, 7, 9 and 12) [57], [74]. An overlay of the rhodamine 123 binding site residues in human Pgp1 and human ABCB5 is shown in Fig. 5c and 5d. Next, we compared these 17 putative interacting residues in human ABCB5 with the 9 residues in human Pgp1 identified by biochemical studies (Table S3). Six (four conserved and two semi-conserved) out of the 17 residues matched the site-directed mutagenesis data in Pgp1, namely Phe330, Phe709, Tyr933, Tyr963, Met966 and Glu970, confirming the high structural and sequence similarity in the binding site between the two proteins. The conserved residues in ABCB5, Phe330, Phe709, Tyr933 and Met966, mainly formed hydrophobic interactions or hydrogen bonds with rhodamine 123, whilst the semi-conserved residues Tyr963 and Glu970 formed hydrogen bonds. The residues that were not identified by molecular docking include Tyr337, Tyr750 (semi-conserved) and Ala967 (conserved) (Table S3). The steric clashes in the model and the possible variation in the orientation of the sidechains of the neighbouring residues of these three residues could potentially have hindered potential interactions with the ligand. Seven additional residues were identified by molecular docking as potential interacting residues in human ABCB5: Met65, Tyr295, Met298, Tyr302, Ile333, His334 and His706. Of these seven residues, Met65 and Tyr302 are conserved whilst Tyr295 are semi-conserved. These seven residues formed weak hydrogen bonds and hydrophobic interactions with rhodamine 123. Two (Met65 and Tyr302) out of seven residues are conserved, while five are non-conserved between Pgp1 and ABCB5. The steric clashes in the model and the possible variation in the orientation of the sidechains of neighbouring residues of these seven residues could potentially lead to these other interactions. However, the prediction of false positive interactions by AutoDock Vina is, of course, possible. Thus, despite finding good agreement in the interactions with the residues mapped from site-directed mutagenesis of human Pgp1, further experimental studies are required to validate these predictions.

Previous molecular docking simulations of rhodamine 123 with mouse Pgp1 (PDB ID 3G5U) using AutoDock Vina showed that rhodamine 123 binds to the R-site and prazosin-site with equal affinity (ΔGb = −9.2 kcal/mol) [22], [19]. However, in our docking simulations, rhodamine 123 was predicted to only interact with the prazosin-site (i.e. site-1; ΔGb = −8.4 kcal/mol). This could be due to differences in the identity of interacting residues and the conformation of sidechains between human ABCB5 and mouse Pgp1. Moreover, the ΔGb values of rhodamine 123 in human ABCB5 was less negative than that reported with Pgp1, suggesting a lower affinity of rhodamine 123 for ABCB5. We used the recently reported crystal structure of human-mouse chimeric Pgp1 co-crystallised with paclitaxel to validate our docking predictions of paclitaxel in human ABCB5 [1]. Our docking simulations predicted that paclitaxel interacts with two binding pockets in ABCB5, namely site-1 (ΔGb = −9.3 kcal/mol) and site-2 (−8.2 kcal/mol) (Fig. 5d). In Pgp1, paclitaxel was found to bind to the prazosin-site (i.e. site-1), in line with our findings (Fig. 6a) because the lower ΔGb value predicted for site-1 indicates preferential binding to site-1 over site-2. The interacting residues in ABCB5 were identified and are listed alongside those of Pgp1 in Table S4 in the Supplementary Information. The putative interacting residues in ABCB5 identified by molecular docking (28 residues) were compared with the residues in the crystal structure of human-mouse chimeric Pgp1. Nine out of the 28 interacting residues predicted by the docking simulations match the residues found to interact with paclitaxel in the crystal structure (Fig. 6, Table S4). Five of these nine residues are conserved and showed hydrophobic interactions and hydrogen bonds with paclitaxel. The non-matching residues are not conserved between the two structures and mostly formed hydrogen bonds, with a few hydrophobic interactions. Nevertheless, the binding region for paclitaxel predicted in human ABCB5 (site-1) is in good agreement with the binding sites in human Pgp1, determined by biochemical and X-ray diffraction crystallography data. Finally, a validation docking simulation of paclitaxel with the human-mouse chimeric structure of Pgp1 (PDB ID 6QEX) using similar parameters resulted in nine binding poses of this ligand within the prazosin binding site, all with very similar interactions to those observed in the crystal structure (Fig. S6).

Fig. 6

Putative paclitaxel binding site in the ABCB5 model. (a) Location of the putative paclitaxel binding pocket in human ABCB5 model. The top ranked docked pose of paclitaxel (cyan sticks) is superimposed with the binding mode of paclitaxel observed in the X-ray structure (pink sticks). The ABCB5 model is shown in golden yellow. (b and c) Side chains interacting with paclitaxel in the human ABCB5 model and human-mouse chimeric Pgp1 residues, parallel to the membrane plane. The interacting sidechains in human ABCB5 and human-mouse chimeric Pgp1 are shown in golden yellow and green, respectively. The heteroatoms in the ligands and the residues such as oxygen, nitrogen and sulphur is shown in red, blue and orange, respectively. Hydrogens are not shown. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The binding sites in human ABCB5 for nine other chemotherapeutic drugs considered in this study were predicted (Table S5). Anthracyclines (doxorubicin and daunorubicin) were predicted to bind to two sites, site-1 and site-2 (Fig. 7a). In Pgp1, doxorubicin had been reported to bind to the R-site and prazosin-site by previous biochemical and in silico studies [22], [19], [79], in agreement with our predictions. Three top-ranked binding poses of doxorubicin were predicted at site 2 (ΔGb = −8.3 to −7.9 kcal/mol), and six at site-1 (ΔGb = −7.7 to −7.4 kcal/mol). The ΔGb values of daunorubicin (−10.3 to −9.3 kcal/mol) were lower than doxorubicin (-8.3 to −7.4 kcal/mol), and seven binding poses were found at site-1 (ΔGb = −10.3 to −9.0 kcal/mol) and two at site-2 (−9.3 kcal/mol). In the case of docetaxel, all nine predicted binding poses were found at site-1 (ΔGb = −9.1 kcal/mol). The ΔGb values of taxanes (paclitaxel and docetaxel) at site-1 were lower than those of anthracyclines, indicating a higher affinity of taxanes for site-1. All nine predicted binding poses for docetaxel were found at site-1 (Fig. 7a), whereas paclitaxel was predicted to interact at both site-1 and site-2. Vincristine was also predicted to bind to site-1 and site-2, similar to the anthracyclines. The ΔGb value of vincristine was −8.5 at site-1 and −8.0 kcal/mol at site-2. Overall, the predicted binding poses of anthracyclines and taxanes revealed two putative binding sites, site-1 and site-2, identified for the first time in our study [44].

Fig. 7

Putative binding sites in full-length human ABCB5 for the chemotherapeutic drugs (a) doxorubicin, daunorubicin, docetaxel and vincristine; (b) Etoposide, camptothecin, mitoxantrone and 5-FU. The image represents the top docking mode of each ligand in each drug binding site (if present). The docking pose of each of these ligands is represented as sticks.

Docking simulations of etoposide and camptothecin predicted that these drugs bind to site-1 and site-2, whilst mitoxantrone and 5-FU were predicted to bind to a cavity closer to the inner leaflet towards the N-terminus, hereafter referred to as site-3 (Fig. 7b). The ΔGb values at site-1 and site-2 for the top-ranked poses of camptothecin (−9.3 to −8.6 kcal/mol) and etoposide (−9.6 to −8.9 kcal/mol) were similar. Mitoxantrone was predicted to bind to site-1 as well as site-3, with equal ΔGb (−7.7 kcal/mol). The location of site-3 overlaps with the H-site in human Pgp1. The two binding sites, H- and R-sites in Pgp1, have been reported to act cooperatively with each other [78]. In human ABCB5, the existence of such cooperativeness between site-2 and site-3 is unknown and requires thorough investigation.

Docking simulations of 5-FU predicted that it binds to all three binding sites but with different affinities. Nevertheless, the ΔGb of 5-FU was substantially lower (−5.8 to −4.8 kcal/mol) than for all other molecules considered in this study. Previous in vitro studies revealed that 5-FU does not act as a substrate for human full-length ABCB5 transporter in HEK293 cells overexpressing the ABCB5 transporter [44], in line with our predictions. Contrary to this finding, however, a positive enrichment of ABCB5-expressing cells in colorectal cancer following neoadjuvant treatment with 5-FU indicated that 5-FU is a substrate of ABCB5 [89]. It is unclear whether these contrasting observations for 5-FU could be attributed to the ABCB5β isoform and whether the two isoforms have different substrate specificity. Nonetheless, the docking simulations of the full-length ABCB5 predict that 5-FU has low affinity, which suggests that 5-FU might exhibit weaker interactions with ABCB5.

Characterisation of the putative drug binding sites

Following the identification of the three binding sites in the internal cavity of ABCB5, the characterisation of these binding sites was carried out. The amino acid residues lining the three putative binding sites were identified and listed in Table 2 as well as Table S6 (provided as a .xlsx file) in the Supplementary Information. Only two out of the ten chemotherapeutic drugs considered were found to bind to site-3, so more drugs that bind to this site need to be identified to better characterise the site. The amino acids found in the three sites are plotted as a bar graph (Fig. 8a). The percentages of hydrophobic residues (Ala, Val, Leu, Ile, Phe, Cys and Met) in site-1, site-2 and site-3 were found to be 45, 40 and 33%, respectively (Fig. 8b), possibly attributed not only to the transmembrane location of these drug binding sites but also to the hydrophobic nature of the drugs that bind to ABCB5. The drug-binding sites in Pgp1 have also been reported to be hydrophobic, and most ligand–protein interactions in Pgp1 involve hydrophobic interactions [80]. Table S7 in the Supplementary Information lists all of the amino acid residues in sites 1–3 of human ABCB5 and the corresponding residues in mouse Pgp1 (PDB ID 4Q9H). Only a small proportion of residues are conserved between ABCB5 and Pgp1, with the majority being semi and non-conserved residues, which are thought to contribute to the substrate specificity of ABCB5.

Table 2

Amino acid residues composing the three putative binding sites in ABCB5 identified by molecular docking simulations.

Site-1		Site-2		Site-3
Leu 62	Ser 710	Ser 225	Glu 762	Ile 125
Met 65	Phe 713	Arg 226	Ile 763	Ile 129
Pro 66	Leu 742	Val 228	Thr 765	Ala 132
Leu 69	Cys 746	Ile 229	Met 766	Ala 133
Leu 184	Tyr 750	Leu 231	Arg 769	Thr 136
Gln 187	Gly 810	Thr 232	Ile 800	Arg 140
Asn 188	Gln 818	Glu 235	Asp 801	Asp 172
Met 218	Asn 819	Leu 273	Ala 803	Ser 175
Ala 221	Asn 822	Ala 276	Gln 804	Asp 176
Ala 222	Leu 848	Lys 277	Gln 806	Asp 180
Leu 291	Thr 851	Gly 280	Gly 807	His 345
Tyr 295	Gly 852	Arg 283	Ala 808	Ile 854
Met 298	Met 853	Ser 287	Thr 809	Glu 855
Asn 299	Tyr 922	Lys 288	Ser 811	Thr 856
Tyr 302	Ser 925	Leu 678	Arg 812	Met 859
Ile 333	His 926	Leu 679	Ile 813	Asp 866
His 334	Ile 929	Leu 682	Gly 814	Lys 867
Tyr 337	Tyr 933	Lys 686	Val 815	Gln 907
Cys 338	Phe 955	Val 692	Leu 816	Thr 911
Gly 340	Thr 959	Gln 753	Val 973	Lys 914
Ala 341	Tyr 963	Gly 754	Leu 974	Ile 918
Asn 702	Met 966	Tyr 757	Ala 975
Val 705	Ala 967	Gly 758	Glu 977
His 706	Gly 969	Arg 759	Tyr 978
Pro 707	Glu 970	Gly 761	Lys 980
Phe 709

Fig. 8

Characterisation of the three putative drug binding sites in human ABCB5 predicted by molecular docking. (a) Amino acid residue distribution (in percentage) in the three putative drug binding sites in full-length human ABCB5, (b) Percentage of hydrophobic residues in each of the three putative binding sites.

The putative site-1 in human ABCB5 is composed of ~20% aromatic residues (namely Phe and Tyr), whilst site-2 is composed of only 4% aromatic residues. No aromatic residues were found in site-3. Site-2 was found to have relatively more aliphatic non-polar residues. The percentage of aliphatic non-polar residues (Gly, Ala, Val, Leu, Ile and Pro) in site-1, site-2 and site-3 were 39, 50 and 29%, respectively. The percentage of polar residues (Ser, Thr, Cys, Met, Asn and Gln) in site-1, site-2 and site-3 were 34, 22 and 29%. Negatively charged amino acids (Asp and Glu) are more abundant in site-3 (24%) compared to site-1 (2%) and site-2 (8%). Positively charged amino acids (Lys, Arg and His) are more abundant in site-2 (18%) and site-3 (20%), compared to site-1 (6%). Overall, site-1 and site-2 are mostly composed of aromatic, aliphatic, non-polar and positively charged residues, whilst site-3 is mostly composed of polar and negatively charged residues.

Conclusions

Full-length human ABCB5 has been reported to mediate chemoresistance by drug effluxion, but its drug binding sites have remained unknown. In this study, an atomistic model of the full-length transporter was constructed with the highest quality using the crystal structure of mouse ABCB1 as a template. The structural model of ABCB5 comprises two TMDs and two NBDs interconnected by a linker region. Molecular dynamics simulation of the model was conducted in two membrane systems, namely POPC-cholesterol and POPC-POPE-cholesterol (to mimic melanoma cell membranes) under similar physiological-like conditions. Larger deviations in the loop regions and NBDs with minimal deviations in the TMDs were observed in both membrane systems. The NBDs of the model in the POPC-POPE-cholesterol membrane were stable, whilst the NBDs of the model in the POPC-cholesterol membrane pivoted towards the central cavity. The final configuration of the ABCB5 model in the POPC-POPE-cholesterol membrane system was then used for subsequent molecular docking simulations.

For the first time, three putative drug binding sites with distinct properties in human full-length ABCB5 transporter were predicted by molecular docking simulations. The putative binding sites of rhodamine 123, paclitaxel, doxorubicin and daunorubicin are in agreement with the biochemical and structural data available for these compounds in human Pgp1. All compounds tested in this study were predicted to interact with site-1 (the equivalent of the prazosin-site in Pgp1) of the human full-length ABCB5 transporter. This suggests that site-1 is possibly more accessible than the other two binding sites. The predicted interactions made by various drugs with these putative binding sites should be validated by future biochemical studies of human ABCB5.

CRediT authorship contribution statement

Lokeswari P. Tangella: Conceptualization, Formal analysis, Investigation, Visualization, Visualization, Writing - original draft. Mahreen Arooj: Data curation, Methodology, Writing - review & editing. Evelyne Deplazes: Data curation, Methodology, Writing - review & editing. Elin S. Gray: Project administration, Supervision, Writing - review & editing. Ricardo L. Mancera: Conceptualization, Resources, Supervision, Writing - review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.