In silico drug repurposing and lipid bilayer molecular dynamics puzzled out potential breast cancer resistance protein (BCRP/ABCG2) inhibitors

Abstract Multidrug resistance (MDR) is a fundamental reason for the fiasco of carcinoma chemotherapy. A wide variety of anticarcinoma drugs are expelled from neoplasm cells through the ATP-binding cassette (ABC) transporter superfamily, rendering the neoplasm cells resistant to treatment. The ATP-binding cassette transporter G2 (ABCG2, gene symbol BCRP) is an ABC efflux transporter that plays a key function in MDR to antineoplastic therapies. For these reasons, the identification of medicaments as BCRP inhibitors could assist in discovering better curative approaches for breast cancer therapy. Because of the deficiency of prospective BCRP inhibitors, the SuperDRUG2 database was virtually screened for inhibitor activity towards the BCRP transporter using molecular docking computations. The most potent drug candidates were then characterized utilizing molecular dynamics (MD) simulations. Furthermore, molecular mechanics-generalized Born surface area (MM-GBSA) binding affinities of the most potent drug candidates were estimated. Based on the MM-GBSA binding affinities throughout 150 ns MD simulations, three drugs—namely zotarolimus (SD002595), temsirolimus (SD003393), and glecaprevir (SD006009)—revealed greater binding affinities towards BCRP transporter compared to the co-crystallized BWQ ligand with ΔG binding values of –86.6 ± 5.6, –79.5 ± 8.0, –75.8 ± 4.6 and –59.5 ± 4.1 kcal/mol, respectively. The steadiness of these promising drugs bound with BCRP transporter was examined utilizing their structural and energetical analyses throughout a 150 ns MD simulation. To imitate the physiological environment, 150 ns MD simulations for the identified drugs bound with BCRP transporter were conducted in the 1-palmitoyl-2-oleoyl-phosphatidylcholine lipid bilayer. These findings identify zotarolimus, temsirolimus and glecaprevir as auspicious anti-MDR drug leads that warrant further experimental assays. Communicated by Ramaswamy H. Sarma


Introduction
Multidrug resistance (MDR) is the main impediment responsible for the cancer chemotherapy fiasco (Szakacs et al., 2006;Ullah, 2008). According to current estimates, MDR is responsible for up to 90% of cancer-linked fatalities (Bukowski et al., 2020). Various mechanisms enhance the MDR, like overexpression of ATP-binding cassette (ABC) transporters (Luqmani, 2005;Simon & Schindler, 1994;Wu et al., 2014). ABC transporters play a crucial role in safeguarding and detoxification towards cytotoxic agents by effluxing xenobiotics outside of the cells (Robey et al., 2009). Due to this, they also play a crucial role in drug ADMET (absorption, distribution, metabolism, excretion and toxicity) (Szakacs et al., 2008). Up to now, 52 human ABC transporters, which can be split into seven subfamilies-namely ABCA through ABCG on the basis of sequence conformity-have been recognized (Dean et al., 2001). Breast cancer resistance protein (BCRP) was characterized by Dean et al. as the second member of the ABCG subfamily (ABCG2) (Dean et al., 2001). Indeed, BCRP is a transmembrane protein involving 655 residues with one membrane-spanning domain and one nucleotide-binding domain (NBD). BCRP works as a homodimer on cellular plasma membranes (Robey et al., 2009). BCRP is expressed in several tissues and located on chromosome 7q22, playing several crucial functions, including medication transmission, protection from toxic substances (Allikmets et al., 1998;Mao & Unadkat, 2005) and creating MDR inside neoplasm cells (Fojo & Coley, 2007).
Hitherto, no clinical trial has been successful in combating the BCRP-mediated MDR (Robey et al., 2018). For this reason, the identification of novel BCRP inhibitors could vitalize the bioavailability of the desirable medicines, attenuate MDR and create a more efficacious remedy (Zhang et al., 2016). It is worth noting that a large and growing body of efforts has been executed to develop BCRP inhibitors (Jackson et al., 2018;Toyoda et al., 2019). Several biological mechanisms were proposed to demonstrate the inhibitory activity of propafenone, tariquidar substructures, BMS-599626, sitravatinib, mitoxantrone, PZ-39 and imatinib against the BCRP transporter (Ashar et al., 2020;Ishikawa et al., 2011;Orlando & Liao, 2020;Peng et al., 2009;Yang et al., 2020). At the time that in silico screening of clinical or investigational drugs unveiled that venetoclax, ledipasvir and pibrentasvir may be auspicious inhibitors (calc. k i ¼ 16.4 nm, 8.19 pm, in addition to 1.14 nM, respectively) against BCRP transporter (Ibrahim, Badr et al., 2021). Different chemical databases were also screened against BCRP transporter to explore putative BCRP inhibitors with the assistance of in silico techniques, and eight prospective inhibitors demonstrated binding energy lower than that of the reference inhibitor (BWQ, calc. -55.8 kcal/mol) (Ibrahim, Badr et al., 2022). Besides, numerous edible plant-derived natural products enticed appreciable interest as prospective BCRP drug candidates (Ibrahim, Abdelrahman, Badr et al., in press;Petersen et al., 2021).
In this work, an ongoing effort was sanctified to the identification of novel drug candidates that can potentially inhibit BCRP transporter. SuperDRUG2 database containing more than 4,600 active ingredients of essential marketed drugs were virtually screened towards BCRP transporter. Binding affinities of the best hits complexed with BCRP transporter were evaluated over the 150 ns molecular dynamics (MD) simulations with the assistance of the molecular mechanics-generalized Born surface area (MM-GBSA) approach. The steadiness of the identified drugs complexed with the BCRP transporter was examined utilizing the structural and energetical analyses during the 150 ns MD course. To mimic physiological conditions, the identified potent drugs complexed with BCRP transporter were simulated in a 1-palmitoyl-2-oleoyl-phosphatidylcholine (POPC) membrane. In summary, these results shine a light on the potency of zotarolimus, temsirolimus and glecaprevir as putative therapies for MDR patients with breast cancer and permit further in vitro/in vivo validation.

BCRP preparation
The cryo-electron microscopy resolved three-dimensional (3D) structure of the BCRP transporter [PDB ID: 6FFC, resolution: 3.56 Å (Jackson et al., 2018)] complexed with BWQ/MZ29 was obtained and used as a template for all in silico computations. The BCRP transporter was prepared by building the missing residues using Modeller software (Marti-Renom et al., 2000). Besides, all heteroatoms, including crystallized water molecules, inhibitors and ions were removed. The Hþþ web server was employed to assign the protonation states of the amino acids (Gordon et al., 2005). All missing hydrogen atoms were also inserted. The pK a for the BCRP amino acids were inspected under the physiologic conditions of 6.5, 80, 10 and 0.15 for pH, external dielectric, internal dielectric and salinity, respectively (Li et al., 2013).

SuperDRUG2 database preparation
To carry out virtual screening (VS), the SuperDRUG2 database was downloaded and prepared (Siramshetty et al., 2018). All compounds were retrieved in two-dimensional (2D) structures stored in SDF format. Thereafter, the 3D chemical structures were generated using Omega2 software (Hawkins et al., 2010;"OMEGA 2.5.1.4," 2013). The energy window of 10 kcal/mol was employed to create all conformers. The lowest energy conformer was subsequently optimized using the MMFF94S force field within SZYBKI software (Halgren, 1999;"SZYBKI 1.9.0.3," 2016). FixpK a and tautomer from the QUACPAC program were applied to scrutinize the ionization states and enumerated the tautomeric forms of the investigated compounds ("QUACPAC," 2016). The Gasteiger-Marsili method was employed to evaluate the charges of each drug inside the SuperDRUG2 database (Gasteiger & Marsili, 1980). Duplicated compounds were removed according to molecular international chemical identifier keys (InChIKey) from the SuperDRUG2 database (Heller et al., 2015). Prepared SuperDRUG2 data files can be obtained from www.compchem.net/ccdb. A workflow diagram of the applied in silico approaches for the screening process of the SuperDRUG2 database is illustrated in Figure 1.

Molecular docking
All molecular docking estimations were conducted with the assistance of AutoDock4.2.6 software (Morris et al., 2009). Prior to docking calculations, relaxation of the BCRP transporter was performed by AMBER16 software (Case et al., 2016) using AMBER force field 14SB (Maier et al., 2015). According to the AutoDock protocol (Forli et al., 2016), the pdbqt file for the BCRP transporter was created utilizing the MGL tools 1.5.7 (www.mgltools.scripps.edu). THR435 and ASN436, key amino acids within the binding side of chains A and B, were placed as flexible amino acids; however, all other amino acids were handled as rigid portions. The Lamarckian genetic algorithm (LGA) was utilized to execute the conformational search (Morris et al., 1998). For all molecular docking computations, the number of generations in each run was adjusted to 27,000 with a population size of 300 individuals. In the current study, fast and expensive docking calculations were performed with number of genetic algorithm (GA) runs of 75 and 250, respectively. In addition, the maximum number of energy evaluations (eval) was 7,500,000, and 25,000,000 for fast and expensive docking computations, respectively. The default values for the other docking parameters were kept in both fast and expensive docking computations. The grid box size was positioned at 80 Å � 80 Å � 80 Å in the x, y, and z dimensions, which is able to accommodate the whole binding pocket of the BCRP transporter. The grid maps with a spacing value of 0.375 Å were created with the assistance of the AutoGrid4.2.6 program (Morris et al., 1998). The docking grid box was centered at 130. 869, 126.675 and 145.206 (in xyz axes, receptively). The predicted docking poses were clustered according to an RMS positional deviation of 1.0 Å (Mansourian et al., 2015). Besides, the pose with the lowest docking score in the greatest cluster was selected as a representative binding pose.

Molecular dynamics simulations
AMBER16 software was applied to execute all MD simulations for the most potent inhibitors bound with the BCRP transporter (Case et al., 2016). General AMBER force field (GAFF2) was adopted for the parameterization of the scrutinized drugs (J. Wang et al., 2004). While, the BCRP transporter was characterized by utilizing AMBER force field 14SB (Maier et al., 2015). The details of the used MD simulations are described by  The investigated drugs were optimized at the HF/6-31G � ab initio level calculation using Gaussian09 software (Frisch et al., 2009). Besides, the charges were computed utilizing the restrained electrostatic potential (RESP) fitting approach (Bayly et al., 1993). As counterions, monovalent ions (Na þ / Cl À ) were inserted to neutralize the unbalanced charges in the investigated systems, and the NaCl salt concentration was maintained at 0.15 M. Thereafter, each investigated complex was immersed into an octahedral TIP3P water box that is protracted 12 Å from any solute atom in all three dimensions (Jorgensen et al., 1983). Initially, all heavy atoms were preserved fixed, and the hydrogen atoms were minimized for 2,500 steps. Then, the entire system was minimized using 5,000 steps. The minimized complexes were thereafter slightly annealed up to 300 K throughout 50 ps utilizing the NVT ensemble with applying a weak restraint of 10 kcal/mol/ Å 2 on the solute atoms. The investigated systems were afterward equilibrated for 10 ns under NPT conditions. Ultimately, each equilibrated complex was submitted to the production stages over 5, 25 and 150 ns, respectively. During the MD simulations, the long-range electrostatic conditions were handled with the help of the Particle Mesh Ewald (PME) method (Darden et al., 1993). The time step of MD simulations was set to 0.002 ps to integrate the equation of motion. As well, a 12 Å cut-off was employed for nonbonded interactions. The Langevin dynamics with the gamma_ln of 1.0 ps À 1 for collision frequency was employed to keep the temperature at 298 K. The Berendsen barostat was employed to control the pressure, with a relaxation time of 2 ps (Berendsen et al., 1984). All bonds consisting of hydrogen atoms were constrained using the SHAKE method (Miyamoto & Kollman, 1992). For binding affinity computations and post-dynamics analyses, the trajectories of the simulated complexes were recorded for each 10 ps. All MD simulations were conducted using the GPU version of PMEMD. BIOVIA Discovery Studio Visualizer 2020 was utilized for all intermolecular interaction visualizations ("Dassault Syst� emes BIOVIA, B.D.S.V.; Version 2019," 2019).

MD simulations in the lipid bilayer
The drug-BCRP complexes were initially incorporated into a POPC membrane bilayer with the assistance of the CHARMM-GUI web interface (Jo et al., 2008). The POPC bilayer was then parameterized utilizing the AMBER Lipid14 force field. Afterwards, the primary axis of the drug-ABCG2 complexes under investigation was placed parallel to the z-axis of the box to be in the center of the POPC membrane. Ultimately, using the identical settings as in the previous section, the MD simulations were conducted.

MM-GBSA binding energy
The MM-GBSA approach was utilized to compute the binding free energy (DG binding ) between the BCRP transporter and the investigated drugs (Massova & Kollman, 2000). GAFF2 and AMBER force field 14SB were utilized to describe drug and BCRP transporter, respectively (Maier et al., 2015;J. Wang et al., 2004). In MM-GBSA, the polar solvation energy (DG pol ) term was calculated by the modified GB model proposed by Onufriev et al. (2004) (igb ¼ 2). The binding free energy (DG binding ) was evaluated on the basis of the uncorrelated snapshots obtained during the MD simulations on the basis of the following equation: where the energy item (G) is estimated as: G SA and G GB denote nonpolar and polar components of desolvation-free energy, respectively. E ele stands for electrostatic energy. E vdw is van der Waals energy. In this work, entropy contribution was ignored because of the high computational cost and time (Hou et al., 2011;E. Wang et al., 2019). The POPC was kept in the MM-GBSA binding energy computations in MD simulations in the lipid bilayer.

Results and discussion
BCRP transporter has a vital role in the expulsion of different endogenous and exogenous substrates involving medicines (Adachi et al., 2005;Ando et al., 2007;Hirano et al., 2005;Jonker et al., 2005;Mizuno et al., 2004;Mizuno et al., 2007). Consequently, the BCRP transporter is known as a substantial factor in the determination of the pharmacokinetic properties of several drugs. In searching for novel molecules to combat MDR, advanced computational approaches were applied to explore the SuperDRUG2 database containing more than 4,600 active ingredients of essential marketed drugs as prospective BCRP inhibitors. For these reasons, a combination of molecular docking and molecular dynamics was utilized, as described in the following sections.

SuperDRUG2 database virtual screening
The outperformance of AutoDock4.2.6 software in exposing the drug-BCRP docking pose and affinity was formerly reported using the same docking parameters (Ibrahim, Badr et al., 2022;Ibrahim, Badr et al., 2021). As documented, AutoDock4.2.6 software predicted the binding mode of the experimental bound MZ29 ligand against BCRP with a RMSD value of 0.25 Å (Ibrahim, Badr et al., 2021). Besides, the employed docking parameters predicted the experimental binding energies with a correlation coefficient (R 2 ) value of �0.94 (Ibrahim, Badr et al., 2021).
In the search for highly potent drug candidates to overcome MDR, AutoDock4.2.6 software was utilized to filtrate the SuperDRUG2 database. The SuperDRUG2 database was initially filtrated against the BCRP transporter with fast docking settings (see Methods section for details). According to the predicted fast docking scores, 131 drug candidates elucidated docking scores lower than the co-crystallized ligand (BWQ ¼ -10.13 kcal/mol). Consequently, those 131 drug candidates were subjected to more delicate docking computations with expensive docking settings. The computed docking scores for the most potent 131 drug candidates are listed in Supplementary Table S1. Thirty-two drug candidates unveiled docking scores lower than the co-crystallized inhibitor (BWQ ¼ À 12.5 kcal/mol). 2D molecular interactions of docking poses within the binding site of the BCRP transporter are illustrated in Supplementary Figure S1. As shown in Supplementary Figure S1, most of the scrutinized drug candidates revealed almost congruous docking poses inside the binding site of the BCRP transporter, exhibiting an essential hydrogen bond with ASN436:A, SER443:A and GLU446:A. The predicted docking scores and 2D chemical structures for those drug candidates are summarized in Table 1.
Zotarolimus (ABT-578/SD002595), an immunosuppressant, disclosed the lowest docking score with a value of -16.3 kcal/mol (Table 1). Investigating the docking pose of zotarolimus with the BCRP transporter demonstrated that the OH group exhibits two hydrogen bonds with the two CO groups of ASN436:A with bond lengths of 2.17 and 2.82 Å, respectively ( Figure 2). Besides, the oxygen atom of the CO group displays a hydrogen bond with the NH 2 group of ASN436:A with a bond length of 2.07 Å ( Figure 2). The OH of the tetrahydro-2H-pyran-2-ol ring contributes to a hydrogen bond with the carboxylate group of GLU446:A with a bond length of 2.06 Å ( Figure 2). Temsirolimus (SD003393), an antineoplastic agent utilized in renal cell cancer (RCC) remediation, revealed the second-lowest docking score with a value of -16.2 kcal/mol (Table 1). Inspecting the docking pose of temsirolimus within the binding site of the BCRP transporter pointed out that the CO and NH 2 groups of ASN436:A exhibit three hydrogen bonds with the two oxygen atoms of OH and CO groups with bond lengths of 2.75, 2.32 and 1.94 Å, respectively ( Figure 2). The OH of the tetrahydro-2H-pyran-2-ol ring participates in a hydrogen bond with the carboxylate group of GLU446:A with a bond length of 2.06 Å ( Figure 2). As well, the OH groups of the two methanol groups demonstrate two hydrogen bonds with the hydroxyl group of THR542:A and CO of PHE439:B with bond lengths of 1.88 and 2.06 Å, respectively ( Figure 2). Glecaprevir (SD006009), a Hepatitis C NS3/4A protease inhibitor, unveiled the third-lowest docking score with a value of -15.9 kcal/mol ( Table 1). As shown in Figure 2, the oxygen atom of sulphur dioxide (SO 2 ) interacts with the NH 2 group of GLN398:A and OH group of SER443:A with bond lengths of 2.21 and 2.31 Å, respectively (Figure 2).
Comparing the BWQ with the investigated drug candidates, it can be seen that BWQ demonstrated a good docking score towards the BCRP transporter with an average value of -12.5 kcal/mol, demonstrating two intrinsic hydrogen bonds with THR435:A (2.14 Å) and ASN436:B (2.16 Å) (Figure 2). A comparison of the investigated drug candidates with BWQ unveiled the great potency of these identified drugs as promising BCRP inhibitors.

Molecular dynamics simulations
The trustworthiness of molecular docking approaches is still questioned, in spite of their significance in anticipating protein-inhibitor binding affinities, because of utilizing unpretentious scoring functions (Guedes et al., 2014). As a result, MD simulations must be executed in respect of credible  revelation of solvent influences, binding energies, conformational pliabilities, the steadiness of the ligand-protein complexes and structural specifics (De Vivo et al., 2016;Kerrigan, 2013). For those reasons, the most potent drug candidates with docking scores less than the co-crystallized inhibitor (BWQ ¼ À 12.5 kcal/mol) were submitted to MD simulations, pursued by binding energy computations. To lessen time and computational costs, short MD simulations of 5 ns were first conducted for the thirty-two drugs in complex with BCRP transporter. The corresponding MM-GBSA binding energies were evaluated and listed in Supplementary Table S2. As shown in Supplementary Table S2, ten drugs manifested binding affinities (DG binding ) less than the co-crystallized BWQ ligand (calc. -60.0 kcal/mol) ( Figure 3). Consequently, these drugs were nominated and subjected to 25 ns MD simulations to fulfill more scrupulous drug-BCRP binding energies. Furthermore, the corresponding binding energies were estimated and depicted in Figure 3. It is worth noting that all scrutinized drugs disclosed less binding energies (DG binding ) than the cocrystallized BWQ ligand (calc. -61.6 kcal/mol) (Figure 3). Notwithstanding, three out of these ten drugs-namely SD002595, SD003393 and SD006009-demonstrated binding energies less than -75.0 kcal/mol. The computed binding affinities for SD002595, SD003393 and SD006009 against the BCRP transporter throughout a 25 ns MD simulation were -85.9, -83.2 and -77.5 kcal/mol, respectively (Figure 3). To achieve more credible binding energies, MD simulations for SD002595-, SD003393-and SD006009-BCRP complexes were prolonged to 150 ns. The corresponding binding affinities were also evaluated and illustrated in Figure 3.
The surpass possibility of SD002595, SD003393 and SD006009 as BCRP inhibitors may be returned to their capacity to demonstrate a variety of hydrogen bonds, van der Waals, p-based and hydrophobic interactions with the key residues within the active site of BCRP transporter. The average structures for SD002595, SD003393 and SD006009 within the BCRP transporter binding pocket throughout 150 ns are shown in Supplementary Figure S2.

Post-dynamics analyses
Energetical and structural analyses were performed over the simulation time of 150 ns and compared to those of the cocrystallized BWQ inhibitor to over and above emphasize the steadiness and stability of SD002595, SD003393, SD006009 and BWQ bound with the BCRP transporter. The structural and energetical stabilization of the identified drugs were attained by investigating hydrogen bond analysis, root-mean-square deviation (RMSD), the binding free energy analysis and center-of-mass (CoM) distance.

Binding free energy analysis
The binding energies of the identified drugs with the BCRP transporter were computed from the entire MD simulation trajectories with the assistance of the MM-GBSA approach. The correlation between MM-GBSA binding free energy and time was measured to evaluate the energetical constancy of SD002595, SD003393, SD006009 and BWQ bound with the BCRP transporter throughout a 150 ns MD simulation ( Figure  4). As shown in Figure 4, the global steadiness of SD002595, SD003393, SD006009 and BWQ was observed with average binding energies (DG binding ) of -86.6 ± 5.6, -79.5 ± 8.0, -75.8 ± 4.6 and -59.5 ± 4.1 kcal/mol, respectively. As a result, it can be concluded that all identified drugs possess a powerful affinity against the BCRP transporter.
The binding affinities were decomposed to identify the dominant forces in the binding of the examined drug  candidates with the BCRP transporter. The MM-GBSA binding energy (DG binding ) and its components are illustrated in Figure 5.
Decomposition of the binding energy unveiled that van der Waals interactions (DE vdw ) were found to dominate the interactions of SD002595, SD003393, SD006009 and BWQ with the BCRP transporter with DE vdw values of -115.2, -107.4, -105.3 and -73.9 kcal/mol, respectively. In addition, the electrostatic attractions (DE ele ) favourably contributed to the total binding energy .1 kcal/mol for SD002595-, SD003393-, SD006009-, and BWQ-BCRP complexes, respectively). Notably, DE vdw is approximately fourfold more vigorous than the DE ele . Together, these findings underpin quantitative data of the binding affinities of SD002595, SD003393, SD006009 and BWQ as potential BCRP inhibitors.
In addition, the binding free energies in the SD002595-, SD003393-, SD006009-, and BWQ-BCRP complexes were decomposed into the participation of residues utilizing the MM-GBSA approach. All amino acids with energetical contributions lower than -0.50 kcal/mol were considered and displayed in Figure 6. It is noted from per residue distribution graph that ASN436:A, SER443:A and GLU446:A are the main contributor to the binding free energy of the scrutinized drugs in the binding pocket of the BCRP transporter ( Figure 6). Particularly, the vitally important participation of the ASN436:A residue in the total binding affinity was noticed with values of -3.6, -2.7, -4.0 and -2.2 kcal/mol for SD002595-, SD003393-, SD006009-and BWQ-BCRP complexes ( Figure 6).

Hydrogen bond analysis
Hydrogen bond analysis was employed to evaluate the steadiness of hydrogen bonding between the identified drugs and BCRP transporter. The analysis was performed for SD002595, SD003393, SD006009 and BWQ bound with the BCRP transporter throughout the 150 ns MD simulations (Table 2). It is apparent from Table 2 that the identified drugs showed steady hydrogen bonds with ASN436:A with the occupancy of 95.6%, 93.3%, 91.1% and 90.1% for SD002595, SD003393, SD006009 and BWQ, respectively. Besides, GLU446:A with SD002595 and SER443:A with BWQ exhibited strong hydrogen bonds with great occupancy ( Table 2).
The number of hydrogen bonds over the simulation time of 150 ns was evaluated and depicted in Figure 7(a). The average number of hydrogen bonds ranged from 2 to 3 for SD002595-, SD003393-, SD006009-, and BWQ-BCRP complexes. The outcomes from the intermolecular hydrogen bonds confirmed the high stability for the SD002595, SD003393, SD006009 and BWQ bound with the BCRP transporter.

Center-of-mass distance
CoM distances between SD002595, SD003393, SD006009 and BWQ and ASN436:A residue were inspected to obtain a deeper comprehension of the immutability of the drug-BCRP complex during the simulation time of 150 ns (Figure 7(b)). What stands out in Figure 7(b) is the noteworthy stabilities for SD002595, SD003393, SD006009 and BWQ inside the binding site of the BCRP transporter with an average value of 8.6, 9.1, 9.1 and 6.5 Å, respectively. This data firmly established the perfect steadiness of the inspected drugs complexed with BCRP transporter throughout a 150 ns MD simulation.

Root-mean-square deviation
To sample the configurational steadiness of the BCRP transporter throughout the MD simulations, the deviation for the backbone of each investigated system was inspected from the starting structures. The RMSD was estimated and presented in Figure 7(c). As shown in Figure 7(c), the global constancy of SD002595-, SD003393-, SD006009-and BWQ-BCRP complexes was monitored with an average RMSD value of 0.29, 0.34, 0.30 and 0.37 nm, respectively. As a consequence, these post-dynamic analyses confirmed that the identified drug candidates do not impact the BCRP overall topology.

MD simulations in the lipid bilayer
The three potential drugs along with BWQ in complex with the BCRP transporter underwent MD simulations in the presence of a POPC bilayer in an effort to simulate the physiological conditions of the human cell. Binding free energies were then estimated for the drug-BCRP complexes employing the MM-GBSA approach throughout the 150 ns MD simulations. The estimated MM-GBSA binding energies proved that the existence of the membrane bilayer had a slight impact on the drug-BCRP complex, leading to a decrease in the total binding affinity. Without the POPC membrane bilayer, the average MM-GBSA binding energies (DG binding ) for SD002595, SD006009, SD003393 and BWQ were -86.6, -79.5, -75.8 and -59.5 kcal/mol, respectively. While the existence of the bilayer demonstrated binding energies (DG binding ) with values of -80.2, -74.3, -71.8 and -53.1 kcal/ mol for SD002595, SD006009, SD003393 and BWQ, respectively ( Figure 8).

Conclusion
Breast cancer resistance protein (BCRP; ABCG2) is an efflux transporter associated with multidrug resistance in carcinoma. Consequently, discovering BCRP inhibitors could assist in struggling with MDR and carcinoma treatment. SuperDRUG2 database was mined in this work to identify prospective BCRP inhibitors using advanced computational techniques. According to the predicted docking scores, MD simulations were executed for the most potent drug candidates, and the corresponding MM-GBSA binding energies were estimated. SD002595, SD003393 and SD006009 manifested promising binding affinity against BCRP transporter with DG binding values of -86.6, -79.5 and -75.8 kcal/mol, respectively. Eventually, the energetical and structural analyses throughout 150 ns emphasized the perfect steadiness for these drugs bound with the BCRP transporter. The lessened effect of lipid bilayer membrane on the drug-BCRP binding affinity was demonstrated. The current findings bring SD002595, SD003393 and SD006009 to light as prospective drug candidates for the sake of BCRP inhibition.

Institutional review board statement
Not applicable.

Informed consent statement
Not applicable.