Discovery of ABCG2/VEGFR2 dual-target inhibitor with anti-drug resistance activity based on the ATP binding site

: Background: ABCG2 plays a critical role in multi-drug resistance in a variety of cancers. Its activity is influenced by ATP hydrolysis and substrate transport, which are synergistic processes but in two independent sites. However, there was no QSAR study of the ABCG2 inhibitors for the two binding sites respectively. Methods: in current study, a QSAR model of ATP binding site was built using DCG scoring strategy, it was used to optimize the structure of a tyrosine kinase inhibitor to avoid the potential drug resistance. Results: a series of novel dual target compounds with arylamide skeleton were obtained, which showed notable activity in enzyme activity test. Compound 7c and 7i showed significant inhibitory activity against drug-resistant cells. Conclusions: the QSAR of ATP binding site of ABCG2 was first discussed, the dual-target compound 7c proposed a new strategy to reverse drug resistance. kinase, ABCG2 ATPase and cell MDR reversal respectively. Molecular docking studies was performed to figure out the interaction patterns between compound 7c and ABCG2 and VEGFR2 respectively. Consequently, compound 7c exhibits select binding affinity towards the two targets.


Introduction
Drug resistance of tumor is one of the most important factors that lead to chemotherapy failure. It depends on several mechanisms, among which protein expression increasing of drug efflux pumps (ATP Binding Cassette-ABC transporters) has been studied over the years. ABCG2, also known as breast cancer resistance protein (BCRP), is one of the most important members of the ABC transporter family, a group of membrane proteins that play a critical role in MDR in a variety of cancers [1][2][3]. Cancer cells would pump these small molecules out, exhibiting drug resistance after exposure to chemotherapeutics. The chemotherapeutic cocktails composed of ABCG2 inhibitor and antitumor agent such as the tyrosine kinase inhibitors (TKIs) were developed to reduce drug resistance. However, combination regimens suffer from drawbacks depending on pharmacokinetic behavior of each drug, unpredictable drug interactions and enhancement of adverse effects [4,5].
Interestingly, many TKIs which take effect by competing with the catalytic domain of oncogenic tyrosine kinases for ATP binding sites, were found to interact with ABCG2 and inhibit its transport activity [6]. Thus, the development of bifunctional agents that simultaneously inhibit tyrosine kinase and ABCG2 represent an intriguing area of research for anticancer chemotherapy. However, the interaction between ABCG2 and its substrates is a more complicated process compared with ligand-TK interaction, the former implicates conformational changes associated with ATP hydrolysis and ligand recognition [7,8]. Thus, the generation of accurate pharmacophore models for identification of substrates/inhibitors of ABCG2 remains a challenge [9,10].
As the cryo-EM structure of ABCG2 has been reported recently, the transmembrane transport mechanism of ABCG2 is becoming clear gradually [11]. The ATP hydrolysis process and substrate transport process of ABCG2 are synergistic, while the ATP binding site and the substrate binding site work independently [12]. The inhibitors such as Ko143 belonging to substrate analogues increase the expression of BCRP gene as a result of the feedback regulation mechanism, and therefore targeting the ATP binding site seems to be a more efficient way to inhibit ABCG2 [13][14][15]. However, for a long time, the specific inhibitors of ABCG2 were found through structural optimization of the known ABCG2 inhibitors [16,17]. There was no organized QSAR studies of the ABCG2 inhibitors at the ATP binding sites individually, due to the lacking of direct detection of ABCG2 ATP hydrolysis activity. The activity of ABCG2 was indirectly measured by detecting the amount of substrate transported, Therefore, it is difficult to distinguish the binding sites of small molecules with ABCG2 inhibitory activity [18,19].
Normalized Discounted Cumulative Gain (NDCG) is a method to measure the quality of rankings. It is often used to measure the effectiveness of network search engine algorithms or related applications in information retrieval. In our previous work, we used the method to define the target of a natural product and clarify its toxic mechanism [20]. Here we distinguish the binding sites of a series reported compounds by NDCG method, present the QSAR model at the ATP binding sites individually. In order to further test the model and develop the bifunctional agents, it served to screen a database of TKIs that we built earlier [21][22][23]. The hit compound and its derivatives were synthesized, with the purpose of achieving the ABCG2/TK dual-target inhibitors with anti-drug resistance activity.

DCG scoring strategy
The activity of ABCG2 is influenced by ATP hydrolysis and substrate transport, which are synergistic processes but in two independent sites. In order to build a QSAR model only for ATP binding site inhibitors, it is necessary to screen the small molecules binding to the site as a training set. The DCG scoring strategy evaluates the ranking similarity of two data sets. If they share a high similarity in ranking, the calculation result of nDCG scoring will be close to 1. In this study, we wrote a script in python to select the compounds randomly, then compare the IC50 ranking with the docking score ranking of the selected compounds to select the nDCG result closest to 1. Seventeen compounds from recent studies were screened out. The DCG score turned out to be 0.98 and 0.56 for the ATP binding site and the substrate site respectively (Table 1). These results indicated that the transport inhibitory activity of these small molecules is only related to the inhibition of ATP hydrolysis. A structure-activity relationship between ABCG2 inhibitors and the ATP binding site was isolated from the complex mechanism, more significantly, the compensatory increase of ABCG2 expression caused by the inhibition of the substrate binding site by small molecules was avoided.  Substrate  43  1  1  8  41  2  2  4  7  3  4  12  40  4  3  1  42  5  7  7  30  6  15  6  37  7  10  3  29  8  6  2  27  9  5  9  35  10  9  17  23  11  13  13  34  12  16  14  28  13  17  11  32  14  11  15  22  15  12  10  33  16  14  16  10  17  8  5  NDCG  1 0.98 0.56

QSAR model of ATP-binding site
The compounds in the training set were described in Table S1. The QSAR model was established by the molecular descriptors method in Molecular Operation Environment (MOE; Chemical Computing Group Inc., Montreal, QC, Canada), which was reported in our previous work [22]. The 435 molecular descriptors of the training set compounds were calculated by the Calculate Descriptors function. Given that too much noise could lead to over fitting of the model, the molecular descriptors implemented in MOE were ranked depending on the importance that represents the significance of correlation with IC50 and the descriptors ranked top 10 were screened out by the SVM method. After filtering the descriptors with SVM method, 10 of 435 molecular descriptors were selected for the further evaluation and validation of the SVM model, and these molecular descriptors were divided into 3 groups in Table 2. The Grid Search result showed that when the value of the parameter C and g were 8 and 0.01 respectively, the LibSVM model demonstrated the best accuracy. The ten-fold cross-validation was utilized to evaluate the constructed SVM model by using the training set. The evaluation result was shown in Table 3. Among the 17 selected compounds, 16 compounds (TP) were correctly predicted, only one compound (FN) was wrongly predicted. The sensitivity (SE) was 94.1%. As for the 281 non-inhibitors, 271 compounds were correctly (TN) predicted and 10 compounds (FP) were incorrectly predicted. The specificity (SP) value was 96.4%. Of all the compounds in the training set, 287 compounds were correctly predicted, with the Q values of 96.3%, which demonstrated that the constructed SVM model was appropriate for the inhibitor and non-inhibitor compounds in training dataset. After screening the molecular descriptors, the normalized linear model manifested significant relation to descriptors (RMSE = 3.59514 R 2 = 0.804431). The equation below showed the generated QSAR model which was comprised by ten descriptors to predict IC50/SD values of ABCG2 inhibitors.
Pred 50 / SD( 50 ) = 1.20479 − 13.50464 * AM1 Eele / SD (AM1 Eele ) + 1.94560 * MNDO E / SD (MNDO E ) + 3.77245 * MNDO Eele / SD (MNDO Eele ) − 0.89409 * PM3 E / SD (PM3 E ) + 9.00507 * PM3 Eele / SD (PM3E ele ) + 0.16430 * pmi / SD (pmi) − 0.25682 * pmi1 / SD (pmi1) + 0.95719 * pmi2 / SD(pmi2) + 0.30571 * pmiX / SD (pmiX) − 0.49045 * pmiY / SD (pmiY) In terms of the equation above, we could learn that the IC50 mainly depends on AM1, PM3 and MNDO, which were three types of general purpose quantum mechanical molecular models. They calculated parameters of a molecule through quantum chemical descriptors, including dipole moment, ionization potential, electron affinity, HOMO and LUMO energy, electronegativity, electrophilicity [24]. It had been confirmed that these descriptors were closely related to the biological activities of drugs by QSAR modeling. PM3E and MNDOE represented the total energy of the molecule, while AM1Eele, MNDOEele and PM3Eele were used to describe electronic energy. The electronic energy of the molecule had a positive correlation with electronegativity and electrophilicity [24], leading to an influence on the biology effect of compounds. The compounds bearing aromatic rings such as phenyl and pyridine as R1 and R2, like compound 40-44 for instance, owned higher electronegativity and electrophilicity comparing to alicyclic rings such as cyclohexyl and adamantyl rings like compound 22-24 due to sp2 hybridization, which resulted in an increase in electronic energy, revealing relatively higher activities over that of alicyclic rings overall according to the statistics. In terms of the substituents, compounds bearing strong electronegativity groups such as trifluoromethyl, nitro and methoxy groups possessed more exceptional activity in contrast to those with no substitution, for which the reason is the same as how aryl ring did. The two factors above were in consistency with the predicted QSAR model. Pmi, pmi1, pmi2, pmiX and pmiY were on behalf of the principal moment of inertia, which was related to the rotational dynamics of the molecule. It was determined by the atomic mass and the number of atoms, as well as the perpendicular distance between the atoms and the rotating center of the molecule. These three parameters have a positive correlation with Pmi and slightly affected the IC50 value since the coefficients of these descriptors in the equation were far less than the ones related to energy, which put electronic energy in the first place as a predominant factor. The structure-activity relationship was showed in Figure  1a. The virtual screening results based on the two models and the interaction between the hit compound and the key residues of the two proteins.

Virtual screening
We used the QSAR model above to screen our previously established fragment based TKIs database, a hit compound 7c with novel skeleton was screened out with potential VEGFR2 inhibitory activity. It was indicated that ABCG2 and VEGFR2 had potential similarity at ATP binding sites, at least they shared the same inhibitor. The compound 7c was docked into the two targets, attempting to explain why it was screened out. As is showed in Figure 1b, it turned out that the results which were filtered based on the descriptors were in highly consistence with the 3D docking results. Each part of the new structure has its own functions: the arylamide part mainly occupied the ATP binding site of VEGFR2, and the structure was highly aligned with Sunitinib. Besides, the related amino acid residues were also in consistence with Sunitinib ( Figure 2e); the part of triazole acetic acid was related to the ATP binding site of ABCG2. It was worth mentioning that the carboxyl group overlapped with the phosphate group of ATP and formed a metal contact with Mg at the bottom of the pocket (Figure 2d), which undoubtedly stabilized the stability of the ligand in the active site. As a consequence, the new structure satisfied the requirements of the two ATP binding sites.

Biological activity of the compounds
The VEGFR2 inhibitory activity of compound 7a-7i was evaluated by cell-free VEGFR2 kinase assay. The ability of MDR reversal was determined through MTT assay. The vanadate-sensitive ATPase assay was conducted to determine the inhibitory activity against ATPase in ABCG2 overexpressing sf9 cells. The concentration-response curves were showed in Figure 2a, 2b, and 2c, respectively. The IC50 values were summarized in Table 4.

Inhibition on ABCG2 and VEGFR2
Compound 7a-7i were detected for the inhibitory activity on stimulated ATPase by vanadatesensitive ATPase assay performed with ABCG2-overexpressing sf9 membranes. Compound 7c and 7i manifested fairly high activity compared to the other compounds by maintaining a low IC50 value of 346 nM and 685 nM. A drastic decrease in activity had been found in the change of substituent position in terms of the fluoro-substituted compound 7h and 7i, since 7i owned a fairly good activity, while the activity of 7h was not determined. As for other substituents, methoxy-substituted compound 7g bore activity that was not prominent compared to the former ones, and the introduction of hydroxyl on the phenyl ring slightly attenuated the activity, which was indicated by the IC50 value of compound 7b and 7f. Electron-withdrawing groups such as nitro groups also lessen the effect of inhibition for compound 7a had a surge in IC50 value compared to 7d. The concentration-response curves of compound 7c and 7i were showed in Figure 2a. The VEGFR2 inhibitory activity of compound 7a-7i was evaluated by cell-free VEGFR2 kinase assay. Compound 7c, 7f and 7i exhibited relatively significant activity (IC50 values of 3.1 nM, 14.4 nM and 36.2 nM) in contrast with others bearing a value of over 700 nM. Compound 7c, with a chlorine atom substituted on the meta-position of the phenyl ring possessed the most superior activity. The other two compounds, 7f and 7i, containing halogen atoms of bromine and fluorine, also revealed a satisfactory activity of inhibition. For structures bearing no halogen atoms, limited inhibitory activity was observed since compound 7d, 7e and 7g, owning nitro, hydroxyl and methoxy as substituents respectively had extremely high IC50 values of above 10μM. It was found that although the compounds act on ATP binding site, the structure-activity relationship of the two targets is not completely consistent. Typical examples are 7e and 7f, which show great differences in the two enzyme activity experiments. The common feature of these compounds was that there were hydroxyl groups on the benzene ring. We speculated that the hydrophilicity of hydroxyl groups made them different from other compounds in this series when they bound with active sites. The concentration response curves of compound 7c and 7i were showed in Figure  2b.
Compound 7c was selected for further molecular dynamics study. The molecular dynamics simulation results of 7c occupied ATP binding site and substrate binding site were shown in Figure 3. The root mean square deviation (RMSD) value indicated that the transformation of the ligand-protein complex backbone had formed the initial structure. In the results of ATP binding site (red), the RMSD value increased to 2.2Å within 1.5 ns, and then remained between 2.2-2.6 Å until 5 ns, and the average RMSD value was 2.1Å. In the result of substrate binding site (blue), the RMSD value continuously fluctuated within 5 ns and remained unstable. Molecular dynamics results showed that 7c had strong binding ability and stability to the ATP binding site, but poor binding ability and stability to another at ABCG2.

Reversal of Drug Resistance
The activity of compound 7a-7i for MDR reversal was detected by MTT assay using S1-M1-80 cells, which was ABCG2-overexpressing cells (selected by mitoxantrone). Compound 7c and 7i had an overall fair activity with IC50 values of 0.92 μM and 2.8 μM. Compound 7f suggested an acceptable activity of 9.6 μM, while some of the compounds like 7b and 7h were above 20 μM. The rest ones were considered to have little effects on the MDR cells since the value of IC50 is more than 100 μM. The concentration cell viability curves of compound 7c and 7i was plotted by the absorbance values against the logarithmic concentrations of the test compounds. The result was showed in Figure 2c.
All compounds with inhibitory activity against dual targets showed activity against drug-resistant cells, which indicated that the target combination was reasonable. However, not all compounds with inhibitory activity against drug-resistant cells need the activity of both two targets. For example, compound 7h had no inhibitory activity against ABCG2 or VEGFR2, but it showed obvious inhibitory activity against drug-resistant cell. We speculated that compound 7h may inhibit other resistance related proteins such as P-gp (ABCG1), suggesting that the selectivity of these compounds needed to be further improved.

NDCG scoring strategy
The DCG (Discounted Cumulative Gain) algorithm was utilized to examine the consistency between the ranking result of compound IC50 values and that of the scores of docking study. A shell script based on Python 3.8 was established, which had the function to select any combination from the 38 compounds (supporting information). Then the similarity between the ranking of IC50 values and docking scores was compared according to the DCG algorithm. The value of NDCG which was closest to 1 was exported, along with the label of the selected compounds. In terms of these compounds, the docking scores share the highest consistency with the actual activity, which referred to the ATP binding site. The training and testing set of the subsequent QSAR model was established only by the selected compounds.

Compounds Collection and Dataset Construction
38 compounds were collected from recent reports about novel ABCG2 inhibitors [25]. A parallel strategy to Shoichet's and Garcia-Vallve's was applied to develop the decoy sets from the NCI database, PubChem and the Binding DB. First, Tanimoto coefficients between a set of 17 selected compounds were calculated based on an ECFP (Extended-connectivity fingerprints) similarity analysis. The compounds with Tanimoto coefficients less than 0.5 were selected. Second, 5 physical properties including molecular weight, number of hydrogen-bond donors and hydrogen-bond acceptors, number of rotatable bonds, logP were defined by using DecoyFinder 2.0 for the ABCG2 inhibitors collected from the previous step. Then, the compounds with five calculated physical properties similar to any of the ABCG2 inhibitors and with a Tanimoto coefficients higher than 0.5 were further selected in a decoys dataset. Then, 79 decoys from NCI, 134 decoys from Pubchem, and 68 decoys from BindingDB were utilized to construct a training set.

SVM model
First, the 435 descriptors for the compounds sets were calculated using the software MOE 2015 after energy minimization. Then, the 435 descriptors were pre-processed by the data filter plugin in python to remove redundant data: (1) descriptors with large quantities of null, constant and zero values were eliminated, (2) descriptors with very small standard deviation values (stdev <0.5%) were eliminated, and (3) descriptors with high correlation to others (correlation coefficient > 0.95) were eliminated. The remained molecular descriptors were normalized to a range of 0 to +1, which were necessary since the different ranges of their values would influence the quality of the SVM model generated. The SVM method was utilized to perform a further feature selection process to the normalized descriptors. Finally, the LibSVM software was used to select descriptors [26].

QSAR model
The compounds database was established in MOE. All the compounds were prepared using ChemDraw, then the structures were converted to 3D in Chem3D and saved as mol2 formats. Then the mol2 files were imported into database and minimized the energy. Furthermore, an extra field was created to accommodate the IC50 data. The database containing compound structures, IC50 values and the top 10 relatively important descriptors was established. The QSAR model was performed using IC50 as activity field and PCA method. The Regression analysis was performed in QuaSAR module. The RMSE and R 2 values were derived from Quasar fit. The basic parent structure molecule was input as a guide and scaffolds were calculated as presets.

Synthesis and characterization of novel compounds
For all synthesized compounds, 1 H-NMR and 13 C-NMR spectra were recorded on a Varian NMR spectrometer operating at 600 MHz for 1 H, and 150 MHz for 13 C. All chemical shifts were measured in DMSO-d6 as solvents. All chemicals were purchased from Sinoreagent Chemical Reagent (Beijing, China) and were used as received, unless stated otherwise. Analytical TLC is performed on silica gel 60 F254 plates (Qingdao Haiyang Chemical Company) and visualized by UV. Flash column chromatography is performed on gel 60 (40-63 mm) (Qingdao Haiyang Chemical Company). Melting points were determined with an Electro thermal melting point apparatus, are uncorrected.
General procedure for the synthesis of the ethyl 2-azidoacetate (2): To a solution of ethyl 2-bromoacetate (1) (10 g, 59.88 mmol) in 30 mL dichloromethane, sodium azide (4.09 g, 62.87 mmol) in 15 mL water and catalyst sodium iodide (0.18 g, 1.2 mmol) were successively added. The reaction was stirred at 40°C for 8 h. After completion of the reaction, the reaction mixture was added to 30 mL water, followed by extraction with 15 mL dichloromethane for three times. Vacuum removal of organic layer after dying with anhydrous sodium sulfate to give compound 2 as colorless transparent oily liquid, yield: 86.90%.
General procedure for the synthesis of N- (3-ethynylphenyl) benzamide derivatives (5a-5i): A mixture of benzoic acid (3a) (3 g, 24.57 mmol) and acetonitrile (20 mL) was stirred at 0°C for 15 min. Then dichlorosulfoxide (3.51 g, 29.48 mmol) was added dropwise. The reaction was stirred at 40°C for 5 h and monitored by TLC. After completion of the reaction, the solvent and dichlorosulfoxide was removed in vacuum. The residual liquid was dissolved in methylene chloride. The reaction was stirred at 0°C for 10 min and 3-ethynylaniline (4), triethylamine was added. After completion of the reaction, adding water and methylene chloride, the organic phase was separated and dried over anhydrous sodium sulfate to give the compound (5a), as a white solid, yield: 67.88%.
Compound 5b to 5i were prepared with the same method by using corresponding raw materials 3b to 3i and 3-ethynylaniline (4).
General procedure for the synthesis of ethyl 2-(4-(3-benzamidophenyl)-1H-1,2,3-triazol-1-yl) acetate derivatives (6a-6i): A mixture of N-(3-ethynylphenyl) benzamide (5a) (5 g, 22.60 mmol) in 20 mL EtOH was stirred at 40°C for 15 min. The water solution (4 mL) of copper sulfate pentahydrate (0.2 g) and ascorbic acid (0.2 g) was added in and stirred for 5 more min. Then ethyl 2-azidoacetate (2) was added. The reaction was stirred at 40°C for 5 h and monitored by TLC. After completion of the reaction, the solvent was removed in vacuum. Adding water and ethyl acetate, the organic phase was separated and dried over anhydrous sodium sulfate. The solvent was removed in vacuum. The residue was purified by silica gel column to give the pure compound (6a), as a pale yellow oily liquid, yield: 68.71%.
Compound 6b to 6i were prepared with the same method by using corresponding raw materials 5b to 5i.

Cell-free detection of VEGFR2 kinase activity
The inhibitory activity of the compound against VEGFR2 was determined by the ADP-Glo assay. The reaction buffer was prepared by 50 mM HEPES, pH 7.5, 0.1 mg/mL BSA, 10 mM MgCl2, 2 mM DTT and 1% DMSO. The recombinant VEGFR2 kinase and ATP were diluted to 2.2 μg/mL and 250 μM respectively with reaction buffer (10 mM), and the compounds were prepared into solutions of different concentration gradients. 2 μL ATP solution, 1 μL drug solution and 2 μL kinase solution were added to the 96-well plate successively and the reaction was incubated at 37°C for 1h. Then 5 μL ADP-Glo reagent was add to each well and the reaction was incubated at room temperature for 40 min. Finally, 10 μL kinase detection reagent was added into each well and incubated at room temperature for 30 min. Luminescence value was measured by microplate reader.

Inhibition on stimulated ATPase activity
The membranes obtained from ABCG2 overexpressing sf9 cells (containing 2mg total protein) was incubated in ATPase assay buffer (50 mmol/L MES (pH 6.8), 50 mmol/L KCl, 5 mmol/L sodium azide, 2 mmol/L EGTA, 2 mmol/L DTT, 1 mmol/L ouabain, and 10 mmol/L MgCl2; 4mL) and divided equally into 2 portions. Vanadate (100mM Na3VO4, PH 10; 50μL) was added to one of the portions, while the other portion was provided with the same amount of purified water, both of which were transferred into a 96-well plate (40 µL/well; 20 µg protein/well) and incubated at 37°C for 5 min. Then different concentrations of compound 7a-7i and Ko143 were added and incubated at 37°C for 3 min. The ATPase reaction was triggered by adding Mg-ATP (20mM, 10μL/well), reaching a final volume of 50μL/well. After incubation at 37°C for 1h, the reaction was terminated by adding 10% SDS solution (30μL/well) to the reaction mix. The amount of inorganic phosphate release was detected by adding colorimetric reagent which made it feasible to observe the absorbance proportional to the amount of phosphate by microplate reader.

Evaluation of MDR reversal activity by MTT assay
S1-M1-80 cells were seeded into 96-well plates at a density of 2000 cells per well in culture medium (160 µL) and kept under a 5% CO2-atmosphere at 37 °C for 4 h. Different concentrations of the compounds 7a-7i and positive drugs were prepared in culture medium. 20 µL of them were added to the corresponding wells containing S1-M1-80 cells (three rows of wells). Additionally, 20 µL of culture medium were added to wells containing S1-M1-80 as well as S1 wild type cells (each three rows of wells). 20 µL of culture medium (negative control) and 20 µL DMSO (positive control, 1% of concentration) were also added to the corresponding wells. Each well contains a final volume of 200 µL. After an incubation period of 72 h, a solution of MTT (20µL) was added to each well, and the reaction was incubated for 1h. After incubation, MTT was reduced to a water insoluble formazan. The supernatants were removed and the cells are lysed with 100 µL DMSO per well. The color intensity of the formed formazan was determined by measuring absorbance at 540 nm. Concentration cell viability curves were generated by plotting the absorbance values against the logarithmic concentrations of the test compounds to calculate IC50.

Molecular docking and dynamics
Molecular docking was performed using the Cryo-EM structure of the ABCG2 E211Q mutant bound to ATP and Magnesium (PDB id: 6HBU) and the crystal structure of VEGFR2 in complex with Sorafinib (PDB id: 4ASD) The PDB file was downloaded from the Protein Data Bank. Structure of the receptor were modified using the Structure Preparation option and the Protonate 3D command was employed to add the missing hydrogens and properly assign the ionization states. The default procedure in the Dock application was used to find the favorable binding configurations of the studied ligands. Placement poses generated by the Triangle Matcher were rescored and filtered by the London dG Scoring method to pick out those exhibiting outstanding hydrophobic, ionic, and hydrogen-bond contacts to the protein. The results of placement were followed by a refinement step. The generated poses were energy minimized using the Amber10: EHT force field. Eventually, the optimized poses were ranked using the GBVI/WSAdG free-energy estimates. Docking poses were observed and interactions with binding pocket residues were analyzed. GROMACS 5.1.4 was used to perform 10 ns molecular dynamics simulation on the obtained compound, the purpose was to insight the analysis of the kinetic properties of the compound and the receptor complex. First, the Amber99SB-ILDN force field in the pdb2gmx module was used to preprocess the receptor protein. The ligand compound was processed using Antechamber and ACPYPE tools to obtain the topological file of the ligand, and the receptor and ligand files are combined into a complex. The solvent system was set to a single-point charge water model, the distance between the complex and the edge of the box was at least 1.0nm, and ions were added to balance the charge of the complex.
In order to ensure the stability of the complex, avoiding steric collisions between ligand molecules and receptors that generates irregular geometric shapes, each system initially adopted the steepest descent algorithm to optimize 5000 steps. Then the system carries out 100ps NVT and NPT system simulation respectively, so that the system reached the equilibrium conditions of p=1.0bar and T=300K. Finally, each system was subjected to molecular dynamics simulation in a 5ns process under the conditions of p=1.0bar and T=300K, with a step size of 2fs. The trace file was recorded every 10ps, and the RMSD of the compounds in these trajectories was calculated to verify the stability of the system. In the entire MD simulation process, the PME (Particle Mesh Ewald) was used to evaluate the long-range electrostatic interaction, and all the bond integrity was constrained by the LINCS (linear constraint solver).

Conclusions
The QSAR model of ABCG2 inhibitors at ATP binding site individually was built using DCG scoring strategy. The model was precise in accuracy, revealing important pharmacophores binding to this site and the results of the structure-activity relationship were slightly different from those in previous studies.
It also indicated that it was inappropriate to study the structure-activity relationship under the circumstance that the inhibitory activity was not partitioned into specific mechanisms depending on the types of targets. A series of dual-target compounds with a phenyl heteroaryl partial structure were screened out by the ABCG2 model and the VEGFR2 model and synthesized. VEGFR2 kinase assay, ATPase assay and MTT assay were carried out for the series of compounds in order to determine the inhibitory activity against VEGFR2 kinase, ABCG2 ATPase and the multi-drug resistance. Compound 7c exhibited the best inhibitory activity among the series of compounds with IC50 values of 3.1 nM, 346 nM and 0.92μM for VEGFR2 kinase, ABCG2 ATPase and cell MDR reversal respectively. Molecular docking studies was performed to figure out the interaction patterns between compound 7c and ABCG2 and VEGFR2 respectively. Consequently, compound 7c exhibits select binding affinity towards the two targets.

Declarations
Availability of data and materials: Availability of data and materials have been submitted along with the manuscript.
Competing interests: The authors declare no conflict of interest.