Soukayna Baammi
Achraf El Allali
Rachid Daoud- 1Chemical and Biochemical Sciences-Green Processing Engineering, Mohammed VI Polytechnic University, Ben Guerir, Morocco
- 2Bioinformatics Laboratory, College of Computing, Mohammed VI Polytechnic University, Ben Guerir, Morocco
More people are being diagnosed with resistant breast cancer, increasing the urgency of developing new effective treatments. Several lines of evidence suggest that blocking the kinase activity of VEGFR-2 reduces angiogenesis and slows tumor growth. In this study, we developed novel VEGFR-2 inhibitors based on the triazolopyrazine template by using comparative molecular field analysis (CoMFA) and molecular similarity indices (CoMSIA) models for 3D-QSAR analysis of 23 triazolopyrazine-based compounds against breast cancer cell lines (MCF -7). Both CoMFA (Q2 = 0.575; R2 = 0.936, Rpred2 = 0.956) and CoMSIA/SE (Q2 = 0.575; R2 = 0.936, Rpred2 = 0.847) results demonstrate the robustness and stability of the constructed model. Six novel compounds with potent inhibitory activity were carefully designed, and screening of ADMET properties revealed their good oral bioavailability and ability to diffuse through various biological barriers. When compared with the most active molecule in the data set and with Foretinib (breast cancer drug), molecular docking revealed that the six designed compounds had strengthened affinity (−8.9 to −10 kcal/mol) to VEGFR-2. Molecular Dynamics Simulations and MMPBSA calculations were applied to the selected compound T01 with the highest predicted inhibitory activity, confirming its stability in the active pocket of VEGFR-2 over 100 ns. The present results provided the basis for the chemical synthesis of new compounds with improved inhibitory properties against the breast cancer cell line (MCF -7).
Introduction
The number of worldwide breast cancer diagnoses and deaths was estimated at 2.3 million and 685,000 successively in 2020 with higher incidence rates in older women age groups (WHO, 2020; Muhammad et al., 2022). Despite the clinical utilization of diverse treatment modalities, the current high mortality rates associated with breast cancer persist, particularly in individuals diagnosed with triple-negative breast cancer (TNBC). As a result, there is currently an intensive research initiative aimed at developing innovative and more effective treatments for breast cancer. This goal can be attained through an improved comprehension of the underlying pathophysiology of breast cancer. It is widely acknowledged that angiogenesis represents one of the fundamental characteristics of cancer, particularly in the context of breast cancer. This process is crucial for endothelial cell proliferation, migration, and survival (Chavakis and Dimmeler, 2002). Stimulation of angiogenesis is one of the hallmarks of tumor growth and malignancy (Ziyad and Iruela-Arispe, 2011) and is significantly influenced by vascular endothelial growth factors (VEGFs) and their receptors, vascular endothelial growth factor receptors (VEGFRs) (Roy and Mitra, 2011). VEGF-A, VEGF-B, VEGF-C, VEGF-D, and placental growth factor (PLGF) contribute to the activation of VEGFRs receptors including VEGFR-1 (Flt-1), VEGFR-2 (Flk-1/kinase domain receptor (KDR), and VEGFR-3 (Flt-4) (Shibuya, 2010; Shibuya, 2011).
Several cancers, including colorectal carcinoma, breast carcinoma, and hepatocellular carcinoma, are associated with overexpression of VEGF (Hicklin and Ellis, 2005), which has been linked to the development of oncogenes, lack of tumor suppressor activity, and fluctuations in glucose or oxygen levels (Wedge et al., 2005). To stimulate proliferation of adjacent endothelial cells, VEGF interacts with one of three tyrosine kinase receptors (VEGFR-1–3), and its overexpression is associated with autophosphorylation of VEGFR-2 in malignancy (Wang et al., 2020). Therefore, several lines of evidence suggest that preventing this pathway by inhibiting VEGFR-2 kinase activity impairs angiogenesis and tumor development (Wang et al., 2018).
Several VEGFR-2 modulators have been discovered, some of which are currently in clinical trials (Aziz et al., 2016). However, these inhibitors fail in clinical trials due to acquired resistance and their side effects, such as receptor redundancy, thromboembolic complications, proteinuria, hemorrhage, anal fistulas, gastrointestinal perforations (GI), posterior reversible encephalopathy syndrome, hand-foot skin reaction, oral mucositis, diarrhea, thyroid dysfunction, bone marrow suppression, emergence of hypoxia-tolerant tumor cells, selection of hypoxia-resistant malignant clones, increase in circulating nontumoral proangiogenic substances, and endothelial cell polymorphisms (Modi et al., 2020).
To develop new compounds with potent inhibitory activity, 3D-QSAR and pharmacophore modeling approaches are commonly used to design new ligand-based drugs based on a series of assays and syntheses of different analogs of the compound’s structure (Nicolaou et al., 2008; Acharya et al., 2011). The 3D-QSAR methods such as CoMFA (Comparative Molecular Field Analysis) and CoMSIA (Comparative Molecular Similarity Indices Analysis) have been used to generate valid and consistent models for the synthesis of different compounds and evaluation of their activities on therapeutic targets (Zhao et al., 2011). The major obstacle in drug discovery is the lack of ADME (absorption, distribution, metabolism, and excretion) features (Vugmeyster et al., 2012). To address this problem, a computational approach for drug design that incorporates ADMET predictions is used to generate new hits during the development process (Wu et al., 2020).
In this study, 3D-QSAR methods are used to model and predict the anti-cancer activities of new Triazolopyrazine analogs. Molecular docking analysis and molecular dynamics simulations show that the proposed molecules exhibit greatly enhanced binding and inhibitory activity against VEGFR-2, the major breast cancer receptor, compared to the leading known breast cancer drug, Foretnib. The present results provided the basis for the synthesis of novel triazolopyrazine analogs as improved drugs against breast cancer.
Materials and methods
Data collection
A database of twenty-three compounds, including triazole-pyrazine molecules, with their inhibitory activity against breast cancer cell lines (MCF-7), was compiled (Liu et al., 2022). This database was divided into two units for the 3D-QSAR study. Nineteen molecules were randomly selected to be used as training sets for the generation of 3D-QSAR models, and the remaining four molecules were used to evaluate the accuracy of the proposed models. This database was selected for the following reasons: (i) The pIC50 values for the compounds ranged from 4.72 to 5.97; (ii) In the field of therapeutic medicine, triazolopyrazine is the active class I pharmacodynamic structure that exhibits good antitumor activity. The IC50 values of the twenty-three 1,2,4-triazole compounds used in this study were converted to pIC50 using the following expression: (pIC50 = −log IC50) and presented in Supplementary Table S1.
Molecular alignment
The molecular modeling studies in this study were performed using the SYBYL-X2.0 package on a Windows 10 64-bit desktop computer (Wang et al., 2018). SKETCH, a module of the SYBYL program, was used to carefully build the 3D scaffolds of the triazolopyrazine derivatives, and the resulting structures were then minimized using the Tripos force field (Waller et al., 1996). The Gasteiger-Hückel atomic partial charges were calculated using the Powell method and a convergence criterion of 0.01 kcal/mol during the minimization phase (Tsai et al., 2010). After structure design and minimization, the distill component of SYBYL-X2.0 performed molecular alignment on the database using molecule 22 (the most active) as a reference (Tabti et al., 2022) (Figure 1).
FIGURE 1. Superposition and alignment of the investigated compounds using molecule 22 as a template.
3D-QSAR modeling
CoMFA and CoMSIA techniques were implemented using the Sybyl X-2.0 program (Vistoli and Pedretti, 2007). The CoMFA method was developed using electrostatic and steric fields as well as the Lennard Jones and Coulomb potentials (Tosco and Mackey, 2017). For the steric and electrostatic energies, we used a sp3 hybridized carbon atom with a Van Der Waals radius of 1.52 and a net charge of + 1.0, using the standard rate of 30 kcal/mol for the power cutoff calculations (Hadni and Elhallaoui, 2020). In addition to steric and electrostatic fields, the CoMSIA method was developed to calculate hydrophobic, hydrogen bonding donor, and acceptor fields. The CoMFA model considers the same factors (Roy et al., 2015).
Partial least square analysis
The 3D QSAR models were generated using the partial least squares (PLS) method, an extension of multiple regression analysis (Vistoli and Pedretti, 2007). For cross-validation, the leave-one-out method (LOO) was first applied (el Mchichi et al., 2022). Here, a single inhibitor was removed from the data set and the derived model was used to make predictions about the activity of that inhibitor. PLS with leave-one-out cross-validation is used to evaluate the accuracy of the model by assessing the optimal number of components (ONC) and the correlation coefficient of cross-validation (Q2) (Consonni et al., 2009). The leave-one-out cross-validation method was used to calculate the coefficient of determination (R2), the F value (F), and the standard error of the estimate (SEE) (Bouamrane et al., 2022). In addition, to further investigate the robustness of the developed models, they were subjected to external validation with a set of four molecules (Hajjo et al., 2010). The following equation is used to determine the coefficient of determination (Rext2) of the test set:
Models are considered acceptable if each of the following criteria is met simultaneously: Q2 > 0.5, R2 > 0.6, Rext2 > 0.6 (Tropsha et al., 2003).
Y-randomization test
The Y-randomization technique was applied to the derived models as an additional check to eliminate the effects of randomness (Rücker et al., 2007). After each iteration, the pIC50 values of the studied compounds are randomly shuffled several times to obtain a new set of Q2 and R2. QSAR models are reliable when their Q2 and R2 values are low, indicating that the excellent calibration result is not due to random correlation (Roy and Mitra, 2011).
Applicability domain (AD)
Since all QSAR models are developed based on a limited number of molecules, there is a specific region of chemical space where the QSAR model can reliably predict new compounds (Roy et al., 2015). This region is referred to as the application domain (AD) (Cherkasov et al., 2014). Therefore, the accurate application of QSAR models requires the calculation of AD (Sahigara et al., 2012). In this work, we applied a method to define AD based on solving the following equation to determine the effects of different leverage values for all compounds
X is the descriptor matrix of the training set, and Xi is the descriptor vector of compound i.
The diagram validates the QSAR model if the calculated leverage value (h) is smaller than the critical value of leverage (h*) (Qin et al., 2017).
P is the number of descriptors and n is the number of compounds.
Molecular docking study
A molecular docking study was performed to predict the molecular interactions between the active site of the VEGFR-2 target protein and the newly designed molecules (Meng et al., 2011). The 3D structure of the VEGFR-2 target protein (PDB code: 4ASD) was obtained from the Protein Data Bank (www.rcsb.org) (Burley et al., 2019). Autodock tools were used to prepare the protein before docking (Hernández-Santoyo et al., 2013). Water molecules and co-crystallized small molecules were removed from the protein structure, the polar hydrogen and Kollmann charges were added to the structure (Shivanika et al., 2020). The grid box spacing was set to 0.375 Å, the center to (−24.611 Å, −0.388 Å, −10.929 Å), and the lattice size to 20 Å × 20 Å × 20 Å where the co-crystallized ligand interacts with the active residues (Baammi et al., 2023a). Nine poses were constructed for each protein-ligand complex based on docking affinity. The Discovery Studio Viewer was used to display and analyze the docking results to find the important interactions between the ligands and the protein binding site (Adeniji et al., 2020; Kesari et al., 2020; Singh et al., 2022).
Redocking
Docking accuracy was determined by comparing the root-mean-square deviation (RMSD) of the heavy atoms between the docked pose and the crystallographic pose of the ligand (sorafenib) (Wen et al., 2019). If the RMSD is below 2 Å, a molecular docking approach is robust (Baammi et al., 2023b).
Molecular dynamics simulation
To determine the structure-function relationship, MD simulations were performed using the GROMACS 2019.3 software program (Abraham et al., 2015) on the docked complexes of compound T01, the active molecule (molecule22), and Foretinib. The CHARMM27 force field was used for the protein (Lindahl et al., 2010), and the topology for the ligands was generated using the Swissparam server (Zoete et al., 2011). Prior to neutralization in the system with counterions, each complex was resolved in a dodecahedron box (1.0 nm) using the TIP3P water model, and counterions (Na+) were added to neutralize the system (Saini et al., 2021). The steepest descent method was used to achieve both the minimum energy and maximum force, with Fmax set at 1000 kJ/mol/nm (Baammi et al., 2022). To equilibrate the system at 300 K and 1 bar, two 100 ps simulations were performed in rapid succession using canonical NVT and isobaric NPT ensembles. Subsequently, 100 ns molecular dynamics simulations were performed for each molecule. The output trajectories were generated, and the data files were analyzed to better understand the behavior of the protein.
Results and discussion
Results of CoMFA and CoMSIA
The main objective of this step is to build powerful CoMFA and CoMSIA models based on the observed and estimated pIC50 values of different models used for training and test sets (Supplementary Table S2) (Sharma et al., 2016). To create the CoMFA model, steric and electrostatic fields were combined. Meanwhile, thirty one different combinations of the five fields: steric, electrostatic, hydrophobic, H-bond donor, and H-bond acceptor were applied to build the CoMSIA models (Supplementary Table S3) (Lu et al., 2011). Table 1 displays the results produced using the CoMFA and CoMSIA models. The Q2 from the cross-validation provides information about the robustness of the CoMFA and CoMSIA models. If Q2 ≥ 0.3, the built model is only significant at the 5% level, while Q2 ≥ 0.5 means that the model is statistically significant (Chu et al., 2020).
TABLE 1. Summary of 3D-QSAR results.
For the CoMFA model, PLS regression analysis yielded a high R2 (0.936), F values of 29.265, a cross-tested coefficient of determination Q2 (0.575), a standard error estimate (SEE) of 0.102, and a small optimal number of components N (6). The CoMFA model was validated using four triazolopyrazine compounds yielding an Rpred2 of 0.956. The steric and electrostatic fields contributed 64% and 36%, respectively, to the CoMFA model, with the electrostatic field being more influential compared to the steric field. For the CoMSIA model, the values of Q2, R2, F, SEE, and N were 0.575, 0.938, 30.437, 0.100, and 6, respectively, and the resulting model exhibited a high level of internal predictability. The same parameters were used to validate the CoMSIA/SE and CoMFA models. The percentage contributions of the steric and electrostatic fields were 64.2% and 35.8%, respectively. The consistency of the CoMSIA models was demonstrated by an external validation correlation coefficient (Rpred2) of 0.845 > 0.6, indicating a high level of accuracy (Chirico and Gramatica, 2011). Furthermore, all these results showed that the steric field is the most important factor in the development of new anticancer drugs. The correlation between observed and predicted activity for the CoMFA and CoMSIA models (Figure 2) exhibits a satisfactory linear correlation for all molecules except molecule 22.
FIGURE 2. Observed versus predicted pIC50 of training and test sets of Triazoloyrazine derivatives inhibitors based on (A) CoMFA and (B) CoMSIA models.
Applicability domain (AD) and Y- randomization test
William’s plot refers to a graph that compares the leverage values and standardized residuals of a particular group of compounds (Beheshti et al., 2016). There are other methods to define the applicability domain models because, without defining the scope, it is not possible to use any of the QSAR models to predict the activities that the new compounds would perform (Roy et al., 2015). Compounds with leverage greater than the threshold are considered potentially disruptive to the performance of the model and are therefore flagged as influential (Moussaoui et al., 2023; Soufi et al., 2023). In this study, the AD of the CoMFA and CoMSIA QSAR models is defined using William’s plots (Figures 3, 4). For the developed CoMFA/ES and CoMSIA/SE models, the warning leverage (h*) was calculated as 0.47. Our result shows that neither the training nor the test sets contain any compounds that are particularly outliers, and all compounds have leverage values lower than the warning h* value (hi < h*). Hence, the model has such good predictive capabilities and robust statistical parameters (Benfenati et al., 2007). In addition, the Y-randomization test was used to examine the robustness of each model. After each iteration, the independent variables of the studied compounds (pIC50) are randomly shuffled (Edraki et al., 2016). The results of 20 random shuffles for the Y-randomization test are shown in Supplementary Table S4 and Supplementary Table S5. After this test, both models have lower values for Q2 and R2 compared to our original models, indicating that the developed models are very robust and reliable and not due to chance.
FIGURE 3. Williams plot for CoMFA/ES model.
FIGURE 4. Williams plot for CoMSIA/ES model.
Contour maps analysis of CoMFA and CoMSIA
The use of 3D contour maps to represent the QSAR equation is an effective way to illustrate the relationships between VEGFR-2 inhibitors and their anti-cancer activities. The 3D contour maps were made using the Stdev* Coeff field type. Due to its high activity, molecule 22 (Figures 5, 6) was selected as a representative compound to analyze the performance of CoMFA and CoMSIA models. Figure 5 shows a contour map of the steric field of CoMFA, with the effect of the steric field on the activity depicted in green and yellow. This steric field accounts for 64% of the total contribution.
FIGURE 5. Electrostatic (A) and steric (B) contour maps of the CoMFA model around molecule 22. When the peripheral part of the molecule is colored green, the activity of the molecule is enhanced by a group with a large connecting space, and when it is colored yellow, the activity of the compound could be reduced by the presence of the group with the large connecting space. The blue colored regions indicate a preference for positive groups that may lead to an increase in anti-VEGFR-2 activity, while the red colored regions indicate a preference for negative groups.
FIGURE 6. Electrostatic (A) and steric (B) contour maps of the CoMSIA model around molecule 22. When the peripheral part of the molecule is colored green, the activity of the molecule is enhanced by a group with a large connecting space, and when it is colored yellow, the activity of the compound could be reduced by the presence of the group with the large connecting space. The blue colored regions indicate a preference for positive groups that may lead to an increase in anti-VEGFR-2 activity, while the red colored regions indicate a preference for negative groups.
Figures 5A shows a green contour dispersed around the R1 substituent, suggesting that adding more groups at the R1 substituent site could increase the activity of the compound. For example, the bioactivity of compound 12 (pIC50 = 5.75) with a methyl ring was significantly higher than that of compound 1 with a hydrogen ring (pIC50 = 4.86). Figures 5B displays the contour map of the electrostatic field of CoMFA. The blue (80%) and red (20%) colors represent the influence of the electrostatic field (36%) on the activity. The compound of an electron-donating group enhances the activity of the compound, as shown by the blue contour around the molecule, while the compound of an electron with drawing group also decreases the activity of the compound, as shown by the red contour. The activity of the molecule increases when electron-withdrawing groups are placed near the X substituent (Supplementary Table S1), which explains the higher activity of compound 9 compared to compound 3 (pIC50 = 5.43 vs. pIC50 = 5.19). The same observation was made when the hydrogen was replaced by a fluorine atom (compounds 16 pIC50 = 4.86 vs. compound 6 pIC50 = 5.03)
The contour maps of the steric field (Figures 6A) or the electrostatic field (Figures 6B) of CoMSIA are nearly identical to those of the CoMFA model. Thus, the CoMSIA fractions (Table 1) indicate that the major fields potentially characterizing inhibitory activity are steric (64.2%) and electrostatic (35.8%), which are the same fields found in COMFA.
Newly designed compounds
Several novel triazolopyrazine derivatives were designed to increase activity considering the 3D QSAR results and the predominant modified region determined by the contour map results, using molecule 22 as a template (Figure 7).
FIGURE 7. The 3D-QSAR analysis revealed new insights into the structure-activity relationship.
In the case of the newly predicted molecules T01, T02, and T03, the replacement of methyl by isopropyl, phenyl, and cyclobutyl, respectively, in R1 increases the activity (Table 2). In addition, simultaneous substitution of R1 and radical Z increases the activity from 5.83 to 6.35, which is the case for T04 and T06; A4 (R1 = isobutyl + Z = fluorite); T05 (Z = pyrroline fluoride); and T06 (R1 = isobutyl + radical Z = fluorine). All these molecules blocked VEGFR-2 much better than compound 22, suggesting that they could be used as new VEGFR-2 inhibitors.
TABLE 2. The 3D-QSAR model’s predicted activity pIC50 (Pred).
ADME/toxicity prediction and analysis
Due to the lack of permeability of the blood-brain barrier, toxicity problems, and lack of efficacy, most drugs fail in the research and development phase. Therefore, prediction and optimization of ADMET properties of novel chemical entities are crucial to avoid potential problems during clinical trials. pKCSM (Pires et al., 2015) and the Swiss ADME (Daina et al., 2017) Predictor are commonly used programs to predict ADMET properties of small molecules. The in silico ADMET results for the six newly designed compounds are given in Table 3.
TABLE 3. ADME parameters of newly designed compounds.
The in silico ADMET results for the six newly designed compounds are given in Table 3. The ability of a drug to be absorbed into the bloodstream after being administered is directly related to its water solubility. Absorption values for all compounds are greater than 70% (100% for T03), indicating a greater absorption potential in the human intestine (Table 3). A drug travels to various organs and systems once it enters the bloodstream. Because of the blood-brain barrier, most substances are unable to cross over from the bloodstream into the central nervous system. This barrier is the primary one that stands between the bloodstream and the CNS. In order to enter the central nervous system (CNS), a drug molecule must therefore satisfy a number of requirements before it can cross the blood-brain barrier (BBB). These prerequisites include blood-brain permeability surface product (logPS) greater than −2 and a logarithmic ratio brain-to-plasma drug concentration (logBB) greater than 0.3. As shown in Table 1, designed compounds meet these values. In addition, the activity of an isoenzyme may be decreased or increased, depending on the drug. In some cases, the metabolism of a drug requires more than one isoenzyme. Approximately 90% of commonly prescribed drugs require the involvement of four isoenzymes for their metabolism. These isoenzymes are referred to as CYP1A2, CYP2C9, CYP2D6, and CYP3A4.
The in-silico Swiss-ADME prediction shows that all expected compounds are both substrates and inhibitors for 3A4, but not for 2D6. The results of this assay showed that none of the chemicals were mutagenic or carcinogenic. In addition, the acute toxicity (LD50) of the new compounds ranged from 2.423 to 2.923 mol/kg, and none of the candidates caused skin sensitization (Table 4). These proposed compounds are reliable candidates for further clinical studies because they exhibit exciting properties such as high intestinal absorption, distribution, permeability, and toxicity across the blood-brain barrier.
TABLE 4. Toxicity profile of newly designed compounds.
Docking results
The binding pattern of the novel designed inhibitors to VEGFR-2 was analyzed using molecular docking and compared to the mode of action of the standard inhibitor (Foretinib). To validate our approach and parameters before docking all chemicals, we redocked the native ligand sorafenib of VEGFR-2 into the binding pocket. The redocked conformation of sorafenib is very similar to its original conformation; the RMSD between the two configurations is only 1.2 to 1.01 (Supplementary Figure S1). Once the binding site was determined, docking was performed for each compound, including the active compound, proposed compounds, and Foretinib using the identical docking settings. The 3D and 2D binding interactions of all compounds are presented in Table 5.
TABLE 5. The 3D and 2D binding interactions showed that all proposed molecules interacted with a higher number of residues via hydrophobic interactions than Foretinib and compound 22.
Foretinib, a receptor tyrosine kinase inhibitor, has the potential for treating breast cancer by acting on multiple kinases involved in cancer growth. Compared to Tamoxifen and Trastuzumab, which target specific pathways, Foretinib offers a broader spectrum of activity. It formed two H-bonds during docking with residues Cys919, and Asp1046 (Table 5), and several hydrophobic and electrostatic interactions with Ala866, Leu840, Val848, Val899, Leu889, Leu1019, His1026, Asp1046 and Leu1035. In addition, Foretinib formed three carbon-hydrogen bonds with Ala881, His1026, and Asp1046. Molecule 22, the most active molecule in the data set, showed two conventional H-bond with Asp 1046, and Cys919, two carbon-H bonds due to the strong interaction of Cys1045 and Glu885, twelve hydrophobic via Pi-sigma, Pi-Pi T-shaped, P-alkyl and alkyl interactions with Cyc919, Leu1035, Leu840, Ala866, Phe918, Val848, Leu889, Leu1019, Cy1024, ILe892, Val899, Val916 and two Fluorine interactions with Cys1045 and Ile1044. By the same token, all complexes formed by designed compounds (T01-T06) and the VEGFR-2 receptor have better binding affinity ranging from −9.6 to −10 kcal/mol (Table 2). Therefore, these complexes are more stable than the complexes generated by Foretinib (−9 kcal/mol) and compound 22 (−8.8 kcal/mol) in the data set, suggesting that these compounds have a higher inhibitory potential against VEGFR-2. Furthermore, the three-dimensional binding interaction of the compounds (T1-T6) showed a similar H-binding interaction profile with Asp1046 and hydrophobic and electrostatic interactions with Val899, Val848, Ala866, Leu1035, Leu889, Leu1019, Leu1035, and Leu840, with compound 22, and Foretinib suggesting that these amino acids play a critical role in enhancing activity as reported in previous studies (Yousef et al., 2022; Parves et al., 2023). These compounds interact with a higher number of residues via hydrophobic interactions than molecule 22 and Foretinib, which increases their stability and affinity in the binding pocket of VEGFR-2.
Molecular dynamics simulation
Molecular dynamics modeling was used to assess the effects of molecule 22, the lead compound (Foretinib), and designed compound T01 with the highest predicted inhibitory activity (Table 2) on the structure of VEGFR-2 and its stability in the binding pocket (Agrahari et al., 2018; Agrahari et al., 2019). Accordingly, Gromacs 2019.3 was used to conduct a molecular dynamics simulation for these complexes (Figure 8), and various features such as root mean square deviation (RMSD), Root-mean square fluctuation (RMSF), and Radius of gyration (Rg) were examined on the trajectories corresponding data files.
FIGURE 8. The results of the molecular dynamics study: (A) Time evolution of the backbone of the target protein; (B) The comparative RMSF values for the target protein with the reference molecule, molecule 22, and designed molecule T01; (C) The comparative Radius of gyration values for the target protein with the reference molecule, molecule 22, and designed molecule T01; (E, F, G) The comparative hydrogen bonds and pairs within 0.35 nm for the target protein with the reference molecule, designed molecule 20, and designed molecule T01; (H) Hydrogen bond distributions of the reference molecule, molecule 22, and the designed molecule T01 with the target protein during the 100 ns. (D): The comparative SASA values for the target protein with reference molecule 22, and designed molecule T01.
Conformation of protein
The RMSD was computed during the simulation to establish the overall stability of the selected systems. The calculated value was considered as the main criterion for measuring the convergence of the system. The RMSD values of the 4ASD, 4ASD-Foretinib, mol22, and T01 complexes were calculated to be 0.17, 0.19, 0.18, and 0.18 nm, respectively. Up to 100 ns in the simulation, none of the three systems had a significant shift in root-mean-square-deviation (RMSD) values, which quantify conformational changes over time (Figures 8A). The resulting RMSD plot for 4ASD-Mol22 and 4ASD-T01 showed an increasing trend with increasing RMSD values between 0 and 25 ns, ranging from 0.10 to 0.30 nm, indicating that the compounds were adapting to a new conformation within the binding pocket (Aljuaid et al., 2022). Thereafter, the plateau continued until it reached a final value of 0.20 nm, which is below the threshold value of 0.3 nm. However, 4ASD-Foretinib has the same profile but not the same RMSD values as apo 4ASD. Finally, the lower RMSD values for all complexes studied indicate that the T01 inhibitor in 4ASD is stable, which provides a good basis for our investigation.
Root mean square fluctuation (RMSF)
Using the RMSF approach, we analyzed how ligand binding alters flexible protein structure and essential amino acid behavior. During simulation, a higher RMSF value indicates greater flexibility, while a lower RMSF value indicates greater rigidity (En-nahli et al., 2022). The RMSFs of the Apo 4ASD, 4ASD-Foretinib, 4ASD-mol22, and 4ASD-T01 complexes were calculated (Figures 8B). Compared with the apo form of 4ASD, the fluctuations of the residues in the ligand-bound complexes are quite stable, especially in the region where the residues participate in ligand binding. Furthermore, the average RMSF values of the Apo 4ASD, 4ASD-Foretinib, 4ASD-mol22, and 4ASD-T01 complexes were 0.10, 0.10, 0.09, and 0.08 nm, respectively. This result indicates that the binding of compound T01 contributes to the structural stability of 4ASD.
Radius of gyration analysis
We calculated the radius of gyration (Rg) as a function of time to investigate how the compactness of the protein structure changes when bound to various ligands. When Rg is sufficiently high, a ligand tends to be flexible, making it unstable. Conversely, conformations with lower Rg values tend to be dense and tightly packed (Naz et al., 2023). The average Rg values of the Apo 4ASD, 4ASD-Foretinib, 4ASD-mol22, and 4ASD-T01 complexes were 1.98, 2.00, 1.98, and 2.00 nm, respectively, suggesting that the binding of Foretinib, mol22, and the designed compound T01 to the 4ASD packing does not cause a significant change. As can be seen in Figure 8C, the Rg of the 4ASD/T01 complex appears to stabilize more rapidly during the 100-ns simulation than that of the 4ASD-Foretinib and 4ASD-mol22 complexes.
Solvent accessible surface area
SASA measures the surface area of the protein that is in direct contact with the solvent. The interpretation is based on the fact that the surface of the macromolecule-ligand complex is in contact with the water molecules surrounding it (Baammi et al., 2023b). The change of SASA for the complex of protein active molecule and for the complex of protein-designed molecules was analyzed during 100 ns (Figures 8D). The average values of SASA for the Apo 4ASD, 4ASD-Foretinib, 4ASD-mol22, and 4ASD-T01 complexes were 168.18 nm2, 169.00 nm2, 168.48 nm2, and 168.57 nm2, respectively. Analysis of these data sets revealed no significant variation in SASA values between complexes.
Hydrogen bonds analysis
Hydrogen bonding is an essential feature that determines binding affinity and contributes to the binding relationship between ligands and proteins. In drug discovery, it is also responsible for drug specificity, metabolization, and adsorption (Bhardwaj et al., 2020). To confirm the stability of all the docked complexes, the hydrogen bonds between 4ASD-Foretinib, 4ASD-mol22, and 4ASD-T01 were estimated in a solvent environment during MD simulations (Figures 8E, F, G). It was found that Foretinib formed an average of 2.01 hydrogen bonds and 3.23 bond pairs within 0.35 nm of the active pocket of 4ASD. Similarly, the designed molecule T01 was linked to 4ASD in the binding site via an average of 3.97 hydrogen bonds, while the average number of pairs within 0.35 nm was 4.15. However, for the 4ASD/molecule 22 complex, the average number of hydrogen bonds was 2.09, and the average number of pairs within 0.35 nm was 3.45. Notably, the H-bonding plot revealed that compound T01 was likely to interact more strongly with the binding pockets of 4ASD throughout the simulation compared with molecule 22 and Foretinib. The Hydrogen bond analysis emphasized the significant roles played by specific amino acid residues, in addition to catalytic residues, in the complexes of Foretinib, Molecule 22, and the developed compound T01 (Table 6). The distribution of hydrogen bond numbers further showed that the complex of designed compound T01 formed hydrogen bonds with affinities ranging from high to low, which is comparable to the distribution of hydrogen bonds in the complexes of Foretinib (Figures 8H).
TABLE 6. The hydrogen bond occupancy of amino acid residues throughout the simulation in various protein-ligand complexes.
Binding free energy analysis
The MM-PBSA method was used to determine the binding free energy (ΔE) between the VEGFR-2-Foretinib, VEGFR-2-molecule 22, and VEGFR-2-T01 complexes using the MmPbStat.py script for whole trajectories (Chen et al., 2016; Gomari et al., 2023). The total nonpolar, polar, and non-bonded interaction energies (electrostatic interaction and Van der Waals) were calculated for each complex and are displayed in Table 7. Foretinib, mol22, and T01 all bind to VEGFR-2 with free energies of −48.084, −34.943, and −59.176 kJ/mol, respectively, proving the validity of the molecular dynamic simulation model used in this study. The non-polar solvation free energy (Enon polar), electrostatic energy (Eele), and van der Waals energy (Evdw) all contributed to the binding energy of the two systems, but the polar energy (E polar) was undesirable, demonstrating the significance of the intermolecular van der Waal contribution. This is consistent with the docking study and MD simulation interactions, where the large interaction of the ligand with the hydrophobic binding pocket was observed.
TABLE 7. MMPBSA calculations of binding free energy for all complexes.
Conclusion
To build the 3D-QSAR model, a series of triazolopyrazine derivatives against the breast cancer cell line MCF-7 were collected, optimized, and calculated. Statistically, both the CoMFA and CoMSIA models provide good results, with R2 > 0.9 and Q2 > 0.5. External validation and the Y-randomization test were used to compare the predictive quality of the 3D-QSAR model. Thus, we used 3D-QSAR to design and predict the properties of 6 novel compounds. The results show that the expected activity and ADME/T curves for these molecules are quite strong. The molecular docking results show better binding affinity in the range of—8.9 to–10 kcal/mol, respectively, and strong binding to VEGFR-2 through several interactions. The MD simulation was used to study the stability of the conformations with the lowest binding value of each complex. 4ASD/T01 are stable based on RMSD, RMSF, Rg and SASA. During this research, the calculation of MM-GBSA confirmed the result of molecular docking by showing that the novel compound Pred T01 is more stable and has the lowest binding energy. Compared to the leading known breast cancer drug (Foretinib) the proposed six molecules exhibit enhanced binding and inhibitory activity against VEGFR-2, the major breast cancer receptor. This result provided the basis for the synthesis of novel triazolopyrazine analogs as improved drugs which expand the number of the new potential agents to face the rising resistant breast cancer. Chemical synthesis and further experimental validation to assess the proliferative activities of the designed molecule represents the main limitation of the present study.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.
Author contributions
SB: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Validation, Writing–original draft, Writing–review and editing, Software. AE: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Validation, Writing–original draft, Writing–review and editing, Resources, Supervision. RD: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Resources, Supervision, Validation, Writing–original draft, Writing–review and editing.
Funding
The authors declare that no financial support was received for the research, authorship, and/or publication of this article.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2023.1288652/full#supplementary-material
References
Abraham, M. J., Murtola, T., Schulz, R., Páll, S., Smith, J. C., Hess, B., et al. (2015). GROMACS: high performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 1–2, 19–25. doi:10.1016/J.SOFTX.2015.06.001
Acharya, C., Coop, A., Polli, J. E., and MacKerell, A. D. (2011). Recent advances in ligand-based drug design: relevance and utility of the conformationally sampled pharmacophore approach. Curr. Comput. Aided Drug Des. 7, 10–22. doi:10.2174/157340911793743547
Adeniji, S. E., Uba, S., and Uzairu, A. (2020). In silico study for evaluating the binding mode and interaction of 1, 2, 4-triazole and its derivatives as potent inhibitors against Lipoate protein B (LipB). J. King Saud. Univ. Sci. 32, 475–485. doi:10.1016/J.JKSUS.2018.07.014
Agrahari, A. K., Doss, G. P. C., Siva, R., Magesh, R., and Zayed, H. (2019). Molecular insights of the G2019S substitution in LRRK2 kinase domain associated with Parkinson’s disease: a molecular dynamics simulation approach. J. Theor. Biol. 469, 163–171. doi:10.1016/J.JTBI.2019.03.003
Agrahari, A. K., Kumar, A., Siva, R., Zayed, H., and George Priya Doss, C. (2018). Substitution impact of highly conserved arginine residue at position 75 in GJB1 gene in association with X-linked Charcot-Marie-tooth disease: a computational study. J. Theor. Biol. 437, 305–317. doi:10.1016/J.JTBI.2017.10.028
Aljuaid, A., Salam, A., Almehmadi, M., Baammi, S., Alshabrmi, F. M., Allahyani, M., et al. (2022). Structural homology-based drug repurposing approach for targeting NSP12 SARS-CoV-2. Molecules 27, 7732. doi:10.3390/molecules27227732
Aziz, M. A., Serya, R. A. T., Lasheen, D. S., Abdel-Aziz, A. K., Esmat, A., Mansour, A. M., et al. (2016). Discovery of potent VEGFR-2 inhibitors based on furopyrimidine and thienopyrimidne scaffolds as cancer targeting agents. Sci. Rep. 6, 24460. doi:10.1038/SREP24460
Baammi, S., Daoud, R., and el Allali, A. (2022). Assessing the effect of a series of mutations on the dynamic behavior of phosphite dehydrogenase using molecular docking, molecular dynamics and quantum mechanics/molecular mechanics simulations. J. Biomol. Struct. Dyn. 41, 4154–4166. doi:10.1080/07391102.2022.2064912
Baammi, S., Daoud, R., and El Allali, A. (2023b). In silico protein engineering shows that novel mutations affecting NAD+ binding sites may improve phosphite dehydrogenase stability and activity. Sci. Rep. 13, 1878–1913. doi:10.1038/s41598-023-28246-3
Baammi, S., El Allali, A., and Daoud, R. (2023a). Unleashing Nature’s potential: a computational approach to discovering novel VEGFR-2 inhibitors from African natural compound using virtual screening, ADMET analysis, molecular dynamics, and MMPBSA calculations. Front. Mol. Biosci. 10, 1227643. doi:10.3389/FMOLB.2023.1227643
Beheshti, A., Pourbasheer, E., Nekoei, M., and Vahdani, S. (2016). QSAR modeling of antimalarial activity of urea derivatives using genetic algorithm–multiple linear regressions. J. Saudi Chem. Soc. 20, 282–290. doi:10.1016/J.JSCS.2012.07.019
Benfenati, E., Chrétien, J. R., Gini, G., Piclin, N., Pintore, M., and Roncaglioni, A. (2007). “Validation of the models,” in Quantitative structure-activity relationships (QSAR) for pesticide regulatory purposes (Elsevier Inc), 185–199. doi:10.1016/B978-044452710-3/50008-2
Bhardwaj, V. K., Singh, R., Sharma, J., Rajendran, V., Purohit, R., and Kumar, S. (2020). Identification of bioactive molecules from tea plant as SARS-CoV-2 main protease inhibitors. J. Biomol. Struct. Dyn. 39, 3449–3458. doi:10.1080/07391102.2020.1766572
Bouamrane, S., Khaldan, A., Hajji, H., El-mernissi, R., Alaqarbeh, M., Alsakhen, N., et al. (2022). In silico identification of 1,2,4-triazoles as potential Candida Albicans inhibitors using 3D-QSAR, molecular docking, molecular dynamics simulations, and ADMET profiling. Mol. Divers 27, 2111–2132. doi:10.1007/S11030-022-10546-X
Burley, S. K., Berman, H. M., Bhikadiya, C., Bi, C., Chen, L., di Costanzo, L., et al. (2019). RCSB Protein Data Bank: biological macromolecular structures enabling research and education in fundamental biology, biomedicine, biotechnology and energy. Nucleic Acids Res. 47, D464–D474. doi:10.1093/NAR/GKY1004
Chavakis, E., and Dimmeler, S. (2002). Regulation of endothelial cell survival and apoptosis during angiogenesis. Arterioscler. Thromb. Vasc. Biol. 22, 887–893. doi:10.1161/01.ATV.0000017728.55907.A9
Chen, D., Oezguen, N., Urvil, P., Ferguson, C., Dann, S. M., and Savidge, T. C. (2016). Regulation of protein-ligand binding affinity by hydrogen bond pairing. Sci. Adv. 2, e1501240. doi:10.1126/SCIADV.1501240
Cherkasov, A., Muratov, E. N., Fourches, D., Varnek, A., Baskin, I. I., Cronin, M., et al. (2014). QSAR Modeling: where have you been? Where are you going to? J. Med. Chem. 57, 4977–5010. doi:10.1021/JM4004285
Chirico, N., and Gramatica, P. (2011). Real external predictivity of QSAR models: how to evaluate it? Comparison of different validation criteria and proposal of using the concordance correlation coefficient. J. Chem. Inf. Model. 51, 2320–2335. doi:10.1021/CI200211N
Chu, H., He, Q. X., Wang, J., Hu, Y., Wang, Y. Q., and Lin, Z. H. (2020). In silico design of novel benzohydroxamate-based compounds as inhibitors of histone deacetylase 6 based on 3D-QSAR, molecular docking, and molecular dynamics simulations. New J. Chem. 44, 21201–21210. doi:10.1039/D0NJ04704J
Consonni, V., Ballabio, D., and Todeschini, R. (2009). Comments on the definition of the Q2 parameter for QSAR validation. J. Chem. Inf. Model. 49, 1669–1678. doi:10.1021/ci900115y
Daina, A., Michielin, O., and Zoete, V. (2017). SwissADME: a free web tool to evaluate pharmacokinetics, drug-likeness and medicinal chemistry friendliness of small molecules. Sci. Rep. 7, 42717. doi:10.1038/SREP42717
Edraki, N., Das, U., Hemateenejad, B., Dimmock, J. R., and Miri, R. (2016). Comparative QSAR analysis of 3,5-bis (Arylidene)-4-Piperidone derivatives: the development of predictive cytotoxicity models. Iran. J. Pharm. Res. 15, 425–437.
el Mchichi, L., Tabti, K., Kasmi, R., El-Mernissi, R., el Aissouq, A., En-nahli, F., et al. (2022). 3D-QSAR study, docking molecular and simulation dynamic on series of benzimidazole derivatives as anti-cancer agents. J. Indian Chem. Soc. 99, 100582. doi:10.1016/J.JICS.2022.100582
En-nahli, F., Baammi, S., Hajji, H., Alaqarbeh, M., Lakhlifi, T., and Bouachrine, M. (2022). High-throughput virtual screening approach of natural compounds as target inhibitors of plasmepsin-II. J. Biomol. Struct. Dyn. 41, 10070–10080. doi:10.1080/07391102.2022.2152871
Gomari, M. M., Tarighi, P., Choupani, E., Abkhiz, S., Mohamadzadeh, M., Rostami, N., et al. (2023). Structural evolution of Delta lineage of SARS-CoV-2. Int. J. Biol. Macromol. 226, 1116–1140. doi:10.1016/J.IJBIOMAC.2022.11.227
Hadni, H., and Elhallaoui, M. (2020). 3D-QSAR, docking and ADMET properties of aurone analogues as antimalarial agents. Heliyon 6, e03580. doi:10.1016/J.HELIYON.2020.E03580
Hajjo, R., Grulke, C. M., Golbraikh, A., Setola, V., Huang, X. P., Roth, B. L., et al. (2010). Development, validation, and use of quantitative structure-activity relationship models of 5-hydroxytryptamine (2B) receptor ligands to identify novel receptor binders and putative valvulopathic compounds among common drugs. J. Med. Chem. 53, 7573–7586. doi:10.1021/JM100600Y
Hernández-Santoyo, A., Tenorio-Barajas, A. Y., Altuzar, V., Vivanco-Cid, V., and Mendoza-Barrera, C. (2013) “Ogawa, Protein-protein and protein-ligand docking” in Protein engineering–technology and applications. (United Kingdom) vol. 196. doi:10.5772/56376
Hicklin, D. J., and Ellis, L. M. (2005). Role of the vascular endothelial growth factor pathway in tumor growth and angiogenesis. J. Clin. Oncol. 23, 1011–1027. doi:10.1200/JCO.2005.06.081
Kesari, P., Pratap, S., Dhankhar, P., Dalal, V., Mishra, M., Singh, P. K., et al. (2020). Structural characterization and in-silico analysis of Momordica charantia 7S globulin for stability and ACE inhibition. Sci. Rep. 10, 1160. doi:10.1038/S41598-020-58138-9
Lindahl, E., Bjelkmar, P., Larsson, P., Cuendet, M. A., and Hess, B. (2010). Implementation of the CHARMM force field in GROMACS: analysis of protein stability effects from correction maps, virtual interaction sites, and water models. J. Chem. Theory Comput. 6, 459–466. doi:10.1021/CT900549R
Liu, X., Li, Y., Zhang, Q., Pan, Q., Zheng, P., Dai, X., et al. (2022). Design, synthesis, and biological evaluation of [1,2,4]triazolo[4,3-a] pyrazine derivatives as novel dual c-met/VEGFR-2 inhibitors. Front. Chem. 10, 815534. doi:10.3389/fchem.2022.815534
Lu, S., Liu, H. C., Chen, Y. D., Yuan, H. L., Sun, S. L., Gao, Y. P., et al. (2011). Combined pharmacophore modeling, docking, and 3D-QSAR studies of PLK1 inhibitors. Int. J. Mol. Sci. 12, 8713–8739. doi:10.3390/IJMS12128713
Meng, X. Y., Zhang, H. X., Mezei, M., and Cui, M. (2011). Molecular Docking: a powerful approach for structure-based drug discovery. Curr. Comput. Aided Drug Des. 7, 146–157. doi:10.2174/157340911795677602
Modi, S. J., Modh, D. H., and Kulkarni, V. M. (2020). Insights into the structural features of anticancer 1,6-naphthyridines and pyridopyrimidines as VEGFR-2 inhibitors: 3D-QSAR study. J. Appl. Pharm. Sci. 10, 001–022. doi:10.7324/JAPS.2020.10101
Moussaoui, M., Baassi, M., Baammi, S., Soufi, H., Salah, M., Daoud, R., et al. (2023). In silico design of novel CDK2 inhibitors through QSAR, ADMET, molecular docking and molecular dynamics simulation studies. J. Biomol. Struct. Dyn. 2023, 1–17. doi:10.1080/07391102.2023.2212304
Muhammad, A., Forcados, G. E., Katsayal, B. S., Bako, R. S., Aminu, S., Sadiq, I. Z., et al. (2022). Potential epigenetic modifications implicated in triple-to quadruple-negative breast cancer transition: a review. Epigenomics 14, 711–726. doi:10.2217/EPI-2022-0033
Naz, F., Khan, I., Baammi, S., and Islam, A. (2023). Investigation of the interactions of HSA and SARS-CoV-2 papain-like protease against eugenol for novel COVID-19 drug discovery: spectroscopic and insilico study. J. Biomol. Struct. Dyn. 41, 10161–10170. doi:10.1080/07391102.2022.2164062
Netzeva, T. I., Worth, A. P., Aldenberg, T., Benigni, R., Cronin, M. T. D., Gramatica, P., et al. (2005). Current status of methods for defining the applicability domain of (quantitative) structure-activity relationships. The report and recommendations of ECVAM Workshop 52. Altern. Lab. Anim. 33, 155–173. doi:10.1177/026119290503300209
Nicolaou, K. C., Wu, T. R., Sarlah, D., Shaw, D. M., Rowcliffe, E., and Burton, D. R. (2008). Total synthesis, revised structure, and biological evaluation of biyouyanagin A and analogues thereof. J. Am. Chem. Soc. 130, 11114–11121. doi:10.1021/ja802805c
Parves, M. R., Riza, Y. M., Alam, S., and Jaman, S. (2023). Molecular dynamics-based insight of VEGFR-2 kinase domain: a combined study of pharmacophore modeling and molecular docking and dynamics. J. Mol. Model. 29, 17–15. doi:10.1007/s00894-022-05427-x
Pires, D. E. V., Blundell, T. L., and Ascher, D. B. (2015). pkCSM: predicting small-molecule pharmacokinetic and toxicity properties using graph-based signatures. J. Med. Chem. 58, 4066–4072. doi:10.1021/acs.jmedchem.5b00104
Qin, L., Zhang, X., Chen, Y., Mo, L., Zeng, H., and Liang, Y. (2017). Predictive QSAR models for the toxicity of disinfection byproducts. Molecules 22, 1671. doi:10.3390/MOLECULES22101671
Roy, K. K., Kar, S., and Das, R. N. (2015). “Validation of QSAR models,” in Understanding the basics of QSAR for applications in pharmaceutical sciences and risk assessment (Elsevier Inc), 231–289. doi:10.1016/B978-0-12-801505-6.00007-7
Roy, K., and Mitra, I. (2011). On various metrics used for validation of predictive QSAR models with applications in virtual screening and focused library design. Comb. Chem. High. Throughput Screen 14, 450–474. doi:10.2174/138620711795767893
Rücker, C., Rücker, G., and Meringer, M. (2007). Y-randomization and its variants in QSPR/QSAR. J. Chem. Inf. Model. 47, 2345–2357. doi:10.1021/ci700157b
Sahigara, F., Mansouri, K., Ballabio, D., Mauri, A., Consonni, V., and Todeschini, R. (2012). Comparison of different approaches to define the applicability domain of QSAR models. Molecules 17, 4791–4810. doi:10.3390/MOLECULES17054791
Saini, G., Dalal, V., Gupta, D. N., Sharma, N., Kumar, P., and Sharma, A. K. (2021). A molecular docking and dynamic approach to screen inhibitors against ZnuA1 of Candidatus Liberibacter asiaticus. Mol. Simul. 47, 510–525. doi:10.1080/08927022.2021.1888948
Sharma, R., Dhingra, N., and Patil, S. (2016). CoMFA, CoMSIA, HQSAR and molecular docking analysis of ionone-based chalcone derivatives as antiprostate cancer activity. Indian J. Pharm. Sci. 78, 54–64. doi:10.4103/0250-474X.180251
Shibuya, M. (2010). Tyrosine kinase receptor flt/VEGFR family: its characterization related to angiogenesis and cancer. Genes. Cancer 1, 1119–1123. doi:10.1177/1947601910392987
Shibuya, M. (2011). Vascular endothelial growth factor (VEGF) and its receptor (vegfr) signaling in angiogenesis: a crucial target for anti- and pro-angiogenic therapies. Genes. Cancer 2, 1097–1105. doi:10.1177/1947601911423031
Shivanika, C., Deepak Kumar, S., Ragunathan, V., Tiwari, P., Sumitha, A., and Brindha Devi, P. (2020). Molecular docking, validation, dynamics simulations, and pharmacokinetic prediction of natural compounds against the SARS-CoV-2 main-protease. J. Biomol. Struct. Dyn. 40, 1–27. doi:10.1080/07391102.2020.1815584
Singh, V., Dhankhar, P., Dalal, V., Tomar, S., and Kumar, P. (2022). In-silico functional and structural annotation of hypothetical protein from Klebsiella pneumonia: a potential drug target. J. Mol. Graph Model. 116, 108262. doi:10.1016/J.JMGM.2022.108262
Soufi, H., Moussaoui, M., Baammi, S., Baassi, M., Salah, M., Daoud, R., et al. (2023). Multi-combined QSAR, molecular docking, molecular dynamics simulation, and ADMET of Flavonoid derivatives as potent cholinesterase inhibitors. J. Biomol. Struct. Dyn. 2023, 1–15. doi:10.1080/07391102.2023.2238314
Tabti, K., Baammi, S., ElMchichi, L., Sbai, A., Maghat, H., Bouachrine, M., et al. (2022). Computational investigation of pyrrolidin derivatives as novel GPX4/MDM2–p53 inhibitors using 2D/3D-QSAR, ADME/toxicity, molecular docking, molecular dynamics simulations, and MM-GBSA free energy. Struct. Chem. 33, 1019–1039. doi:10.1007/s11224-022-01903-5
Tosco, P., and Mackey, M. (2017). “Lessons and successes in the use of molecular fields, comprehensive medicinal chemistry III,” in Reference module in chemistry, molecular sciences and chemical engineering. Rotella, and D. Ward (Elsevier (GmbH) 3–8, 253–296. doi:10.1016/B978-0-12-409547-2.12353-4
Tropsha, A., Gramatica, P., and Gombar, V. K. (2003). The importance of being earnest: validation is the absolute essential for successful application and interpretation of QSPR models. QSAR Comb. Sci. 22, 69–77. doi:10.1002/QSAR.200390007
Tsai, K. C., Chen, Y. C., Hsiao, N. W., Wang, C. L., Lin, C. L., Lee, Y. C., et al. (2010). A comparison of different electrostatic potentials on prediction accuracy in CoMFA and CoMSIA studies. Eur. J. Med. Chem. 45, 1544–1551. doi:10.1016/J.EJMECH.2009.12.063
Vistoli, G., and Pedretti, A. (2007). Molecular fields to assess recognition forces and property spaces. Compr. Med. Chem. II 5, 577–602. doi:10.1016/B0-08-045044-X/00142-5
Vugmeyster, Y., Harrold, J., and Xu, X. (2012). Absorption, distribution, metabolism, and excretion (ADME) studies of biotherapeutics for autoimmune and infla mmatory conditions. AAPS J. 14, 714–727. doi:10.1208/S12248-012-9385-Y
Waller, C. L., Oprea, T. I., Chae, K., Park, H. K., Korach, K. S., Laws, S. C., et al. (1996). Ligand-based identification of environmental estrogens. Chem. Res. Toxicol. 9, 1240–1248. doi:10.1021/TX960054F
Wang, B., Shen, J., Wang, Z., Liu, J., Ning, Z., and Hu, M. (2018a). Isomangiferin, a novel potent vascular endothelial growth factor receptor 2 kinase inhibitor, suppresses breast cancer growth, metastasis and angiogenesis. J. Breast Cancer 21, 11–20. doi:10.4048/JBC.2018.21.1.11
Wang, M. Y., Liang, J. W., Olounfeh, K. M., Sun, Q., Zhao, N., and Meng, F. H. (2018b). A comprehensive in silico method to study the QSTR of the aconitine alkaloids for designing novel drugs. A J. Synthetic Chem. Nat. Prod. Chem. 23, 2385. doi:10.3390/MOLECULES23092385
Wang, X., Bove, A. M., Simone, G., and Ma, B. (2020). Molecular bases of VEGFR-2-mediated physiological function and pathological role. Front. Cell. Dev. Biol. 8, 599281. doi:10.3389/fcell.2020.599281
Wedge, S. R., Kendrew, J., Hennequin, L. F., Valentine, P. J., Barry, S. T., Brave, S. R., et al. (2005). AZD2171: a highly potent, orally bioavailable, vascular endothelial growth factor receptor-2 tyrosine kinase inhibitor for the treatment of cancer. Cancer Res. 65, 4389–4400. doi:10.1158/0008-5472.CAN-04-4409
Wen, C., Yan, X., Gu, Q., Du, J., Wu, D., Lu, Y., et al. (2019). Systematic studies on the protocol and criteria for selecting a covalent docking tool. Molecules 24, 2183. doi:10.3390/MOLECULES24112183
WHO (2020). Breast cancer. Avaliable at: https://www.who.int/news-room/fact-sheets/detail/breast-cancer (Accessed September 29, 2023).
Wu, F., Zhou, Y., Li, L., Shen, X., Chen, G., Wang, X., et al. (2020). Computational approaches in preclinical studies on drug discovery and development. Front. Chem. 8, 726. doi:10.3389/fchem.2020.00726
Yousef, R. G., Elwan, A., Gobaara, I. M. M., Mehany, A. B. M., Eldehna, W. M., El-Metwally, S. A., et al. (2022). Anti-cancer and immunomodulatory evaluation of new nicotinamide derivatives as potential VEGFR-2 inhibitors and apoptosis inducers: in vitro and in silico studies. J. Enzyme Inhib. Med. Chem. 37, 2206–2222. doi:10.1080/14756366.2022.2110868
Zhao, X., Chen, M., Huang, B., Ji, H., and Yuan, M. (2011). Comparative molecular field analysis (CoMFA) and comparative molecular similarity indices analysis (CoMSIA) studies on α(1A)-adrenergic receptor antagonists based on pharmacophore molecular alignment. Int. J. Mol. Sci. 12, 7022–7037. doi:10.3390/IJMS12107022
Ziyad, S., and Iruela-Arispe, M. L. (2011). Molecular mechanisms of tumor angiogenesis. Genes. Cancer 2, 1085–1096. doi:10.1177/1947601911432334
Keywords: breast cancer, VEGFR-2, 3D-QSAR, ADMET, triazolopyrazine, molecular docking, molecular dynamic simulations, MMPBSA calculation
Citation: Baammi S, El Allali A and Daoud R (2023) Potent VEGFR-2 inhibitors for resistant breast cancer: a comprehensive 3D-QSAR, ADMET, molecular docking and MMPBSA calculation on triazolopyrazine derivatives. Front. Mol. Biosci. 10:1288652. doi: 10.3389/fmolb.2023.1288652
Received: 04 September 2023; Accepted: 06 November 2023;
Published: 22 November 2023.
Edited by:
Vikram Dalal, Washington University in St. Louis, United StatesReviewed by:
Poonam Dhankhar, Cornell University, United StatesGourav Choudhir, Chaudhary Charan Singh University, India
Gunjan Saini, Purdue University, United States
Ravi Rathi, Christ University, India
Copyright © 2023 Baammi, El Allali and Daoud. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Achraf El Allali, achraf.elallali@um6p.ma; Rachid Daoud, rachid.daoud@um6p.ma