Integration of Hybridization Strategies in Pyridine–Urea Scaffolds for Novel Anticancer Agents: Design, Synthesis, and Mechanistic Insights

Annually, millions of new cancer cases are reported, leading to millions of deaths worldwide. Among the newly reported cases, breast and colon cancers prevail as the most frequently detected variations. To effectively counteract this rapid increase, the development of innovative therapies is crucial. Small molecules possessing pyridine and urea moieties have been reported in many of the currently available anticancer agents, especially VEGFR2 inhibitors. With this in mind, a rational design approach was employed to create hybrid small molecules combining urea and pyridine. These synthesized compounds underwent in vitro testing against breast and colon cancer cell lines, revealing potent submicromolar anticancer activity. Compound 8a, specifically, exhibited an impressive GI50 value of 0.06 μM against the MCF7 cancer cell line, while compound 8h displayed the highest cytotoxic activity against the HCT116 cell line, with a GI50 of 0.33 ± 0.042 μM. Notably, compounds 8a, 8h, and 8i demonstrated excellent safety profiles when tested on normal cells. Molecular docking, dynamic studies, and free energy calculations were employed to validate the affinity of these compounds as VEGFR2 inhibitors.


Introduction
Cancer is a serious public health concern around the world, and it is the second leading cause of mortality after cardiovascular diseases [1]. Cancer patients suffer from high mortality rates due to late diagnosis and delayed medical treatment [2]. According to recent estimates, there were 10.9 million newly reported cancer diagnoses, resulting in 6.7 million deaths, and 24.6 million people were still living with cancer 3 years after their diagnosis. Breast cancer and colorectal cancer were two of the most commonly diagnosed cancers [3,4]. While the overexpression of various proteins has been linked to the development and growth of cancer, VEGFR2 has emerged as a key therapeutic target due to its role in the formation of new blood vessels as well as the regulation of endothelial differentiation of colon cancer cells. Moreover, VEGFR2 inhibition has been found to hamper breast cancer cell proliferation via enhanced mitochondrial biogenesis.
Despite the availability of approved therapy for VEGFR2, there is a significant rise in multidrug resistance to the current drugs, necessitating the development of novel therapeutic options [5,6]. In order to accelerate drug discovery, several techniques have emerged, among which hybridization is one of the most frequently used. Molecular hybridization is a promising strategy that involves the combination of two or more pharmacophoric moieties in a single molecular architecture with improved activity and affinity with reference to their parent molecules [7]. It has been employed in the design of anticancer drugs by linking one anticancer drug with another anticancer drug or a carrier molecule strategically. This results in increased selectivity, fewer side effects, reduction in drug resistance, better patient compliance, and more biological potency [8,9].
Most of the drugs that have been developed to combat cancer and its spread possess heterocycles in their backbone with non-fused pyridines emerging as one of the main building blocks for the development of anticancer drugs (Figure 1) [10][11][12]. Consequently, many of the drugs developed to combat cancer share the common feature of heterocyclics in their backbone, with non-fused pyridines playing a prominent role [13,14]. This class of heterocycles is known to possess diverse biological activities, particularly anticancer activities.
Molecules 2023, 28, x FOR PEER REVIEW 2 o differentiation of colon cancer cells. Moreover, VEGFR2 inhibition has been found to ha per breast cancer cell proliferation via enhanced mitochondrial biogenesis. Despite the availability of approved therapy for VEGFR2, there is a significant rise multidrug resistance to the current drugs, necessitating the development of novel the peutic options [5,6]. In order to accelerate drug discovery, several techniques ha emerged, among which hybridization is one of the most frequently used. Molecular bridization is a promising strategy that involves the combination of two or more pharm cophoric moieties in a single molecular architecture with improved activity and affin with reference to their parent molecules [7]. It has been employed in the design of an cancer drugs by linking one anticancer drug with another anticancer drug or a carr molecule strategically. This results in increased selectivity, fewer side effects, reduction drug resistance, be er patient compliance, and more biological potency [8,9].
Most of the drugs that have been developed to combat cancer and its spread poss heterocycles in their backbone with non-fused pyridines emerging as one of the m building blocks for the development of anticancer drugs (Figure 1) [10][11][12]. Consequen many of the drugs developed to combat cancer share the common feature of heterocyc in their backbone, with non-fused pyridines playing a prominent role [13,14]. This cl of heterocycles is known to possess diverse biological activities, particularly antican activities. Several pyridine-based small molecule VEGFR2 inhibitors have been approved anticancer drugs, such as sorafenib, axitinib, and regorafenib ( Figure 1) [13,15,16]. Th compounds also share the presence of a urea or amide linker (urido), with the urea moi being an essential pharmacophoric feature in several anticancer drugs such as sorafe and regorafenib [17]. With consideration of these pharmacophoric features, a hybridi tion approach was adopted in this study, utilizing pyridine-urea as the foundational c for the development of novel hybrid compounds.
Quinazolines are versatile heterocyclic compounds that have shown signific promise in anticancer research due to their unique chemical structure and potent cytoto activity against various cancer types [18,19]. Several compounds derived from quina line have gained approval for the treatment of various conditions, including benign pr tatic hyperplasia, pos raumatic stress disorder, as well as lung and pancreatic canc {Jafari, 2016 #115}. Based on the successful incorporation of the quinazoline ring in pre ous studies, as well as the compelling evidence supporting its efficacy, the pres Several pyridine-based small molecule VEGFR2 inhibitors have been approved as anticancer drugs, such as sorafenib, axitinib, and regorafenib ( Figure 1) [13,15,16]. These compounds also share the presence of a urea or amide linker (urido), with the urea moiety being an essential pharmacophoric feature in several anticancer drugs such as sorafenib and regorafenib [17]. With consideration of these pharmacophoric features, a hybridization approach was adopted in this study, utilizing pyridine-urea as the foundational core for the development of novel hybrid compounds.
Quinazolines are versatile heterocyclic compounds that have shown significant promise in anticancer research due to their unique chemical structure and potent cytotoxic activity against various cancer types [18,19]. Several compounds derived from quinazoline have gained approval for the treatment of various conditions, including benign prostatic hyperplasia, posttraumatic stress disorder, as well as lung and pancreatic cancers {Jafari, 2016 #115}. Based on the successful incorporation of the quinazoline ring in previous studies, as well as the compelling evidence supporting its efficacy, the present investigation chose to incorporate the quinazoline ring as the "head" moiety in the new scaffold attached to the pyridine ring [20].
While various aliphatic, cyclic, and aromatic substitutions were used as the "tail" moiety ( Figure 2). Unfortunately, the quinazoline-based hybrids did not exhibit satisfactory anticancer activity against breast and colon cancer cell lines, likely due to the bulky nature of the quinazoline ring. Predicting that the bulky nature of the quinazoline ring was the cause of the unsatisfactory anticancer activity, the quinazoline moiety was replaced with a less bulky one. Inspired by a recent study that reported the successful replacement of a benzene ring with an indazole to produce new anticancer derivatives and given that the indazole and quinazoline moieties contain the same number of nitrogen atoms, the quinazoline ring was replaced with an indazole one in the new hybrids ( Figure 2) [21]. This substitution led to a significant increase in cytotoxic activity against both breast and colon cancer cell lines.
investigation chose to incorporate the quinazoline ring as the "head" moiety in the new scaffold a ached to the pyridine ring [20].
While various aliphatic, cyclic, and aromatic substitutions were used as the "tail" moiety ( Figure 2). Unfortunately, the quinazoline-based hybrids did not exhibit satisfactory anticancer activity against breast and colon cancer cell lines, likely due to the bulky nature of the quinazoline ring. Predicting that the bulky nature of the quinazoline ring was the cause of the unsatisfactory anticancer activity, the quinazoline moiety was replaced with a less bulky one. Inspired by a recent study that reported the successful replacement of a benzene ring with an indazole to produce new anticancer derivatives and given that the indazole and quinazoline moieties contain the same number of nitrogen atoms, the quinazoline ring was replaced with an indazole one in the new hybrids ( Figure  2) [21]. This substitution led to a significant increase in cytotoxic activity against both breast and colon cancer cell lines.  The synthesized compounds were evaluated for their cytotoxic activity against HCT-116 and MCF-7 cancer cell lines, a method that has been employed to assess the cytotoxic activity of various VEGFR2 inhibitors [22,23]. Furthermore, the most active compounds were subjected to molecular docking and molecular dynamics (MD) studies as well as free energy calculations to ensure their binding affinity toward VEGFR2.

Chemistry
The synthesis of compounds 4a-n was carried out as illustrated in Scheme 1. The palladium-catalyzed coupling reaction of bis(pinacolato)diboron with commercially available 4-chloropyridin-2-amine (1) in the presence of potassium acetate and SPhos afforded boronic acid 2. The synthesis of boronic acid was afforded by following the reported procedure [24,25]. The installation of quinazoline was achieved using the Suzuki cross-coupling reaction of boronic acid 2 with 6-bromoquinazoline in the presence of Pd(dppf)Cl 2 .DCM to afford the amine 3 [26]. Finally, the target compounds 4a-n were synthesized by treating the precursor amine 3 with the corresponding isocyanate in the presence of N,N-diisopropylethylamine (DIPEA).
The synthesized compounds were evaluated for their cytotoxic activity against HCT-116 and MCF-7 cancer cell lines, a method that has been employed to assess the cytotoxic activity of various VEGFR2 inhibitors [22,23]. Furthermore, the most active compounds were subjected to molecular docking and molecular dynamics (MD) studies as well as free energy calculations to ensure their binding affinity toward VEGFR2.
3-Methoxy-1H-indazole derivatives 8a-k were synthesized using a similar approach with a slightly modified reaction order and procedure in the first and last steps, as shown.
In scheme 2, the synthesis started with the protection of 5-bromo-3-methoxy-1H-indazole (5) with Boc anhydride in the presence of triethylamine (Et3N) to afford 6. Unfortunately, the Suzuki coupling methodology employed in scheme 1 failed to afford compound 7. Accordingly, Xphos-pd-G2 and K2CO3 was used. The Boc-protected indazole 6 was converted to 7 using a Suzuki cross-coupling reaction with (2-aminopyridin-4-yl)boronic acid (2). The amine in 7 was treated using the appropriate isocyanate in the presence Scheme 1.
3-Methoxy-1H-indazole derivatives 8a-k were synthesized using a similar approach with a slightly modified reaction order and procedure in the first and last steps, as shown.
In Scheme 2, the synthesis started with the protection of 5-bromo-3-methoxy-1Hindazole (5) with Boc anhydride in the presence of triethylamine (Et 3 N) to afford 6. Unfortunately, the Suzuki coupling methodology employed in Scheme 1 failed to afford compound 7. Accordingly, Xphos-pd-G2 and K 2 CO 3 was used. The Boc-protected indazole 6 was converted to 7 using a Suzuki cross-coupling reaction with (2-aminopyridin-4yl)boronic acid (2). The amine in 7 was treated using the appropriate isocyanate in the presence of N,N-diisopropylethylamine (DIPEA) followed by removal of the Boc-protecting group with 4 N HCl, affording compounds 8a-k.

Biological Evaluation
To evaluate the anticancer activity of the synthesized compounds, the inhibitory effect of the synthesized compounds on cell proliferation in the form of GI50 was measured and compared to the positive control drug, irinotecan. Irinotecan (Supporting Information, Figure S68) was chosen as the standard drug, as it is a widely used FDA-approved chemotherapy drug for treating various types of cancer, including colorectal cancer, small cell lung cancer, pancreatic cancer, and gastric cancer, among others. For this study, the HCT116 and MCF7 cancer cell lines were selected due to their well-established nature and extensive use in evaluating various breast and colon cancer candidates. Furthermore, it was discovered that VEGFR-2 plays a critical role in cell survival and regulates endothelial differentiation in both MCF-7 breast cancer cells and HCT-116 human colorectal carcinoma cells [27][28][29][30].
The results showed that irinotecan exhibited GI50 values of 0.35 ± 0.069 and 0.62 ± 0.003 µM against the MCF7 and HCT116 cancer cell lines, respectively. Generally, the quinazoline-based compounds (4a-n) exhibited lower cytotoxic activity compared to the 3-methoxy-1H-indazole compounds (8a-k). Among the synthesized quinazoline-based compounds, only compounds 4i, 4k, and 4l showed a strong inhibitory effect (Table 1) on the breast cancer cell line (MCF7) with GI50 values of 3.35 ± 0.101, 3.03 ± 0.061, and 9.00 ± 0.347 µM, respectively. However, these three compounds did not exhibit any potent effect on the colon cancer cell line, HCT116, indicating that the quinazoline-based compounds are more specific to breast cancer.
On the other hand, with the exception of 3-fluoro-substituted derivative (compound 8b), the unsubstituted 3-methoxy-1H-indazole derivatives, namely compounds 8a, 8h, and 8i, displayed the most potent activity against both MCF7 and HCT116 cancer cell lines. Compound 8a exhibited GI50 values of 0.06 ± 0.014 µM and 1.19 ± 0.035 µM against MCF7 and HCT116 cell lines, respectively. This result means that compound 8a displayed 6-fold potency when compared to the positive control drug Irinotecan against the MCF7

Biological Evaluation
To evaluate the anticancer activity of the synthesized compounds, the inhibitory effect of the synthesized compounds on cell proliferation in the form of GI 50 was measured and compared to the positive control drug, irinotecan. Irinotecan (Supporting Information, Figure S68) was chosen as the standard drug, as it is a widely used FDA-approved chemotherapy drug for treating various types of cancer, including colorectal cancer, small cell lung cancer, pancreatic cancer, and gastric cancer, among others. For this study, the HCT116 and MCF7 cancer cell lines were selected due to their well-established nature and extensive use in evaluating various breast and colon cancer candidates. Furthermore, it was discovered that VEGFR-2 plays a critical role in cell survival and regulates endothelial differentiation in both MCF-7 breast cancer cells and HCT-116 human colorectal carcinoma cells [27][28][29][30].
The results showed that irinotecan exhibited GI 50 values of 0.35 ± 0.069 and 0.62 ± 0.003 µM against the MCF7 and HCT116 cancer cell lines, respectively. Generally, the quinazoline-based compounds (4a-n) exhibited lower cytotoxic activity compared to the 3-methoxy-1H-indazole compounds (8a-k). Among the synthesized quinazolinebased compounds, only compounds 4i, 4k, and 4l showed a strong inhibitory effect (Table 1) on the breast cancer cell line (MCF7) with GI 50 values of 3.35 ± 0.101, 3.03 ± 0.061, and 9.00 ± 0.347 µM, respectively. However, these three compounds did not exhibit any potent effect on the colon cancer cell line, HCT116, indicating that the quinazoline-based compounds are more specific to breast cancer.
On the other hand, with the exception of 3-fluoro-substituted derivative (compound 8b), the unsubstituted 3-methoxy-1H-indazole derivatives, namely compounds 8a, 8h, and 8i, displayed the most potent activity against both MCF7 and HCT116 cancer cell lines. Compound 8a exhibited GI 50 values of 0.06 ± 0.014 µM and 1.19 ± 0.035 µM against MCF7 and HCT116 cell lines, respectively. This result means that compound 8a displayed 6-fold potency when compared to the positive control drug Irinotecan against the MCF7 cancer cell line. Compounds 8h and 8i exhibited GI 50 values of 0.23 ± 0.015 µM and 0.16 ± 0.013 µM against the MCF7 cancer cell lines (Table 1). Furthermore, compound 8h showed a GI 50 value of 0.33 ± 0.042 µM against HCT116 cell line, which is twice the potency displayed by the FDA-approved drug irinotecan against the HCT116 cancer cell line (GI 50 = 0.62 ± 0.003). Overall, the 3-methoxy-1H-indazole compounds demonstrated potent cytotoxic activity against both MCF7 and HCT116 cell lines, particularly the unsubstituted derivatives 8a, 8h, and 8i. These findings suggest that 3-methoxy-1H-indazole compounds could be promising candidates for the development of novel anticancer drugs targeting breast and colon cancers. Table 1. GI 50 of the synthesized compounds over MCF7 and HCT116 cancer cell lines.

Compd.
GI 50  In order to assess the safety profile of the synthesized compounds, the compounds with the best cytotoxic activity (compounds were 8a, 8h, and 8i) were tested against the normal cell line, HEK293. In HEK293 cells, it was confirmed that the compounds exhibited lower cell growth inhibitory effects than the anticancer effects shown in MCF and HCT116, where the cell growth rate did not decrease by less than 50% compared to the control group as illustrated in Figure 3.

Molecular Docking
Molecular docking plays a significant role in computer-aided drug design and has diverse applications. One such application involves predicting the binding modes between a ligand and its target, enabling the ranking of compounds based on their docking scores and establishing a relationship between these scores and potential activity [32]. The visualization of interactions obtained from docking studies can further assist in improving the affinity characteristics of the ligands under investigation [33]. Consequently, a thorough molecular docking analysis was conducted to assess the binding affinity of the best cytotoxic compounds, compounds 8a, 8h, and 8i, towards VEGFR2.
The active site of VEGFR-2 features straightforward architecture consisting of a front pocket and a back pocket. The front pocket, responsible for ATP binding, is influenced by two critical residues, namely Glu917 and Cys919. On the other hand, the back pocket, which is hydrophobic in nature, includes Glu885 and Asp1046. Glu885 is positioned along the αC helix, while Asp1046 forms an integral component of the triad [34]. The compounds displayed a similar binding pattern withing the binding cavity highlighted by occupying the same position within the binding site. Moreover, all three compounds established two hydrogen bonds with Glu885 and Cys919 amino acid residues of the binding cavity. Moreover, the compounds shared a common π-sulfur interaction with Cys1045 amino acid residues ( Figure 4). Interestingly, in the crystallized ligand, Tivozanib established three hydrogen bonds with Glu885, Cys919, and Asp1046, which further validates the ability of the tested compounds as VEGFR2 inhibitors.

Molecular Docking
Molecular docking plays a significant role in computer-aided drug design and has diverse applications. One such application involves predicting the binding modes between a ligand and its target, enabling the ranking of compounds based on their docking scores and establishing a relationship between these scores and potential activity [32]. The visualization of interactions obtained from docking studies can further assist in improving the affinity characteristics of the ligands under investigation [33]. Consequently, a thorough molecular docking analysis was conducted to assess the binding affinity of the best cytotoxic compounds, compounds 8a, 8h, and 8i, towards VEGFR2.
The active site of VEGFR-2 features straightforward architecture consisting of a front pocket and a back pocket. The front pocket, responsible for ATP binding, is influenced by two critical residues, namely Glu917 and Cys919. On the other hand, the back pocket, which is hydrophobic in nature, includes Glu885 and Asp1046. Glu885 is positioned along the αC helix, while Asp1046 forms an integral component of the triad [34]. The compounds displayed a similar binding pa ern withing the binding cavity highlighted by occupying the same position within the binding site. Moreover, all three compounds established two hydrogen bonds with Glu885 and Cys919 amino acid residues of the binding cavity. Moreover, the compounds shared a common π-sulfur interaction with Cys1045 amino acid residues ( Figure 4). Interestingly, in the crystallized ligand, Tivozanib established three hydrogen bonds with Glu885, Cys919, and Asp1046, which further validates the ability of the tested compounds as VEGFR2 inhibitors.

MD Simulations
The synthesized compounds exhibited potent cytotoxic activity and demonstrated favorable interactions as revealed by molecular docking studies. However, it is essential to acknowledge that this method solely considers the flexibility of ligand conformations while keeping the protein structure rigid. Moreover, MD simulation plays a critical role in assessing theoretical accuracy, necessitating significant computational resources. Various theoretical methods are employed during MD analysis to select conformations from simulations, enabling the identification of relevant molecular conformations and aiding in the interpretation of simulation data [35,36]. To evaluate the stability of the binding poses and assess the dynamics of the protein conformation, molecular dynamics (MD) simulations were conducted for 100 ns on protein-ligand complexes involving compounds 8a, 8h, and 8i.
The stability of the protein-ligand complexes and the reliability of the simulation protocol were evaluated by calculating the root mean square deviation (RMSD). The RMSD values provide a quantitative measure of the structural similarity between the initial and final states of the protein-ligand complexes, enabling an assessment of the stability of these complexes throughout the simulation [37]. Moreover, by comparing the simulation results with the unbound state of the protein and a known complex, the validity of the simulation protocol can be assessed [38]. This approach offers a more accurate understanding of the binding interactions and the stability of the complexes, accounting for the dynamic behavior of the protein and the ligand in solution.  . Depiction of the 3D binding orientation of compounds 8a (orange), 8h (blue), and 8i (magenta) overlayed within the binding cavity of VEGFR2 and their corresponding 2D interactions established with the binding cavity's amino acid residues. Favorable interactions represented as dashed lines: green-hydrogen bonds; yellow-π-sulfur; dark pink-π-π stacking interactions; light pink-hydrophobic interactions; turquoise-halogen interaction.

MD Simulations
The synthesized compounds exhibited potent cytotoxic activity and demonstrated favorable interactions as revealed by molecular docking studies. However, it is essential to acknowledge that this method solely considers the flexibility of ligand conformations while keeping the protein structure rigid. Moreover, MD simulation plays a critical role in assessing theoretical accuracy, necessitating significant computational resources. Various theoretical methods are employed during MD analysis to select conformations from simulations, enabling the identification of relevant molecular conformations and aiding in the interpretation of simulation data [35,36]. To evaluate the stability of the binding poses and assess the dynamics of the protein conformation, molecular dynamics (MD) simulations were conducted for 100 ns on protein-ligand complexes involving compounds 8a, 8h, and 8i.
The stability of the protein-ligand complexes and the reliability of the simulation protocol were evaluated by calculating the root mean square deviation (RMSD). The RMSD values provide a quantitative measure of the structural similarity between the initial and final states of the protein-ligand complexes, enabling an assessment of the stability of these complexes throughout the simulation [37]. Moreover, by comparing the simulation results with the unbound state of the protein and a known complex, the validity of the simulation protocol can be assessed [38]. This approach offers a more accurate understanding of the binding interactions and the stability of the complexes, accounting for the dynamic behavior of the protein and the ligand in solution.
After initial fluctuations in the first 50 ns, the complexes of compounds 8a, 8h, and 8i displayed a similar stability profile with average RMSD of 2.1 Å, which is well below the accepted range of 3 Å, indicating the formation of stable complexes with VEGFR2. Moreover, the ligands in all three complexes showed RMSD of less than 1 Å (raw data file), Figure 4. Depiction of the 3D binding orientation of compounds 8a (orange), 8h (blue), and 8i (magenta) overlayed within the binding cavity of VEGFR2 and their corresponding 2D interactions established with the binding cavity's amino acid residues. Favorable interactions represented as dashed lines: green-hydrogen bonds; yellow-π-sulfur; dark pink-π-π stacking interactions; light pink-hydrophobic interactions; turquoise-halogen interaction.
After initial fluctuations in the first 50 ns, the complexes of compounds 8a, 8h, and 8i displayed a similar stability profile with average RMSD of 2.1 Å, which is well below the accepted range of 3 Å, indicating the formation of stable complexes with VEGFR2. Moreover, the ligands in all three complexes showed RMSD of less than 1 Å (raw data file), indicating stability throughout the simulation and a proper fit inside the binding cavity (Supporting Information, Figure S70).
The affinity of a molecule for a binding site is greatly influenced by its capacity to establish and sustain hydrogen bonds with residues within the binding site [39]. To further evaluate the stability of ligand binding modes, an examination of the intermolecular hydrogen bonds formed within the protein-ligand complexes was conducted. Hydrogen bond analysis revealed that only compound 8a was able to maintain a hydrogen bond with Cys919 amino acid residue of the hinge area, which is reported to be a crucial hydrogen bond formation with the adenine ring of ATP in the front pocket. The hydrogen bond between compound 8a and Cys919 was maintained for more than 50% of the simulation time, indicating its stability. Meanwhile, all three compounds were able to establish stable hydrogen bonds with Glu885 (located in the αC-helix) and Asp1046 (which is a part of the DFG motif) indicating a shared binding mode of action against the hydrophobic back pocket, which is formed by the DFG-out conformation (Supporting Information, Figure S71) [40][41][42][43]. Compound 8h maintained the hydrogen bonds with Glu885 and Asp1046 for around 60% of the simulation time. Meanwhile, compound 8i was able to maintain the same hydrogen bonds for almost 100% of the simulation time, indicating greater stability. The ability to consistently maintain stable hydrogen bonds throughout the simulations confirms that the compounds share a common binding mode. Additionally, it highlights the crucial role played by these particular amino acid residues in facilitating the binding of inhibitors to VEGFR2. This finding aligns with previous reports that have also emphasized the essentiality of these amino acids in the binding process [44,45].

MM-GBSA
The MM-GBSA (molecular mechanics-generalized Born surface area) method is a robust and widely utilized approach for predicting binding free energy (∆GBind) after conducting simulations, considering protein flexibility, entropy, and polarizability [46]. It effectively identifies ligands that exhibit efficient binding to receptors and plays a vital role in biomolecular research for comprehending molecular activities [47]. To validate the compounds identified using docking and molecular dynamics simulations, calculations of MM-GBSA binding free energy were employed. At intervals of 10 frames from frame 0 to 1001, post-simulation MM-GBSA was computed, resulting in a total of 100 conformations for each simulated complex. The MM-GBSA binding energy statistics provided in Table 2 indicate the average cumulative contributions of coulombic, hydrogen bonding, solvation energy, and lipophilic interactions, all of which exert a significant influence on ∆GBind. Out of the three complexes under investigation, it was observed that the complex formed by compound 8i exhibited the highest binding free energy (∆GBind) with a value of −59.55 kcal/mol. This superior binding affinity can be attributed to its notable ability to establish both hydrogen bonds and lipophilic interactions, surpassing the complexes formed by the other two compounds. Notably, throughout the entire duration of the 100% MD simulation, compound 8i consistently maintained its hydrogen bonds, as depicted in Figure S71 as shown in Supporting Information. This persistence in forming and preserving hydrogen bonds further reinforces the prediction that compound 8i exhibits the most potent activity towards VEGFR2. Moreover, compound 8i, which exhibited the highest binding free energy, was the most potent compound against the MCF-7 cell line, indicating the presence of a correlation between the free binding energy and the cytotoxicity of the compounds.

In Silico Pharmacokinetic Study
In the development of active therapeutic agents, pharmacokinetic properties such as absorption, metabolism, excretion, and toxicity (ADMET) play a critical role. It is important to note that a successful antagonistic interaction of inhibitors with a receptor protein or enzyme does not guarantee their suitability as drugs [48]. The presence of poor ADME characteristics and unfavorable toxicology are common reasons for the failure of drug candidates in clinical experiments [49]. Consequently, conducting ADME analysis is vital in the process of drug development. The freely accessible web server SwissADME was utilized to predict the pharmacokinetic properties of the most active compounds (compounds 8a, 8h, and 8i). This machine learning platform employs graph-based signatures encoded as distance/pharmacophore patterns to predict small-molecule pharmacokinetic properties [50].
According to the calculations, all three compounds (Table 3) were found to be watersoluble and demonstrated high gastrointestinal (GI) absorption. Additionally, they were found to comply with Lipinski's rule, suggesting their potential for oral activity. Furthermore, the predictions indicated that none of the three compounds were able to penetrate the blood-brain barrier (BBB), indicating a low likelihood of causing central nervous system (CNS) side effects.

Conclusions
In this study, rational drug design was employed to design a novel series of pyridineurea hybrid compounds. Using a stepwise optimization process, a group of 3-methoxy-1Hindazoles that possesses potent anticancer activity was identified. These novel compounds demonstrated significant efficacy against both breast (MCF7) and colon (HCT116) cancer cell lines. Among the synthesized compounds, compounds 8a, 8h, and 8i exhibited the highest cytotoxicity against both MCF7 and HCT116 cell lines. Compound 8a demonstrated exceptional potency against MCF7 cells, with a GI 50 value of 0.06 ± 0.014 µM. This GI 50 value indicates around sixfold greater potency compared to the positive control drug irinotecan. Furthermore, when tested on normal cells, compounds 8a, 8h, and 8i exhibited low cytotoxicity, indicating a promising safety profile. In silico investigation of the compounds predicted their action as VEGFR2 inhibitors and explained their mechanism of action. This information provides valuable insights into their mode of action and further supports their potential as anticancer agents.

Chemistry
General Procedures: The commercial chemicals, reagents, isocyanates, and solvents used were purchased from Sigma-Aldrich (Seoul, Republic of Korea), TCI (Seoul, Republic of Korea), Alfa Aesar (Seoul, Republic of Korea), and U CHEM (Seoul, Republic of Korea). All reactions were carried out in oven-dried glassware under a nitrogen atmosphere. Reactions were monitored using analytical thin-layer chromatography on 0.25 mm silica plates (E. Merck; silica gel 60 F254). Flash column chromatography was performed on silica gel 60 (230-400 mesh Kieselgel 60) using a Biotage system. The mass spectra were recorded using high-resolution mass spectrometry (HRMS) (electron ionization MS) on a Waters G2 QTOF mass spectrometer. The chemical shifts were reported in parts per million (ppm) and the coupling constants in hertz (Hz) (proton magnetic resonance spectra were recorded using a Varian 400 MHz spectrometer and carbon magnetic resonance spectra were recorded using a Varian 100 MHz spectrometer (Varian Medical Systems, Inc., Palo Alto, CA, USA)). RP-HPLC data analysis of the final products was performed using a Waters Corp. HPLC system equipped with an ultraviolet (UV) detector set at 254 nm and a YMC C18 (HS12S05-1564WT) column (5 µM, 12 nm) that was 4.6 mm × 150 mm in size was also employed. The RP-HPLC mobile phases used were (A) water containing 0.05% TFA and (B) acetonitrile. The purity of compounds was assessed using a gradient of 5% to 100% of acetonitrile in 40 min (method A) or a gradient of 5% to 60% of acetonitrile in 40 min (method B) with a flow rate of 1.0 mL/min. Melting points were measured using a Thermo Scientific 9200 melting point apparatus. For all biologically tested compounds, the purity was >95%. The spectral data for all synthesized compounds can be found in the Supplementary Materials.

SRB Assay of Antiproliferative Activity
The SRB assay was performed following the methodology of previously published studies [51]. In brief, the cells (100 µL medium containing 5000 HCT116 or 10,000 MCF7 cells/well) were incubated in 96-well microtiter plates. After 24 h, the drug was added (100 µL) to each well and the cultures were incubated for 48 h at 37 • C. The cells were fixed in 50% TCA (50 µL per well) for 1 h at 4 • C. The liquid was removed from the plate, which was then rinsed five times with water and allowed to dry at room temperature (RT). Washed cells were stained for 10 min at RT with 0.4% SRB (100 µL per well). After staining, the plate was washed 3 times with 1% glacial acetic acid and dried at RT. The SRB stain was then solubilized in 10 mM Tris and absorbances were read at 515 nm. The effect of the drugs was expressed in terms of the GI 50 (the concentration resulting in 50% maximal inhibition of cell proliferation). Each experiment was repeated at least three times, and the GI50 value calculated based on the results was expressed as mean ± standard deviation.

SRB Assay of Antiproliferative Activity on Normal Cells
The antiproliferative activity of compounds 8a, 8h, and 8i was evaluated in a sulforhodamine B (SRB) assay using normal cell line HEK293 cells. Cells were treated for 48 h under the same conditions used for the assessment of the antiproliferative activity on the cancer cells.

Molecular Docking
The crystal structure of VEGFR2 (juxtamembrane and kinase domains) in complex with Tivozanib (Supporting Information, Figure S69) was downloaded from the protein data bank (PDB ID: 4ASE; resolution: 1.83 Å) [52]. It is interesting to note that Tivozanib possessed a urea linker as well as a methoxy group on its "tail moiety". During the protein preparation phase, several actions were taken. Firstly, water molecules and crystallized inhibitor (Tivozanib) were removed, while any missing residues and hydrogen atoms were added. Subsequently, the prepared protein structure was protonated to attain a state suitable for physiological pH conditions. Finally, energy minimization was performed using the OPLS-2005 force field and a conjugate gradient method. This optimization process aimed to achieve a unique low-energy minimum for the protein structure. The Ligprep module of Maestro Schrodinger 2021.2 was implemented to prepare the proteins, involving the removal of water molecules and co-crystals as well as the inclusion of any missing residues or hydrogen atoms. The Glide extra precision module of Maestro Schrodinger was used for docking the compounds with each protein structure, producing 32 poses for each ligand. To establish the center of the Glide grid, the position of the co-crystallized ligand was utilized. Default settings were employed for both the grid generation and the docking processes. Using BIOVIA Discovery Studio Visualizer 2021 software, the pose with the lowest energy score was chosen for visualization [53,54].

MD Calculations
DESMOND software (version 2021.2) was utilized to perform molecular dynamics (MD) simulations with specific parameters. The simulation system was constructed by creating an orthorhombic box, which was solvated using single-point charge (SPC) water molecules. The temperature of the system was maintained at 300 K, and the pressure was set to 1 bar throughout the MD run. The simulation duration was set to 100 ns, and a relaxation time of 1 ps was applied for the selected positions. To minimize the energy of the solvated system, OPLS 2005 force parameters were employed. Electrostatic interactions were handled using the particle mesh Ewald (PME) method, while periodic boundary conditions (PBC) and a nonbond interaction cutoff distance of 9.0 Å (Angstrom) were implemented [53].
For each MD simulation, a six-step relaxation protocol was followed. This involved 2000 steps of LBFGS minimization, with the solute being restricted, and a loose convergence criteria of 50 kcal mol −1 Å −1 . Two preliminary simulations were conducted: a 12 ns simulation to equilibrate temperature (T = 10 K) with a thermostat relaxation constant of 0.1 ps, and a 12 ns simulation to equilibrate pressure (P = 1 atm) with a barostat relaxation constant of 50 ps. In both preliminary simulations, the nonhydrogen solute atoms were restrained. Subsequently, a 24 ps NPT ensemble simulation was performed to complete the preparation processes, with a temperature of 300 K, a thermostat relaxation constant of 0.1 ps, and a pressure of 1 atm. After the relaxation phase, the systems underwent a 5 ns MD simulation in the NPT ensemble. This simulation utilized a Nose-Hoover thermostat and Martyna-Tobias-Klein barostat, maintaining a temperature of 300 K with a thermostat relaxation time of 1.0 ps, a pressure of 1 atm, and a barostat relaxation time of 2.0 ps. The resulting MD simulation data were analyzed using QtGrace (version 2009) and Microsoft Excel [53].

Molecular Mechanics-Generalized Born Surface Area (MM-GBSA) Calculations
The MM-GBSA method, which combines molecular mechanics and the generalized Born surface area method, was used to calculate the binding free energy of ligand-protein complexes [55]. The thermal mmgbsa.py script was used to perform the calculations, which employed the 0-100 ns MD simulation trajectories, the VSGB solvation model, the OPLS3e force field, and a sampling rate of 10 steps per ns [56]. The binding free energy was calculated using the law of additivity, which took into account a variety of energy components such as hydrogen bonding, van der Waals interactions, columbic interactions, lipophilic interactions, covalent interactions, solvation, stacking, and self-contact of the ligand and protein [57].

In Silico Pharmacokinetic Study
The in silico pharmacokinetics of the most active compounds were predicted using the SwissADME server [58].