A network pharmacological approach to reveal the multidrug resistance reversal and associated mechanisms of acetogenins against colorectal cancer

Abstract Multidrug Resistance (MDR) in tumors is caused by the over-expression of ATP Binding Cassette transporter proteins such as Multidrug Resistance Protein 1 and Breast Cancer Resistance Protein 1. This in silico study focuses on identifying a MDR inhibitor among acetogenins (AGEs) of Annona muricata and also aims at predicting colorectal cancer (CRC) core targets of AGEs through a network pharmacological approach. Twenty-four AGEs were initially screened for their ADME properties. Molecular interaction studies were performed with the two proteins MRP1 and BCRP1. As the structure of MRP1 was not available, an inward-facing conformation of MRP1 was modeled. A Protein-protein interaction network was constructed for the correlating targets of CRC. KEGG pathway and Gene Ontology analysis were performed for the predicted CRC targets. We identified four lead AGEs: Muricatocin B, Annonacinone, Annonacin A and Annomuricin E having a higher binding affinity towards MDR proteins. MD simulation studies performed with the three lead AGEs and the MDR proteins showed that MRP1(DBD): Annomuricin E complex was stable throughout the simulation. Our analysis revealed ABCG2, ERBB2, STAT3, AR, SRC and ABCC1 as CRC targets of the lead molecules. The top 10 signaling pathways and functions of correlative CRC targets were also predicted. We conclude that the identified lead molecules might act as competitive inhibitors for reversing MDR in CRC. Additionally, network pharmacological studies established the correlative CRC targets and their mechanisms of action. Further experimental studies are needed to validate our findings. Communicated by Ramaswamy H. Sarma


Introduction
Colorectal cancer (CRC) is the third most commonly diagnosed cancer for which the mortality rates are higher throughout the world (Rao et al., 2017). The estimated cases of CRC are about 1.4 million every year (Siegel et al., 2017). Owing to the complexity of this disease, multidisciplinary treatment approaches including surgery, radiation and chemotherapy are generally recommended for treating tumor cells (Brouwer et al., 2018;Patil et al., 2017;Rady et al., 2017). Chemotherapy is often advised during pre and post-surgery of the polyps. Some common drugs used for treating CRC are 5-fluorouracil (5-FU), irinotecan (Camptosar) and oxaliplatin (Chau & Cunningham, 2009;Xie et al., 2020).
The development of resistance (multidrug resistance) against these drugs leads to the failure of chemotherapy among a majority of the individuals (Krauze et al., 2014).
Multidrug resistance (MDR) is a phenomenon where the tumor exhibits cross-resistance on exposure to one drug. It may be acquired or inherited within the tumor cells (Liu, 2009). The factors influencing the mechanism of multidrug resistance following chemotherapy are increased drug efflux, mutation of cellular targets, inactivation of apoptotic pathways, altered drug targets, drug compartmentalization, and so on Leary et al., 2018). Increased drug efflux by ATP-dependent pumps constitute for a majority of MDR in cancer cells (Robey et al., 2018). Few membrane transporter proteins (P-gp, MRP1, BCRP1) belonging to the ABC transporter family were reported to be highly involved in multidrug resistance in colon carcinoma (Ceballos et al., 2019;Choi, 2005). The substrate (chemotherapeutic drug) binds at the drug-binding domain (DBD) and triggers the attachment of ATP molecules at the nucleotide-binding domain (NBD) of the protein. As a result, ATP is hydrolyzed leading to the release of energy which is highly responsible for the conformational change of ABC proteins from the inward-facing state to the outward-facing state (Pinkett et al., 2007). This conformational change makes the protein efflux a large concentration of drug to an extracellular environment and consequently, the concentration of drug inside the tumor cell gets reduced (Bar-Zeev et al., 2017) leading to MDR in tumor cells. Overexpression of MDR proteins such as MRP1 and BCRP1 of the ABC transporters was highly reported for their increased drug efflux mechanism and inefficient chemotherapy in CRC (Choi & Yu, 2014;Joyce et al., 2015). Our current study focuses on approaches against the drug efflux processes in these proteins.
Multidrug resistance protein 1 (MRP1), a 190 kDa protein is encoded by the ABCC1 gene (Bakos & Homolya, 2007). The protein contains 1531 amino acids, 17 transmembrane helices, 3 MSDs such as MSD0 (TMs 1-5), MSD1 (TMs 6-11), MSD2 (TMs 12-17) and two NBD (Hipfner et al., 1999). The protein's NBD contains the loops and conserved motifs that are required for ATP binding and hydrolysis. Also, walker A, B and signature sequence (C motifs) are present in the NBD (Cole, 2014). A and B walkers are responsible for H-bond formation and stabilization of ADP, respectively (Lu et al., 2015;Ramaen et al., 2006). The main function of this protein is to confer resistance to a broad range of anti-cancer agents such as anthracyclines, vinca alkaloids and epipodophyllotoxins (He et al., 2011;Robey et al., 2011). MRP1 is capable of transporting a variety of organic anions (Conrad et al., 2002) and also plays a major role in the detoxification mechanism (Suzuki & Sugiyama, 1998).
Breast cancer resistance protein (BCRP1) is the second member of the 'G' subfamily of the ABC transporter superfamily (Mao, 2008). The human BCRP1 is designated as ABCG2 (Ni et al., 2010) that contains 655 amino acids (Nakanishi & Ross, 2012) with a molecular weight of approximately 72 kDa (Diop & Hrycyna, 2005). The protein is predicted to have only six transmembrane helices and only one NBD located at the N terminus (Lai, 2013). The human ABCG2 transporter is a reverse half transporter, with the TMD-NBD arrangement reversed (Horsey et al., 2016) which aids dimerization for MDR expression (Ferreira et al., 2017).
MDR in tumors could be controlled by targeting these ABC transporter proteins either by inhibiting the ATP binding site (non-competitive inhibition) or the substrate-binding site (competitive inhibition) (Zloh et al., 2004). Several generations of inhibitors have been identified over the past few decades (Choi, 2005;Choi & Yu, 2014;Dantzic et al., 2018;Hamed et al., 2019;Mohammad et al., 2018;Ye et al., 2019). Moreover, few researchers have also suggested that a combination of inhibitors of ABC proteins along with the conventional chemotherapeutic drugs could be a promising approach for the reduction of MDR (Bukowski et al., 2020;Chen et al., 2016). Hence, identifying a potent inhibitor from natural products and establishment of the biological mechanisms towards MDR reversal is lacking till date.
Annona muricata (soursop) belonging to the family Annonaceae is being widely used as a natural medicine in treating several diseases (Dilipkumar & Agliandeshwari, 2017). The AGEs (bioactive compounds) present in the leaves and stems of graviola have also shown active cytotoxicity against CRC cells (Oberlies et al., 1995;Qazi et al., 2018;Rady et al., 2018). AGEs are a unique class of C-35/C-37 secondary metabolites (a, b-unsaturated c-lactone) derived from longchain (C-32/C34) fatty acids in the polyketide pathway (Agu et al., 2018;Alali et al., 2010;Moghadamtousi et al., 2015). The AGEs were found to bind and block the activity of ubiquinone-linked NADH oxidase in the plasma membranes of cancer cells (Degli Esposti et al., 1994). This mechanism further leads to the inhibition of ATP production and causes the depletion of ATP (Mulia et al., 2015). AGEs are reported to be neurotoxic at increased concentrations of 10 mg/ml (Coria-T ellez et al., 2018;Gavamukulya et al., 2017) and hence it is not recommended for use as a chemotherapeutic molecule. Results from our previous studies have shown that lead compounds from this plant could be effective against blocking the drug efflux pump (ABCB1) and hence could be used in combination therapy against CRC (Jeevitha Priya et al., 2020).
Our current study aimed at using a novel network pharmacology based approach to identify MDR reversal molecules among AGEs against CRC. We have also attempted at employing extensive computational approaches such as homology modeling, molecular docking and molecular dynamic simulation to validate the molecular mechanisms behind this activity. The study is organized under the following phases: (i). Identification of lead AGEs by molecular docking and MM-GBSA. ii) Prediction of potential targets of lead AGEs based on its association with CRC genes through retrieval of necessary information from validated databases. (iii). Establishment of the target's mechanism based on interacting network using network pharmacology approach (iv). Investigation of the roles of the potential targets of lead AGEs through functional enrichment pathway analysis and (v). Comprehensive stability validation of lead AGEs by MD simulation studies. The outcome of this study may provide insights on recognizing selected AGEs as MDR reversal agents in CRC.

Materials and methods
All molecular interaction studies were carried out using various modules of the Schrodinger molecular modeling software suite (Schr€ odinger Release 2019-1, 2019e). Network pharmacology analysis was performed using relevant databases and online tools. A flowchart summarizing the entire methodology was illustrated in Figure 1.

Retrieval of ligands
AGEs that are known to be rich in the leaves of A.muricata were selected for this in silico analysis. The selected AGEs have been shown to possess antitumor activity (Chen et al., 2012;Dilipkumar & Agliandeshwari, 2017;Liaw et al., 2010). Commonly employed drugs like 5-fluorouracil, irinotecan and oxaliplatin used in the treatment of CRC were also included in the study. Four inhibitors were considered as reference molecules. The ligands selected for the current work are listed in Supplementary Table 1. The two-dimensional structures of AGEs, drugs and inhibitors were retrieved from the ChemSpider database (http://www.chemspider.com/) (Ayers, 2012) and DrugBank (https://go.drugbank.com/) (Wishart, 2006).

Drug likeliness of the AGEs
The drug likeliness of the selected AGEs were assessed using an online server, SwissADME (http://swissadme.ch/) (Daina et al., 2017). AGEs were checked for the Lipinski rule of five parameters (RO5) (Lipinski et al., 1997) to differentiate the drug and non-drug-like molecules. The ADME (Absorption Distribution Metabolism Excretion) of the phytocompounds (AGEs) were also analyzed. The compounds that abide by the RO5 were only considered for molecular docking studies.

Homology modeling of MRP1:
The inward-facing structure of MRP1 was modeled using the online tool SWISS-MODEL (Waterhouse et al., 2018) (https:// swissmodel.expasy.org/). The target sequence (UniProt ID -P33527) used for homology modeling of human MRP1 was retrieved from the UniProt database (The UniProt Consortium, 2019). Blast analysis was done for the target sequence against the protein database and the top hit (PDB ID: 5UJ9: Bovine multidrug resistance protein 1) (Johnson & Chen, 2017) was chosen as the template for modeling. Ramachandran plot was generated using an online server PROCHECK (https://www.ebi.ac.uk/thornton-srv/software/ PROCHECK/) (Laskowski et al., 1993) for validating the constructed model. Further, the generated model was given as input in the PDBsum server for verification of secondary structure elements of the modeled protein (https://www.ebi. ac.uk/pdbsum/) (Laskowski et al., 2018). The built model was taken for molecular docking and the interaction of selected ligands at the DBD of MRP1 was analyzed.

Preparation of ligands
The LigPrep module of the Schrodinger suite (Schr€ odinger Release 2019-1, 2019d) was used for the ligand preparation (Dizdaroglu et al., 2020). Missing hydrogen atoms were added to the structure and the chiralities were retained. The water molecules and the counter ions were removed. Charged groups of the ligands were neutralized and the entire structure was ionized. Epik module was specified and the tautomers were generated up to 32 poses per ligand. The structure was finally refined, optimized and energy minimized using the OPLS3e force field. The generated output file was used for the molecular docking studies.

Preparation of protein
The structure of human BCRP1 was retrieved from protein data bank (PDB ID: 6ETI) with inward conformation and used for docking studies (Jackson et al., 2018). The inward facing structure of MRP1 which was modeled as above was used further. The protein structures were prepared using the protein preparation wizard of the Schrodinger suite (Schr€ odinger Release 2019-1, 2019g) (Sastry et al.,2013). The structures were checked for any persisting problems and missing hydrogen atoms were added. The bond orders were refined at the HET groups. The co-crystallized water molecules were removed and the metal ions present in the structures were deleted. Finally, the protein structures were minimized using the OPLS3e force field and the heavy atoms were converged to an RMSD value of 0.30 Å. These minimized protein structures were considered for performing molecular docking analysis.

Grid generation
The receptor grid was generated using the Glide (v8.2) receptor grid generation module of the Schrodinger (Schr€ odinger Release 2019-1, 2019c; Friesner et al., 2006). The Van der Waals radius scaling factor was set as 1.0 and partial charge cut-off was maintained at 0.25 with no constraints. Residue-specific grids were generated for both the NBD and DBD of the MDR proteins and are given in Supplementary  Figure 1 (Johnson & Chen, 2017;Mitomo et al., 2003;Ramaen et al., 2006;Wang et al., 2019). Supplementary Table  2 gives the grid parameters used for molecular docking study.

Docking of ligands
Ligand docking analysis was performed using the Glide (v8.2) module of the Schrodinger suite (Schr€ odinger Release 2019-1, 2019c) (Friesner et al., 2006). The standard Precision (SP) mode of docking was performed within the generated grids of the receptor. The ligands were docked with the NBD and DBD domains of both the MDR proteins (BCRP1 and MRP1). The OPLS3e forcefield was used for calculating the docking scores (Harder et al., 2016) of the ligands. Similarly, the scaling of Van der Waals radii and partial charge cut-off was set as 0.80 and as 0.15, respectively.

MM-GBSA scoring for protein: Ligand complexes
Molecular Mechanics-Generalized Born model and Solvent Accessibility (MM-GBSA) scores were calculated for the top-scoring (docking score) complexes. The Prime module v3.0 of the Schrodinger suite (Schr€ odinger Release 2019-1, 2019f) was used for rescoring (dG bind) of the complexes (Choudhary et al., 2020). The Prime module was computed with OPLS-AA force field and VSGB solvation model . The docked poses were taken as inputs for the energy minimization of the protein-ligand complexes (E complex ), free protein (E protein ) and free ligands (E ligand ). The DG bind (binding free energy) was determined according to the following equation: The complexes with good relative binding energies (highest negative scores) were shortlisted for molecular stability analysis.

PPI network and cluster analysis
The identified overlapping targets from AGEs and CRC were given as input in the STRING database (https://string-db.org/) (v11.0) (Szklarczyk et al., 2019) and a protein-protein interaction network was constructed for those targets. The generated network was then exported for cluster analysis using ClusterONE algorithm of Cytoscape tool (https://cytoscape.org/) (v3.8.2) (Shannon et al., 2003). Few significant parameters such as Degree of Freedom (DOF), average shortest path length, betweenness similarity, clustering coefficient and closeness centrality were determined with a network analyzer using cytoscape (Pan et al., 2019).

Functional and pathway enrichment analysis of the identified CRC targets
Enrichr tool was used for gene ontology enrichment analysis to classify genes based on their molecular functions, biological processes and cellular components. KEGG pathway diagrams were extracted for the core targets using the Enrichr tool (https://maayanlab.cloud/Enrichr/) (Kuleshov et al., 2016) to determine the genes involved in the important cancer pathways. The signaling and the biological pathways that were involved in CRC core targets of AGEs were visualized with the bar diagrams.

Molecular dynamic simulation (MDS)
MDS was carried out using the Desmond package of the Schrodinger suite (Schr€ odinger Release 2019-1, 2019a, 2019b) for the three complexes BCRP1 (DBD): Annonacinone, MRP1(NBD): Annonacin A and MRP1(DBD): Annomuricin E. (Bowers et al., 2006). OPLS3e force field was chosen for protein interactions and the entire system was solvated with a simple point charge water model. The orthorhombic water box was allowed for a 10 Å buffer region between protein atoms and box sides. The system was neutralized with the addition of 17 cl À ions and the radius of the periodic boundary were maintained by enabling the NPT ensemble. The parameters such as temperature (300 K), atmospheric pressure (1.01325 bar) were specified. The simulation process was carried out for 100 ns for the complexes and intervals for stabilization were maintained at 4.8 ps.

Drug likeliness properties
The drug likeliness properties of the selected twenty-four AGEs were assessed by the SwissADME online tool. Lipinski rule of five parameters including molecular weight ( 500Da), logP (Partition coefficient) ( 5), no. of. rotatable bonds ( 10), hydrogen bond donors ( 5) and hydrogen bond acceptors ( 10) were determined for the AGEs. Out of 24 AGEs, solamin violated RO5 and hence was rejected for further analysis. The ADME properties of the selected AGEs were predicted with the pkCSM online tool and listed in Table 1. The water solubility and the percentage of intestinal absorption of the AGEs were predicted to be optimum. Besides, the phytocompounds had a greater HIA (human intestinal absorption). Except for epomuricenin A and epomuricenin B, all other compounds were predicted to be P-gp substrates. The negative values of BBB indicated that all selected AGEs were not readily permeable across BBB. Additionally, the CNS of the compounds were also predicted to be negative. As the cytochrome P450 is responsible for the metabolism of drugs, the predictor assessed the metabolic profile of compounds. The AGEs included in our study were not inhibitors of CYP3A4 but they were CYP3A4 substrates except muricoreacin. All the selected AGEs were represented to possess higher clearance to achieve steadystate concentrations.

Homology modeling of human MRP1
The homology modeling of human MRP1 (hMRP1) was performed using the SWISS-MODEL tool. The template (PDB: 5UJ9 Bovine multidrug resistance protein 1) had a sequence identity of 90.80% and a query coverage of 0.87 with the target sequence (Human MRP1 UniProt ID: P33527). The modeled hMRP1 is shown in Figure 2(A). The QMEAN (À3.27) and GMQE scores (0.54) indicated the high reliability of the built model. The residues involved in the poorly modeled regions of the protein can be noticed from the local quality estimate plot depicted in Figure 2(B). The residues showed a local similarity score of 0.8 representing the acceptable quality of the model. These parameters confirm that the built model is in good agreement with the native structure of the hMRP1. The structural quality of the built model was also assessed using an online structure verification tool PROCHECK. Ramachandran plot shows the stereochemical properties of the modeled hMRP1. The plot revealed that the backbone dihedral angles of modeled hMRP1 were accurate. The percentage of residues located in the most favored region (90.5%), additionally allowed regions (8.4%), generously allowed regions (1%) and disallowed regions (0.2%) can be observed from Figure 3(A). This plot also confirms that the outlier regions contained no residues, suggesting that there is no amino acid with non-favorable dihedral angles. Moreover, the secondary structural elements (SSE) of the generated model was characterized by the PDBsum server and the topology graph is shown in Figure 3(B). Thus, the structural elements of the built model were validated and the modeled hMRP1 was further utilized for docking studies.

Molecular docking analysis and prime MM-GBSA calculation
The docking of ligands (twenty-three AGEs, three drugs and four inhibitors) with the NBD and DBD of ABC transporter proteins (MRP1 and BCRP1) was performed at standard precision mode (SP) using the generated residue-specific grids in the Glide (v8.2) module of the Schrodinger suite (Schr€ odinger Release 2019-1, 2019c). The best ligand protein pose was obtained from the least Glide score and binding energy was calculated by using the Prime module of the Schrodinger suite v3.0 (Schr€ odinger Release 2019-1, 2019f) from MM-GBSA evaluation.

Docking of ligands with BCRP1
The docking scores of the AGEs at the NBD of BCRP1 ranged from À4.5 to À6.7 kcal/mol. Muricatocin B (À6.7 kcal /mol and dG bind of À64.01 kcal/mol), annonacinone (À6.5 kcal/ mol and dG bind À67.54 kcal/mol) exerted a higher binding affinity with the NBD of BCRP1. Similarly, the range of docking scores at DBD of BCRP1 is observed to be À4.4 to 6.2 kcal/mol. Annonacinone (À6.2 kcal/mol and dG bind of À48.38 kcal/mol), gigantetronenin (À6 kcal/mol and dG bind of À58.57 kcal/mol) showed a high affinity towards the DBD of BCRP1. The docking scores of AGEs at both the domains of proteins are represented in frequency distribution graph (Figure 4). The Prime MM-GBSA dG bind scores of top-scoring AGEs are listed in Table 2. The inhibitors selected for the study also exerted similar binding scores as AGEs when docked with the domains (NDB and DBD) of BCRP1 (Table 3).
On the other hand, the binding affinity of the selected chemotherapeutic drugs was found to be lesser (À4.3 to À4.7 kcal/mol) than the other ligands when docked with BCRP1 (NBD and DBD) listed in Table 3.

Docking of ligands with MRP1
The docking scores of AGEs at NBD of MRP1 was around À4.6 to À7.2 kcal/mol. Annonacin A (docking score À7.2 kcal/mol, dG bind À68.93 kcal/mol) and annomuricin B (docking score À7.05 kcal/mol, dG bind À64.64 kcal/mol) showed a greater affinity at the NBD of MRP1. Similarly, the binding scores of AGEs at the DBD of MRP1 ranged between À4.4 to À6.4 kcal/mol. Annomuricin E had the maximum docking score of À6.4 and with a dG bind of À66.35 kcal/mol. Annopentocin A also showed a binding affinity and dG bind of À6 kcal/mol and À49.41 kcal/mol, respectively. Figure 5 represents the docking scores of AGEs at NBD and DBD of MRP1 in frequency distribution graph. MM-GBSA dG bind scores of top-scoring AGEs are listed in Table 4. The docking scores and dG bind of the inhibitors and chemotherapeutic drugs at NBD and DBD of MRP1 is mentioned in Table 5. The binding affinity of the chemotherapeutic drugs at both the domains of MRP1 ranged from À3.9 to À6.7 kcal/mol. The four top-scoring AGEs binding at both NBD and DBD of their corresponding proteins were subjected to network pharmacological analysis.

Network interaction and clustering analysis of the CRC target genes
The four lead AGEs (Muricatocin B, Annonacinone, Annonacin A, Annomuricin E) with favorable docking and MM-GBSA dG bind scores at the NBD and DBD of MRP1 and BCRP1 were taken for network pharmacology analysis. PPI network analysis was performed to predict the core target genes of the four top scoring lead AGEs in CRC. Further, the interaction between the predicted targets and MRP1 and BCRP1 (coded by ABCC1 and ABCG2 genes, respectively) were also explored. A total of 5948 (DisGeNET-2832, CoReCG-  2056, eDGAR-30) CRC targets and AGEs predicted CRC target genes were retrieved from publicly available databases (ChEMBL-244, SwissTargetPrediction-430). 137 overlapping gene targets (CRC genes and AGEs targeting CRC genes) were identified using the Funrich software and the Venn diagram is represented in Figure 6. Duplicate entries were removed. The interaction network was constructed with 137 nodes and 991 edges with a confidence score of 0.4 and enriched p value of < 1.0e À16 by the STRING database. The node proteins interacting with ABCG2 (BCRP1) and ABCC1 (MRP1) were pooled out (Supplementary Figure 3). The PPI network was reconstructed with 17 node proteins and 49 edges with an enriched p value of 1.11e À16 for MDR transport proteins (Figure 7). However, there was no co-expression observed among the interacting genes from the network with ABCC1 and ABCG2 (Supplementary Figure 4). The highly interactive genes (hub genes) were identified from this network based on the number of interactions. The hub genes identified for ABCG2 were PDGFRB, MET, SRC, STAT3 and ERBB2 and FASN was the identified hub gene for ABCC1.

Cluster and network analysis
The constructed interacting network was subjected to clustering analysis with ClusterOne algorithm using Cytoscape (Shannon et al., 2003). This study resulted in three cluster groups (Figure 8). It was observed that the cluster for ABCG2 and ABCC1 interactive genes (second cluster) was observed with 5 nodes (p value <0.05). Network topology parameters such as DOF (which measures maximum number of interactions between interactive nodes), the shortest path length (denotes the path length between two nodes) the clustering coefficient (represents the degree to which our gene targets tend to cluster), betweenness centrality (the number of times a node lies on the shortest path between other nodes) and closeness centrality (scores each node based on their closeness to all other nodes in the network) were determined for the nodes using network analyzer in Cytoscape (Table 6). The DOF ranged between 7 À 13 and the core targets (ERBB2, STAT3, AR, SRC) interacted with ABCG2 and ABCC1 were identified from the cluster analysis.

Functional and pathway enrichment analysis
The biological functions and signaling pathways for the identified core CRC targets of AGEs were subjected to enrichment analysis using the Enrichr online tool. A total of 690 Gene Ontology (GO) annotations consisting of 564 biological functions, 98 molecular functions and 28 cellular components were obtained. GO annotations on AGEs related to CRC included biological functions such as negative regulation of programmed cell death and apoptotic process, phosphotidylinositiol-4,5-bisphosphate 3 kinase activity and protein kinase activity. The top 10 functions were represented in Figure  9(A-C). From KEGG analysis, the mechanism of action of AGEs in CRC were predicted to be enriched in 111 signaling pathways including MicroRNAs in cancer, pathways in cancer, focal adhesion, prostate cancer, proteoglycans in cancer and central carbon metabolism (Figure 9(D)).

Toxicity analysis and P450 site of metabolism prediction
The toxicity profiling of the four AGEs was performed using pkCSM, ToxiM, Pro Tox II and RS-WebPredictor. Table 7 list the various parameters related to the toxicity of compounds. The lead molecules were predicted to be non-mutagenic (Ames toxicity). Moreover, all the four leads were found to be non hepatotoxic. ToxiM score of muricatocin B (0.89) was higher when compared with the other lead molecules indicating a slightly toxic effect. All the four lead molecules were non carcinogenic. Further, the lead compounds were predicted to be highly cytotoxic in nature. The possible sites (primary, secondary and tertiary) of P450 metabolism in AGEs were predicted using RS-WebPredictor 1.0 tool and the results are displayed in Table 8.

Molecular dynamic simulation (MDS)
Out of the four AGEs which were validated to interact with CRC target genes, one of the AGEs muricatocin B was predicted to be minimally toxic at physiological concentration. Hence, MD simulations were performed with the three lead AGEs and their corresponding proteins (BCRP1 and MRP1).
The following three protein: ligand complexes i.BCRP1 (DBD): Annonacinone, ii. MRP1(NBD): Annonacin A, iii. MRP1(DBD): Annomuricin E were simulated and analyzed further. Simulations were performed using the DESMOND package for analyzing the stability of the docked protein-ligand complexes. The entire simulation run was subjected to 100 ns and the most significant parameters such as RMSD, RMSF, SSE, various protein-ligand interactions, SASA, MolSA, rGyr and intraHB were monitored throughout the simulation period.

RMSD of protein: ligand complexes
The stability of the protein:ligand complexes are monitored by the RMSD graphs. BCRP1(DBD): Annonacinone -The protein RMSD was maintained around 2.5 Å and the ligand fit  protein RMSD was maintained at 4.5 Å. The RMSD of ligand fit protein shown in Figure 10(A) was observed to maintain a stable configuration between 10 ns to 78 ns and later experienced a huge fluctuation towards the end of the simulation. Therefore, it was evident that the complex was unstable. MRP1 (NBD): Annonacin A-the RMSD of the protein fit ligand was portrayed in Figure 10(B). It was seen that the complex showed a varied fluctuation at the every ns of the entire simulation. Moreover, the RMSD value of the protein and protein fit ligand was maintained around 2 Å and 9 Å, respectively. The complex was stable between 30 to 40 ns and 85 to 90 ns. The greater RMSD value of ligand fit protein clearly explains the diffusion of the ligand from its binding site towards the end of the simulation. MRP1(DBD): Annomuricin E -The RMSD of the protein and ligand fit protein was maintained at 7.5 Å and 7.2 Å, respectively. The RMSD plot shown in Figure 10C had less fluctuations indicating that the ligand was retained majorly in its binding site throughout the entire simulation of 100 ns. Additionally, the RMSD value (0.3 Å) of the complex was maintained in an acceptable range and equilibrated at the end of the simulation. Out of the three complexes analyzed, MRP1 (DBD): Annomuricin E complex was stable throughout the entire course of the trajectory and hence other significant MD parameters were discussed further for MRP1(DBD): Annomuricin E in detail.

Protein-ligand (MRP1(DBD): annomuricin E) contacts
The interactions between annomuricin E and the DBD of MRP1 are shown in Figure 11. The ligand forms two hydrogen bonds with the residues ASN 1196, ARG 1197 and ASN 1245. It also interacted with PHE 594 and TRP 1246 through hydrophobic interaction. GLU 1204 was highly involved in forming water bridges. No ionic interactions were found between the protein and the ligand. Supplementary Figure 5 portrays the total number of specific contacts of the protein and the ligand throughout the trajectory. A timeline representation of interactions of MRP1 and Annomuricin E in each trajectory frame is shown in Supplementary Figure 5(A). The interactions of ligand atoms with protein that occur more  than 30% of the simulation time in the chosen trajectory (0 to 100 ns) are depicted in Figure 12. The polar residue ASN 1196 at DBD of MRP1 interacted almost 45% of simulation time with (-OH) atom of annomuricin E. Various contacts between the protein: ligand complexes (BCRP1(DBD): Annonacinone, MRP1 (NBD): Annonacin A) can be observed from Supplementary Figure 6. Similarly, GLN 437 of BCRP1 (33%) and THR 660 of MRP1 (30%) interacted with annonacinone and annonacin A, respectively, throughout the entire simulation (Supplementary Figure 7).

Protein (MRP1) RMSF and SSE
The flexibility of the protein was assessed with the RMSF plot (Supplementary Figure 8) generated for the entire run.
The peaks in the RMSF plot indicate the areas of the protein that fluctuate the most during the simulation. The peaks with lesser fluctuations denote the presence of alpha or beta strands for the particular residue. The N and C terminal of the protein fluctuates higher than the rest of the region as they are flexible. The higher fluctuations were observed at residue numbers 98-100 and 680 À 720 indicate the presence of loop regions in the backbone of protein. Further, secondary structural elements of the protein were monitored throughout the simulation period of 100 ns. The alpha helices and beta strands are represented in red and blue color, respectively, in Supplementary Figure 9. Nearly, 49.68 % of helices and 8.02 % of strands (Supplementary Figure 9(A)) were reported as SSE for the protein (MRP1(DBD)). Similarly, Supplementary Figure 9(B,C) summarizes the SSE  (1 Primary metabolic site; 2 Secondary metabolic site; 3 Tertiary metabolic site. The numbers represent the carbon position (C) of the ligands).

Properties of the ligands
The ligand RMSF plot sketched in Supplementary Figure 10 represents the conformational changes of ligand upon protein binding.
The C 1 , C 7 , C 9 , C 14 , C 28 , C 35 , C 39 and C 42 atoms of ligand (annomuricin E) are exposed to the solvent region and, therefore, experiences an internal fluctuation, whereas the remaining atoms are tightly bound to the binding pocket of the protein and hence no deviations are observed during the trajectory (Supplementary Figure 10(A)). Likewise, atoms of Annonacinone and Annonacin A showed a higher fluctuation from C 1 to C 7 upon protein binding

Discussion
ABC transporters play a major role in the multidrug resistance of cancer cells. The reversal of MDR through the inhibition of drug efflux mechanism by natural compounds is reported to be a promising approach (Hamed et al., 2019;Ozben, 2006). Therefore, our present work focused on identifying MDR inhibitors among a group of AGEs by an in silico analysis and validation using a network pharmacology approach. Several generations of inhibitors are being tried to block ABC transporters, and most of these were found to be toxic and non-specific (Tamaki et al., 2011). The conformational change of most of the MDR proteins depends upon the energy released by ATP hydrolysis. This helps in forming an active channel through which the chemotherapeutic drugs are effluxed out of the cells (Pan et al., 2016). MDR inhibitors identified so far have been targeted to block either the ATP hydrolysis (non-competitive inhibition) or the substrate-binding site (competitive inhibition) of the MDR proteins (Jaramillo et al., 2018). Therefore, docking of AGEs was performed with these domains (NBD and DBD) of the ABC transporter proteins (MRP1 and BCRP1).  The docking of AGEs with the substrate-binding site of ABC transporters has been predicted to be more effective when compared with that of the drugs themselves or their reported inhibitors (Figures 4 and 5 and Tables 3 and 5). Also, the AGEs were observed to possess a greater binding affinity towards the DBD than the NBD of the proteins. The compounds muricatocin B (NBD) and annonacinone (DBD) had the highest affinity for BCRP1. In particular, Annonacin A and Annomuricin E bind more efficiently at the NBD and the DBD of MRP1, respectively. These AGEs were found to possess a higher affinity with their corresponding proteins as compared with the chemotherapeutic drugs and inhibitors. Moreover, the hot spot binding residues of the proteins were also explored by the 2D interaction analysis. Residues of MRP1 such as GLU 1204 (3.01 Å) and ASN 1245 (2.27 Å) interacted through Hbond with annomuricin E at C 28 and C 42 atoms , respectively. Annomuricin E interacted hydrophobically with the residues TYR 440 (3.15 Å), THR 439 (3.61 Å), LEU 381 (3.79 Å) and PHE 385 (3.86 Å) of MRP1 (Supplementary Figure 15(A)). The C 5 and C 6 atoms of annonacin A forms a strong H-bond with LYS 684 (2.12 Å), GLY 681 (2.24 Å) and SER 686 (2.86 Å), respectively (Supplementary Figure 15(B)). The hydroxyl groups present at C 5 of annonacinone interacted with the residue THR 435 (3.86 Å) of BCRP1 through H-bond formation. Correspondingly, hydrophobic interactions were also observed with the residues PRO 485 (3.18 Å), PHE 431(3.71 Å), ALA 394 (3.78 Å), PHE 439 (3.94 Å) of BCRP1 when docked with annonacinone ( Supplementary Figure 15(C)). Moreover, the dG bind of the complexes was also significant. Previous studies have affirmed that the elacridar (third-generation inhibitor) inhibited ABC transporters to a greater extent and thereby increased the intracellular concentration of chemotherapeutic drugs in cancer cells. We have reported in our study that the binding scores of the four lead AGEs and elacridar towards the ABC transporters (MRP1 and BCRP1) were similar.
Besides, the ADME properties of the four lead compounds were also observed to be in an acceptable range (Table 1). The lead compounds solubility (log mol/L) values indicated that they are easily soluble in water at 25 C. In addition, all the four AGEs were highly absorbed through the human small intestine (% intestinal absorption). All of these parameters indicate that the lead compounds have a good absorption property. The BBB (log BBB < À1) and the CNS (log PS < À3) of the leads were predicted to not affect the central nervous system. Furthermore, it was also seen that the lead compounds were non-cytochrome P450 inhibitors. The toxicity of the lead compounds given in Table 7 indicated that these natural phytocompounds (AGEs) are safer than chemically synthesized compounds.
Further, the network pharmacological study explored six AGEs targets (ABCC1, ERBB2, STAT3, AR, SRC, ABCG2) for CRC. Researchers have reported that the expression of ABCG2 may be important in the progression and the metastasis of CRC . Besides, studies have suggested that the higher expression of ERBB2 and STAT3 gene was observed to play an important role in developing resistance towards CRC (Lassmann et al., 2006;Pectasides & Bass, 2015;Spitzner et al., 2014). The genes AR, ABCC1 and SRC were found to be involved in controlling cellular proliferation (Huang et al., 2015) and hence considered to be attractive candidates for tailored chemotherapy in CRC (Kunick a & Sou cek, 2014;Jin, 2020). Our study explored the interaction between the ABC transporters (ABCC1 and ABCG2) and their corresponding significant core targets of CRC. It has been observed that only the FASN gene is predicted to interact with ABCC1 with lesser DOF (5) and hence, it was not determined to be a potential CRC target. Moreover, studies have reported that FASN can inhibit cell migration and invasion in HT29 CRC cell lines and could attenuate signaling pathways to decrease CRC . Our network pharmacology analysis leads to the identification of the genes: STAT3, ERBB2, SRC, MET and PDGFRB as potential AGEs targets which were also found to be interlinked with ABCG2. MET activation has been implicated as an important mechanism for the treatment of chemoresistance (Wood et al., 2021). Several studies suggested that MET is associated with the aggressiveness of CRC and hence it may act as a prognostic target (P erez- Vargas et al., 2013). Various experimental studies have validated the role of PDGFRB in cancer cell progression and proposed it as an essential biomarker in characterizing CRC (Sillars-Hardebol et al., 2010;Steller et al., 2013).
Identification of molecular pathways (negative regulation of programmed cell death, apoptosis and microRNA-related cancer pathways, etc.) is a crucial approach towards understanding the mechanisms of drug resistance (Akao et al., 2011;Zhang et al., 2016). We have attempted to study the effect of AGEs on molecular targets of CRC, with an emphasis on the genes involved in drug resistance. The genes ERBB2, STAT3 and SRC identified from our analysis were found to be involved in the negative regulation of programmed cell death and apoptosis. It has been predicted in GO studies, that ERBB2 and STAT3 are involved in phosphatidylinositol-4,5-bisphosphate 3 kinase activity and nuclear transcription factor complex of cellular component, respectively. The nuclear transcription factor has also been reported as a significant target for CRC treatment (Lambert et al., 2018). Furthermore, ABCC1, ERBB2, STAT3 were highly involved in cancer pathways (especially, MicroRNAs in the cancer pathway). The genes reported in our study could probably act through multiple mechanisms to alleviate drug resistance. Additionally, on assessing the stability of the three lead complexes, the MRP1(DBD): Annomuricin E complex was found to be stable than the other complexes such as BCRP1(DBD): Annonacinone, MRP1(NBD): Annonacin A throughout the entire simulation period of 100 ns. It can also be noticed that the residues LEU 381 and PHE 385 at DBD of MRP1 are found to be interacting similarly with both elacridar and annomuricin E (Supplementary Figure 15(D)). Furthermore, the simulation profile and ligand characteristics of annomuricin E with MRP1 (DBD) revealed that the complex was stable during the entire run. Hence, it is proposed that annomuricin E might be an effective inhibitor of MRP1, however, functional studies are needed to confirm our findings.

Conclusion
In conclusion, this systematic computational investigation taken up by our study has shown that the AGEs muricatocin B, annonacinone, annonacin A, annomuricin E are predicted to be capable of reversing MDR caused by ABC transporters. We have identified the core targets of AGEs and their molecular mechanisms, biological functions and pathways. The molecular dynamic studies suggest that the complex MRP1 (DBD): Annomuricin E was stable throughout the entire simulation period. Hence the co-administration of annomuricin E with conventional chemotherapeutic drugs could inhibit MRP1 competitively and possibly evade MDR in tumors. Further, the network pharmacological analysis identified few core targets of AGEs which may help in improving the efficacy of targeted chemotherapy in CRC patients. In addition, this study has explored the interactions between the ABC transporters and few significant genes of CRC which would help in understanding the mechanism of action of AGEs on CRC.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Funding
The author(s) reported there is no funding associated with the work featured in this article.