Exploring Pyrrolo-Fused Heterocycles as Promising Anticancer Agents: An Integrated Synthetic, Biological, and Computational Approach

Five new series of pyrrolo-fused heterocycles were designed through a scaffold hybridization strategy as analogs of the well-known microtubule inhibitor phenstatin. Compounds were synthesized using the 1,3-dipolar cycloaddition of cycloimmonium N-ylides to ethyl propiolate as a key step. Selected compounds were then evaluated for anticancer activity and ability to inhibit tubulin polymerization in vitro. Notably, pyrrolo[1,2-a]quinoline 10a was active on most tested cell lines, performing better than control phenstatin in several cases, most notably on renal cancer cell line A498 (GI50 27 nM), while inhibiting tubulin polymerization in vitro. In addition, this compound was predicted to have a promising ADMET profile. The molecular details of the interaction between compound 10a and tubulin were investigated through in silico docking experiments, followed by molecular dynamics simulations and configurational entropy calculations. Of note, we found that some of the initially predicted interactions from docking experiments were not stable during molecular dynamics simulations, but that configurational entropy loss was similar in all three cases. Our results suggest that for compound 10a, docking experiments alone are not sufficient for the adequate description of interaction details in terms of target binding, which makes subsequent scaffold optimization more difficult and ultimately hinders drug design. Taken together, these results could help shape novel potent antiproliferative compounds with pyrrolo-fused heterocyclic cores, especially from an in silico methodological perspective.


Introduction
Pyrrole and its hetero-fused derivatives have sparked great interest in the field of medicinal chemistry as valuable scaffolds for generating new drugs with biological activity, including anticancer, antimicrobial, or antiviral compounds [1]. Semisynthetic pyrroloquinolines are currently used as first-line agents in cancer therapy, including camptothecin analogs irinotecan and topotecan, with many others currently being investigated as potent antiproliferative drugs [2]. Pyrrolo [2,1-a]isoquinoline is a common motif in various bioactive alkaloids such as lamellarin D, which is known for its topoisomerase I inhibitory properties [3], or the antiproliferative alkaloid crispine A [4]. Pyrrolopyrimidine-based derivatives such as pemetrexed, immucillin H, or variolin B are also used as anticancer  Figure 1. Rational design strategy of pyrrolo-heterocyclic target compounds [6,11,14,35].
The tubulin polymerization inhibition mechanism can be partly responsible for the anticancer activity of these compounds [6,11,34,35]. The best compound in the series, pyrrolo[1,2-a]quinoline 10a, was further chosen for a series of in silico investigations, including global and local docking, molecular dynamics simulations, and configurational entropy calculations. This thorough in silico evaluation was conducted in order to describe the interaction of this compound with the colchicine binding site in detail, as well as investigate the stability of identified ligand-protein interactions through molecular docking experiments.
In the above context, the herein study was designed in order to identify novel anticancer tubulin-targeting agents that bind to the colchicine site, as well as describe the molecular details of hit binding through extensive in silico investigations. The broader goal of our investigation is to deepen the molecular understanding of biological activity and help shape future rational drug design campaigns aimed at the colchicine site of tubulin.

Chemistry
The strategy to build the five series of target phenstatin analogs ( Figure 1) consisted of two main synthetic steps, starting from the desired heterocycle to be fused to the pyrrole ring.
The tubulin polymerization inhibition mechanism can be partly responsible for the anticancer activity of these compounds [6,11,34,35]. The best compound in the series, pyrrolo[1,2-a]quinoline 10a, was further chosen for a series of in silico investigations, including global and local docking, molecular dynamics simulations, and configurational entropy calculations. This thorough in silico evaluation was conducted in order to describe the interaction of this compound with the colchicine binding site in detail, as well as investigate the stability of identified ligand-protein interactions through molecular docking experiments.
In the above context, the herein study was designed in order to identify novel anticancer tubulin-targeting agents that bind to the colchicine site, as well as describe the molecular details of hit binding through extensive in silico investigations. The broader goal of our investigation is to deepen the molecular understanding of biological activity and help shape future rational drug design campaigns aimed at the colchicine site of tubulin.

Chemistry
The strategy to build the five series of target phenstatin analogs ( Figure 1) consisted of two main synthetic steps, starting from the desired heterocycle to be fused to the pyrrole ring.
The fused target compounds were obtained in moderate yields ranging from 30% to 80%. The structures of all intermediate and final compounds were fully confirmed by spectral and physicochemical data (Supplementary Figures S1-S19). The structures of the cycloadducts presented in Schemes 1-5 were confirmed by NMR and IR spectra. These spectra are consistent with those of similar compounds previously reported in the literature [6,9], where clear structural evidence was emphasized, including X-ray diffraction.
We employed the 1,3-dipolar cycloaddition method to synthesize pyrrolopyrazine/ pyrrolopyrimidine/indolizine rings. In this method, we generated heterocyclic ylides in situ from the monoquaternary salts 3, 6, 9, 12, and 15 and reacted them with ethyl propiolate in a basic medium. The reaction schemes for each series are illustrated in Schemes 1-5. As expected, cycloadditions occurred with the selective formation of the regioisomers 4, 7, 10, and 13, which is in agreement with the electronic effects within ethyl propiolate and ylide species.
The fused target compounds were obtained in moderate yields ranging from 30% to 80%. The structures of all intermediate and final compounds were fully confirmed by spectral and physicochemical data (Supplementary Figures S1-S19). The structures of the cycloadducts presented in Schemes 1-5 were confirmed by NMR and IR spectra. These spectra are consistent with those of similar compounds previously reported in the literature [6,9], where clear structural evidence was emphasized, including X-ray diffraction.

Anticancer Activity
All forty synthesized compounds (monoquaternary salts and cycloadducts) were electronically submitted to the National Cancer Institute (NCI) platform, and twenty-five compounds (4b-d, 7b-c, 9a, 9d, 10a-d, 12a-d, 13a-d, 15a, 15d, and 16a-d) were selected for single dose (10 −5 M) screening against a panel of 60 human tumor cell lines, representing leukemia, melanoma, and cancers of the lung, colon, central nervous system, ovary, kidney, prostate, and breast [41]. Selected representative results for best hits and phenstatin are summarized in Table 1, all the other tested compounds being inactive against the NCI cancer cells at 10 −5 M (results provided by NCI for all tested compounds can be found in Supplementary Figures S20-S44).
Compound 10a, which contains a pyrrolo[1,2-a]quinoline scaffold, showed excellent growth inhibitory properties against almost all tested cell lines, with an average GP (growth percent-growth of treated cultures relative to untreated cultures) for all tested cell lines of 25.7% (74.3% GP inhibition). This compound exhibited the best inhibitory activity on the growth of leukemia cell lines, with an average GP inhibition of 87.33%, prostate cancer lines (78.0% GP inhibition), and breast cancer lines (76.4% GP inhibition). Compound 10a also had a cytotoxic effect on several cell lines, most notably on melanoma MDA-MB-435 cells, where 54.59% of cultured cells were killed at a 10 −5 M concentration. Even if phenstatin was superior in terms of average GP inhibition, compound 10a showed similar or even better inhibitory properties against several cell lines. Interestingly, the other three compounds from this series, 10b-10d, did not exert any relevant inhibitory activity, suggesting that  Compounds that fulfilled predetermined threshold inhibition criteria were then chosen by the NCI for the second testing step to determine GI 50 , TGI, and LC 50 parameters. Showing the most significant growth inhibition, compound 10a was the only one selected for evaluation against the sixty tumor cell lines in five-dose assays [41][42][43]. The most representative results from the NCI-60 five-dose screen are shown in Table 2 (full results can be found in Figure S45 from the Supplementary Materials).  (Table 2). However, LC 50 values (the molar concentration of tested compound causing 50% death of tumor cells) for compound 10a and phenstatin, respectively, were found to be >100 µM against all tested lines.
NCI's in silico platform, PRISM, enables investigators to compare their active molecules with other molecules that have been previously tested through NCI from both a biological and chemical perspective [44]; this evaluation is performed with two components: COMPARE and PILOT. Interestingly, when analyzing the most active compound 10a with PRISM, we found a best-fitting profile with N-(4-bromophenyl)-4-(2-cyclopropyloxazol-5yl)benzenesulfonamide [45], which is a compound from a series of reported anticancer derivatives and is also able to inhibit tubulin polymerization, as we expect for our compounds.

In Vitro Tubulin Polymerization Inhibition
To confirm that the observed anticancer activity of compound 10a is conferred by a microtubule-targeting mechanism, we evaluated its effect on the assembly of tubulin, using paclitaxel (a tubulin stabilizer) and phenstatin (a tubulin polymerization inhibitor) as controls. As expected and shown in Figure 2, paclitaxel was found to stimulate tubulin polymerization, while phenstatin and compound 10a inhibited tubulin polymerization, although 10a had a slightly weaker effect. The obtained data indicate that compound 10a effectively inhibits tubulin polymerization in vitro, suggesting that the observed cytotoxicity of this compound is related to a microtubule-targeting mechanism. tion range. This pyrrolo[1,2-a]quinoline derivative outperformed control phens several cases, including renal cancer cell line A498, colon cancer cell line COLO 2 breast cancer cell line T-47D. Furthermore, compound 10a demonstrated good mance against various cancer cell lines (including melanoma MDA-MB-435, rena A498, ovarian cancer OVCAR-3, breast cancer MDA-MB-468, and non-small cell lu cer NCI-H522 cell lines), exhibiting complete growth inhibition at submicromolar trations (Table 2). However, LC50 values (the molar concentration of tested com causing 50% death of tumor cells) for compound 10a and phenstatin, respectivel found to be >100 µM against all tested lines. NCI s in silico platform, PRISM, enables investigators to compare their activ cules with other molecules that have been previously tested through NCI from biological and chemical perspective [44]; this evaluation is performed with two nents: COMPARE and PILOT. Interestingly, when analyzing the most active com 10a with PRISM, we found a best-fitting profile with N-(4-bromophenyl)-4-(2-cyc yloxazol-5-yl)benzenesulfonamide [45], which is a compound from a series of r anticancer derivatives and is also able to inhibit tubulin polymerization, as we ex our compounds.

In Vitro Tubulin Polymerization Inhibition
To confirm that the observed anticancer activity of compound 10a is conferr microtubule-targeting mechanism, we evaluated its effect on the assembly of tubu ing paclitaxel (a tubulin stabilizer) and phenstatin (a tubulin polymerization inhib controls. As expected and shown in Figure 2, paclitaxel was found to stimulate polymerization, while phenstatin and compound 10a inhibited tubulin polymer although 10a had a slightly weaker effect. The obtained data indicate that compou effectively inhibits tubulin polymerization in vitro, suggesting that the observed c icity of this compound is related to a microtubule-targeting mechanism.

Blind Docking
Since compound 10a is similar to phenstatin in inhibiting tubulin polymerization in vitro, we hypothesize that it binds to the colchicine site of the α,β-tubulin heterodimer, the known site for phenstatin [46]. In addition, recent tubulin crystal structures cocrystalized with quinazolinone and tetrahydroisoquinoline anticancer agents reveal the preference of these heterocyclic compounds to the colchicine binding site [47,48] and support our previous hypothesis that quinoline derivatives could exert their anticancer activity by binding to the colchicine site of tubulin [6]. However, tubulin binding sites other than colchicine have been described, including the vinca, taxol, and peloruside/laulimalide Pharmaceuticals 2023, 16, 865 9 of 30 sites [49], and other novel sites are being actively discovered and/or investigated for the rational design of tubulin modulators [50,51] or compounds targeting drug resistant cancer phenotypes [19,52]. Thus, we performed blind docking experiments for compound 10a, in order to investigate its relative preference towards the colchicine binding site, while validating our docking protocol using colchicine and phenstatin ( Figure S46 from Supplementary Material).
The generated poses for compound 10a had estimated binding energies between −8.2 and −2.0 kcal/mol, with an overall score distribution and RMSD from best-scoring conformation for all docked poses resembling both colchicine and phenstatin (Supplementary Figure S47). Clustering with a 2.0 Å RMSD tolerance gave 127 unique cluster representatives for 10a, 13 of which scored lower than −7 kcal/mol. Out of these, the lowest-scoring conformation and an additional 10 cluster representatives were positioned in the colchicine binding site. All colchicine cluster representatives with estimated binding energies lower than −7 kcal/mol were positioned in the colchicine binding site of tubulin, as expected. RMSD between co-crystallized colchicine and lowest-scoring cluster representative from the blind docking was 1.082 Å, which is lower than the 2.0 Å cutoff generally used as a criterion for correct bound structure prediction [53], suggesting that the used docking protocol is suitable for the studied system. None of the phenstatin conformations scored less than −7 kcal/mol, but the four lowest-scoring cluster representatives were also positioned in the colchicine binding site. These results suggest that compound 10a prefers the colchicine binding site of tubulin and is likely able to inhibit tubulin polymerization by binding to the colchicine binding site of tubulin while having low affinity for other tubulin sites.

Local Docking
As the blind docking experiments revealed that the most energetically favorable poses for compound 10a are positioned in the colchicine binding site of tubulin, we performed local docking on this site in order to investigate the molecular nature of these preferential conformations. A clear improvement in the overall docking score distribution can be observed when compared to blind docking due to the decrease in conformational search space, as expected (Supplementary Figure S47). Interestingly, three low scoring clusters of comparable energies (−8.95 kcal/mol, −8.77 kcal/mol, and −8.69 kcal/mol-conformations I, II, and III, respectively) were revealed through local docking and were accommodated in the colchicine binding site in three geometrically different modes (BM I, BM II, and BM III). Of note, BM III greatly overlapped with the lowest-scoring cluster representative from blind docking for compound 10a (RMSD 1.034 Å). These three lowest-scoring cluster representatives were further chosen for analysis, as there have been cases of false positives in top-scoring poses from docking experiments, and subsequent analysis of multiple low-scoring docking poses is strongly recommended to correctly identify the most likely ligand conformation [54,55]. The RMSD between these modes was 7.51 Å and 6.408 Å (using BM I as a reference, as it is the lowest-scoring cluster), and the distances between their centers of mass were 3.197 Å (BM I/BM II), 3.294 Å (BM I/BM III), and 5.145 Å (BM II/BM III), indicating that the three conformations greatly differ from one another. The three best-scoring docking solutions from local docking also varied in terms of key amino acids involved in ligand stabilization at the binding site ( Figure 3). For docking protocol validation, colchicine was re-docked in the same manner as 10a. RMSD between co-crystallized colchicine and the lowest-scoring cluster representative after colchicine re-docking was 1.110 Å.
pocket lined by βCys241, βLeu242, βLeu252, βLeu255, βMet259, βVal315, and βAla316, while the trimethoxy-substituted ring is positioned towards the α subunit, in a manner similar to what we have previously seen for moderately active isoquinoline cycloadducts [6], yet very dissimilar to colchicine (Supplementary Figure S48). This conformation is further stabilized through hydrogen bonding with the sidechain of βLys254, an amino acid that has been previously identified through molecular docking experiments as an interaction partner for other tubulin polymerization inhibitors with anticancer activity [56][57][58]. However, to our knowledge, none of the co-crystallized tubulin inhibitors that bind to the colchicine site have been shown to interact with this residue, although it has been highlighted that this amino acid is likely involved in microtubule assembly through interaction with the γ-phosphate group of the N-site GTP [59]. Furthermore, alanine scan mutations have shown that D251A-R253A-K254A is lethal in yeast [60]. In addition, a recent molecular dynamics study focused on the dynamics of the βT7 loop concluded that an interaction between the backbone N atom of βLys254 and the side chain O atom of βAsp251 is essential for the structural rearrangement of the βT7 loop [61], which has been shown to play an important role in colchicine binding [62] and prevents the 'curved-tostraight structural transition of tubulin from its free form, a process that is necessary for microtubule formation [63]. Therefore, βLys254 could be a likely interaction partner for compound 10a.
BM II was stabilized in the colchicine binding site exclusively through hydrophobic interactions, having an orientation similar to what we have previously predicted for active isoquinoline derivatives [6]. However, no hydrogen bond interaction with βCys241, the presumed anchor point for the trimethoxy-substituted ring, was observed. Nevertheless, the trimethoxy-substituted ring roughly occupies the same hydrophobic pocket as the trimethoxy-substituted moiety of colchicine ( Figure S48). Since some colchicine site binders, including 2-aroylindoles, have been shown to be stabilized in the colchicine site exclusively through hydrophobic interactions and water-mediated polar interactions that do not include βCys241 [64,65], such a conformation could also be probable for compound Residues are annotated with a three-letter amino acid code and corresponding chain and colored as follows: hydrophobic-green; polar uncharged-light blue; negatively charged-red; and positively charged-dark blue; hydrogen bonds are indicated as pink arrows.
In the case of BM I, the pyrrolo[1,2-a]quinoline ring is buried in the hydrophobic pocket lined by βCys241, βLeu242, βLeu252, βLeu255, βMet259, βVal315, and βAla316, while the trimethoxy-substituted ring is positioned towards the α subunit, in a manner similar to what we have previously seen for moderately active isoquinoline cycloadducts [6], yet very dissimilar to colchicine (Supplementary Figure S48). This conformation is further stabilized through hydrogen bonding with the sidechain of βLys254, an amino acid that has been previously identified through molecular docking experiments as an interaction partner for other tubulin polymerization inhibitors with anticancer activity [56][57][58].
However, to our knowledge, none of the co-crystallized tubulin inhibitors that bind to the colchicine site have been shown to interact with this residue, although it has been highlighted that this amino acid is likely involved in microtubule assembly through interaction with the γ-phosphate group of the N-site GTP [59]. Furthermore, alanine scan mutations have shown that D251A-R253A-K254A is lethal in yeast [60]. In addition, a recent molecular dynamics study focused on the dynamics of the βT7 loop concluded that an interaction between the backbone N atom of βLys254 and the side chain O atom of βAsp251 is essential for the structural rearrangement of the βT7 loop [61], which has been shown to play an important role in colchicine binding [62] and prevents the 'curved-to-straight' structural transition of tubulin from its free form, a process that is necessary for microtubule formation [63]. Therefore, βLys254 could be a likely interaction partner for compound 10a.
BM II was stabilized in the colchicine binding site exclusively through hydrophobic interactions, having an orientation similar to what we have previously predicted for active isoquinoline derivatives [6]. However, no hydrogen bond interaction with βCys241, the presumed anchor point for the trimethoxy-substituted ring, was observed. Nevertheless, the trimethoxy-substituted ring roughly occupies the same hydrophobic pocket as the trimethoxy-substituted moiety of colchicine ( Figure S48). Since some colchicine site binders, including 2-aroylindoles, have been shown to be stabilized in the colchicine site exclusively through hydrophobic interactions and water-mediated polar interactions that do not include βCys241 [64,65], such a conformation could also be probable for compound 10a. In fact, this conformation occupies hydrophobic center II according to a previously described structure-based colchicine binding-site inhibitor model [21] and greatly overlaps with 2-aroylindole tubulin polymerization inhibitor D64131 [65], with an additional hydrophobic extension towards the α subunit ( Figure S49).
In BM III, the trimethoxy-substituted ring was oriented towards the α subunit and engaged in hydrogen bonding with many amino acid sidechains known or thought to be involved in colchicine binding site inhibitor interaction, including αAsn101 [51], αThr179 [65], αSer178 [31], βLys352 [66], and βGln247 [67]. However, to our knowledge, no other colchicine site inhibitors possessing a trimethoxy-substituted ring have been co-crystallized in a similar conformation [18]. This conformation is also the most dissimilar from colchicine ( Figure S48).
Importantly, all three lowest cluster representatives of 10a had similar theoretical binding energies. While scoring functions take into account a wide range of contributions, including electrostatic interactions, van der Waals contacts, and desolvation effects [68], they are limited in reflecting the conformational changes induced by ligand binding, as well as the stability of identified amino acid interactions in time.

Molecular Dynamics Simulations
To investigate the detailed dynamics and interaction stability of ligand-tubulin complexes, the tubulin heterodimer and the three best-scoring local docking poses were subjected to 10 ns MD simulations.
Overall, all three systems reached equilibrium, according to root mean square deviation (RMSD) analysis (Figure 4a), even though RMSD was higher for all investigated binding modes when compared to reference colchicine. In addition, BM I shifted from the initial frame more than BM II or BM III (Figure 4b), having a mean RMSD of 1.154 Å, while BM II and BM III had mean RMSD of 0.721 Å and 0.66 Å, respectively. All binding modes had a higher RMSD than colchicine (mean deviation 0.39 Å), which would be expected since the starting tubulin heterodimer structure is co-crystallized with colchicine. The 'Lig fit Prot' (Figure 4c) represents the RMSD of ligand-heavy atoms after first aligning the protein-ligand complex on the protein backbone. If the observed values differ greatly from the protein backbone RMSD, it is likely that the ligand diffuses away from the binding site [69]. In the case of BM I, this parameter shows drift during the simulation, particularly at the beginning of the trajectory, but is stabilized until the end of the simulation. This suggests that the initially found local docking conformation does not maintain close contact with surrounding amino acids and that another conformation is stabilized during the trajectory.
10a. In fact, this conformation occupies hydrophobic center II according to a previously described structure-based colchicine binding-site inhibitor model [21] and greatly overlaps with 2-aroylindole tubulin polymerization inhibitor D64131 [65], with an additional hydrophobic extension towards the α subunit ( Figure S49).
In BM III, the trimethoxy-substituted ring was oriented towards the α subunit and engaged in hydrogen bonding with many amino acid sidechains known or thought to be involved in colchicine binding site inhibitor interaction, including αAsn101 [51], αThr179 [65], αSer178 [31], βLys352 [66], and βGln247 [67]. However, to our knowledge, no other colchicine site inhibitors possessing a trimethoxy-substituted ring have been co-crystallized in a similar conformation [18]. This conformation is also the most dissimilar from colchicine ( Figure S48).
Importantly, all three lowest cluster representatives of 10a had similar theoretical binding energies. While scoring functions take into account a wide range of contributions, including electrostatic interactions, van der Waals contacts, and desolvation effects [68], they are limited in reflecting the conformational changes induced by ligand binding, as well as the stability of identified amino acid interactions in time.

Molecular Dynamics Simulations
To investigate the detailed dynamics and interaction stability of ligand-tubulin complexes, the tubulin heterodimer and the three best-scoring local docking poses were subjected to 10 ns MD simulations.
Overall, all three systems reached equilibrium, according to root mean square deviation (RMSD) analysis (Figure 4a), even though RMSD was higher for all investigated binding modes when compared to reference colchicine. In addition, BM I shifted from the initial frame more than BM II or BM III (Figure 4b), having a mean RMSD of 1.154 Å, while BM II and BM III had mean RMSD of 0.721 Å and 0.66 Å, respectively. All binding modes had a higher RMSD than colchicine (mean deviation 0.39 Å), which would be expected since the starting tubulin heterodimer structure is co-crystallized with colchicine. The 'Lig fit Prot (Figure 4c) represents the RMSD of ligand-heavy atoms after first aligning the protein-ligand complex on the protein backbone. If the observed values differ greatly from the protein backbone RMSD, it is likely that the ligand diffuses away from the binding site [69]. In the case of BM I, this parameter shows drift during the simulation, particularly at the beginning of the trajectory, but is stabilized until the end of the simulation. This suggests that the initially found local docking conformation does not maintain close contact with surrounding amino acids and that another conformation is stabilized during the trajectory.   Indeed, for BM I, the 10 ns MD simulation revealed that the interaction with βLys254, which was the main anchor point for this docking conformation, was only transient, being maintained in less than 1% of the simulation time (Figure 5a). However, BM I maintained contacts in more than 25% of simulation time with αTyr224 (44.96%), βAsn249 (58%), βLys352 (53.35%), βAla250 (59.44%), βLeu255 (61.49%), βAla316 (37.91%), and the cocrystalized GTP molecule (39.51%). The ligand RMSD with regards to the initial docking conformation had the highest value of all simulated binding modes (3.543 Å). Superimposition between the docked solution and the final frame of the simulation revealed major differences in terms of amino acid contacts within the colchicine binding site ( Figure S50). Of note, the carboxylate moiety flips after 3.57 ns and remains flipped throughout the simulation. The flexibility of this group is also evident in the ligand root mean square fluctuation-RMSF ( Figure S51). Overall, the BM I-containing simulated system converged to the least similar conformation from the starting docking solution.
differences in terms of amino acid contacts within the colchicine binding site ( Figure S50). Of note, the carboxylate moiety flips after 3.57 ns and remains flipped throughout the simulation. The flexibility of this group is also evident in the ligand root mean square fluctuation-RMSF ( Figure S51). Overall, the BM I-containing simulated system converged to the least similar conformation from the starting docking solution.
BM II maintained all the initially identified hydrophobic contacts, while slightly drifting from the initial docked conformation (RMSD 1.746 Å). In addition, a cation-pi interaction with βLys352 was observed in 27.32% of the simulation time. This binding mode was also the only one to engage in favorable interactions with βCys241, maintaining contact with this amino acid in 20.53% of simulation time (Figure 5b), mainly through the central methoxy moiety from its trimethoxy-substituted ring. The flexibility of this substituent in BM II is most evident from the abrupt increase or decrease in the main RMSD profile with respect to the first frame of the simulation (Figure 4b), as well as from the ligand root mean square fluctuation ( Figure S51). Of note, RMSF of the same substituent in colchicine shows a similar behavior throughout the simulation ( Figure S51). In the case of BM III, the H-bonds with αAsn101, αSer178, and βGln247 identified from the docking solution were mostly not maintained during the MD simulation (Figure 5c), except for the initial frames ( Figure S52). The H-bond with αThr179 was maintained in 15.03% of simulated time, either directly or mediated through a water molecule, while interaction with βLys352 was maintained throughout the entire simulation (Figure 5c). RMSD between the initial docked solution and the final frame of the MD simulation was 3.525 Å.
It is important to note that from all final frame MD simulation conformations, BM II was the only one that largely overlapped with co-crystallized colchicine binding site inhibitors with a trimethoxy-substituted ring [25,65,70], with most hydrophobic contacts being maintained throughout the simulation, even though the number of polar contacts was the smallest from all simulations ( Figures S53 and S54). In addition, all simulated BMs for compound 10a maintained various kinds of contacts with βLys352 in more than 25% of the simulated time, suggesting that this amino acid could be relevant for the binding of BM II maintained all the initially identified hydrophobic contacts, while slightly drifting from the initial docked conformation (RMSD 1.746 Å). In addition, a cation-pi interaction with βLys352 was observed in 27.32% of the simulation time. This binding mode was also the only one to engage in favorable interactions with βCys241, maintaining contact with this amino acid in 20.53% of simulation time (Figure 5b), mainly through the central methoxy moiety from its trimethoxy-substituted ring. The flexibility of this substituent in BM II is most evident from the abrupt increase or decrease in the main RMSD profile with respect to the first frame of the simulation (Figure 4b), as well as from the ligand root mean square fluctuation ( Figure S51). Of note, RMSF of the same substituent in colchicine shows a similar behavior throughout the simulation ( Figure S51).
In the case of BM III, the H-bonds with αAsn101, αSer178, and βGln247 identified from the docking solution were mostly not maintained during the MD simulation (Figure 5c), except for the initial frames ( Figure S52). The H-bond with αThr179 was maintained in 15.03% of simulated time, either directly or mediated through a water molecule, while interaction with βLys352 was maintained throughout the entire simulation (Figure 5c). RMSD between the initial docked solution and the final frame of the MD simulation was 3.525 Å.
It is important to note that from all final frame MD simulation conformations, BM II was the only one that largely overlapped with co-crystallized colchicine binding site inhibitors with a trimethoxy-substituted ring [25,65,70], with most hydrophobic contacts being maintained throughout the simulation, even though the number of polar contacts was the smallest from all simulations (Figures S53 and S54). In addition, all simulated BMs for compound 10a maintained various kinds of contacts with βLys352 in more than 25% of the simulated time, suggesting that this amino acid could be relevant for the binding of this compound in the colchicine site, regardless of the initially identified docking solution. The importance of this residue in 10a binding could be further investigated through alanine scanning, as has been previously performed for the anticancer natural product pironetin [71].
Since receptor-ligand binding does not rely solely upon specific interactions between amino acids and particular chemical moieties but also on reduction of molecular flexibility upon binding [72], we further investigated the configurational entropy changes of all three binding modes upon interacting with the colchicine binding site of tubulin. This gave us an estimation of the entropic cost upon target binding. Since longer MD simulations are preferred in order to fully assess the stability of the investigated systems [30], and sufficient sampling should be reached to estimate configurational entropies [73], we extended our simulations to 25 ns, for which a good convergence of the configurational entropy profile was achieved. We found that the internal configurational entropy is reduced upon binding from 156.02 J/mol K to 99.52 J/mol K, 99.55 J/mol K, and 98.28 J/mol K, for BM I, BM II, and BM III, respectively (∆S conf = −56.50 J/mol K; −56.47 J/mol K; −57.74 J/mol K, respectively). Of note, BM II exhibited a markedly lower internal configurational entropy than BM I and BM III in the first quarter of the simulation but steadily increased to comparable values as the other two binding modes throughout the rest of the simulated time ( Figure S55). Our calculated values are similar in magnitude to other receptor-ligand systems [73].
Overall, the molecular modeling outcomes enabled us to predict, from a structural perspective, the interactions between our novel compound and tubulin with atomic resolution. Our results suggest that all three binding modes to tubulin described here are characterized by both positive enthalpic and entropic contributions, highlighting that all three binding modes are equally probable in a biological context. The molecular mechanics energies combined with the Poisson-Boltzmann or generalized Born and surface area continuum solvation (MM/PBSA and MM/GBSA) methods [74] could also have been of particular use in assessing ligand binding affinities for each of the three poses, but they have not been explored in the current study.

In Silico ADME and Toxicity Predictions
The result of the predicted parameters, including molecular properties, pharmacokinetics, drug-likeness, and medicinal chemistry are presented in Table 3. Compound 10a follows both Lipinski's rule of five (without any violations) and Veber's rule, having less than 10 rotatable bonds and a topological polar surface area smaller than 140 Å 2 . In addition, compound 10a displays a moderately soluble behavior, with a Log S (ESOL) value of −6.00 and no PAINS or Brenk alerts in its structure. Derivative 10a also scores well in terms of bioavailability and ease of synthetic accessibility (Table 3). This compound is predicted to have a high gastrointestinal absorption but is unable to permeate through the blood-brain barrier (BBB). Compound 10a does not appear to be a P-glycoprotein (P-gp) substrate.
The chart generated from the Swiss ADME QSAR web tool, regarding the accessibility of compound 10a to be orally bioavailable, is presented in Figure 6. Compound 10a follows both Lipinski s rule of five (without any violations) and Veber s rule, having less than 10 rotatable bonds and a topological polar surface area smaller than 140 Å 2 . In addition, compound 10a displays a moderately soluble behavior, with a Log S (ESOL) value of −6.00 and no PAINS or Brenk alerts in its structure. Derivative 10a also scores well in terms of bioavailability and ease of synthetic accessibility (Table 3). This compound is predicted to have a high gastrointestinal absorption but is unable to permeate through the blood-brain barrier (BBB). Compound 10a does not appear to be a P-glycoprotein (P-gp) substrate.
The chart generated from the Swiss ADME QSAR web tool, regarding the accessibility of compound 10a to be orally bioavailable, is presented in Figure 6. This radar involves six parameters, lipophilicity (LIPO), size, polarity (POLAR), insolubility (INSOLU), unsaturation (INSATU), and flexibility (FLEX), of the tested compound and is represented by a red line integrated into a pink area. Molecules that fall within the pink region of the radar are considered drug-like. Compound 10a exhibited compliance to only four of the six rules, with violations to INSATU (ratio of hybridized sp3 atoms to the total number of C atoms) and LIPO (XLOGP3 between −0.7 and +5.0). Taken together, these predicted data show a promising ADME and drug-likeness profile for compound 10a.
The predicted toxicity spectrum is represented by a list of activities with probabilities "to be active" (Pa) and "to be inactive" (Pi). The obtained results presented in Table 4 show predicted cytotoxicity (Pa > Pi and Pa > 0.3) against several cancer cell lines, including four of the tested NCI cell lines: T-47D, HT-29, DU-145, and MCF7 (Table 3). The fact that no normal human cell lines appeared on the list could be an indication for a good selectivity of compound 10a against cancer cell lines.  This radar involves six parameters, lipophilicity (LIPO), size, polarity (POLAR), insolubility (INSOLU), unsaturation (INSATU), and flexibility (FLEX), of the tested compound and is represented by a red line integrated into a pink area. Molecules that fall within the pink region of the radar are considered drug-like. Compound 10a exhibited compliance to only four of the six rules, with violations to INSATU (ratio of hybridized sp3 atoms to the total number of C atoms) and LIPO (XLOGP3 between −0.7 and +5.0). Taken together, these predicted data show a promising ADME and drug-likeness profile for compound 10a.
The predicted toxicity spectrum is represented by a list of activities with probabilities "to be active" (P a ) and "to be inactive" (P i ). The obtained results presented in Table 4 show predicted cytotoxicity (Pa > Pi and Pa > 0.3) against several cancer cell lines, including four of the tested NCI cell lines: T-47D, HT-29, DU-145, and MCF7 (Table 3). The fact that no normal human cell lines appeared on the list could be an indication for a good selectivity of compound 10a against cancer cell lines.

Chemistry
All commercially available reagents and solvents employed were used without further purification. Melting points were recorded on an A. Krüss Optronic Melting Point Meter KSP1 and are uncorrected. Analytical thin-layer chromatography was performed with commercial silica gel plates 60 F254 (Merck Darmstadt, Germany) and visualized under UV light (λ max = 254 or 365 nm). The NMR spectra were recorded on a Bruker Avance III 500 MHz spectrometer or a Bruker Avance 400 DRX (400 MHz) (Bruker, Vienna, Austria). Chemical shifts are reported in delta (δ) units, part per million (ppm), and coupling constants (J) in Hz. The following abbreviations are used to designate chemical shift multiplicities: s = singlet, d = doublet, t = triplet, q = quartet, m = multiplet, and bs = broad singlet. Infrared (IR) spectra were recorded as films on potassium bromide (KBr) pellets on an FTIR Prestige 8400 s spectrophotometer (Shimadzu, Kyoto, Japan) or a Jasco 660 FTIR spectrophotometer. Elemental analyses indicated by the symbols of the elements were within ± 0.4% of the theoretical values. HR-MS experiments were recorded on a HESI Orbitrap Exploris 120 Mass Spectrometer (Thermo Fisher, Walthan, MA, USA) in positive mode. Monoquaternary Salts 3, 6, 9, 12, and 15 The corresponding heterocycle (pyrazine 1, 4-(4-chlorophenyl)pyrimidine 5, quinoline 8, benzo[f ]quinoline 11 or isoquinoline 14) (1 mmol, 1 equiv.) was dissolved in 5-7 mL acetone. Then, the corresponding reactive 2-bromoacetophenone 2a-d (1.1 mmol, 1.1 equiv.) was added and the resulting mixture was stirred overnight at room temperature (r.t.). The formed precipitate was filtered and washed with diethyl ether to obtain the desired product, which was used in the next reaction without any further purification.

General Procedure for Compounds 4, 7, 10, 13, and 16
The cycloimmonium salt (3, 6, 9, 12, and 15) (1 mmol, 1 equiv.) and ethyl propiolate (1.1 mmol, 1.1 equiv.) were added to dichloromethane (DCM). Then, a solution of triethylamine (TEA) (3 mmol, 3 equiv.) in DCM (3 mL) was added dropwise over 1 h (magnetic stirring), and the resulting mixture was stirred overnight at r.t. Methanol (10 mL) was added, and the resulting solid was collected by filtration and washed with 5 mL methanol. The product was then purified by crystallization from dichloromethane/methanol (1/1, v/v) and/or column chromatography using dichloromethane/methanol (99.5/0.5, v/v). Compounds 4, 7, 10, 13, and 16 were the only pure compounds obtained from the reaction mixture fractions. While spectral evidence for some decomposition products of the ylides was observed in the other fractions, no evidence was found for the presence of the other possible regioisomer cycloadducts.   13  was initiated by the addition of tubulin to a final concentration of 3.0 mg/mL. Paclitaxel and phenstatin were used as positive controls under the same conditions. Absorbance was measured at 340 nm for 1 h at 1 min intervals at 37 • C.

Molecular Docking
Flexible-ligand docking experiments were performed in two steps. First, blind docking was performed using Autodock 4.2 [77] using a 90 × 118 × 114 gridbox with 0.825 Å spacing (total search volume of 998,811 Å 3 , as to include the entire heterodimer, PDB ID: 4O2B [63]) and was centered on the entire structure (x = 18.997, y = 68.705, z = 46.491). Co-crystallized GTP, GDP, Ca 2+ , and Mg 2+ ions were kept during receptor preparation, and the target protein was kept rigid during all docking experiments. The 3D structures of colchicine, phenstatin, and compound 10a were constructed in Avogadro v1.2.0 [78] and were energetically optimized in the MMFF94 force field until a local energy minimum was reached. A total of 8 jobs of 2000 runs were performed for each of the 3 ligands using the Lamarckian Genetic Algorithm (LGA) for a total of 16,000 runs per ligand. Conformations having theoretical binding energies lower than −5.5 kcal/mol were clustered and ranked using an RMSD tolerance of 2.0 Å in order to select the lowest-binding representatives for further binding pose refinement. Cluster representatives with binding energies lower than −7.0 kcal/mol were visually inspected for binding orientation.
A smaller gridbox centered on the lowest binding cluster of compound 10a from blind docking (58 × 58 × 58 points, 0.375 Å spacing, x = 14.316, y = 67.184, z = 44.464) was chosen for local docking. In this case, 1000 poses were generated per ligand in a single run, and results were clustered and ranked based on theoretical binding energy. RMSD between re-docked and co-crystallized colchicine ligand was used as quality control for the docking protocol. Visual inspection and molecular graphics were made in the PyMOL Molecular Graphics System, Version 2.5.0. (Schrödinger, LLC, New York, NY, USA). Furthermore, 2D ligand interaction diagrams were generated in Maestro, release 2018-4 (Schrödinger, LLC, New York, NY, USA).

Molecular Dynamics Simulations
Each of the three lowest-scoring cluster representatives of compound 10a obtained from local docking experiments was subjected to MD simulations in the colchicine binding site of the α,β-tubulin heterodimer (chains A and B of PDB ID: 4O2B). The BioLuminate ® (Schrödinger, Inc., New York, NY, USA) graphical interface [79] was used for initial system construction, while the subsequent MD simulations were performed with the Desmond (D.E. Shaw Research) software [80]. The missing loop of tubulin chain B (residues 276-281) was modeled using MODELLER [81], as implemented in UCSF Chimera release 1.15 [82]. Proteins, ligands, and ions were described using the OPLS-AA force field [83,84]. The 10a/tubulin complexes were solvated using a TIP3P water model [85] in orthorhombic simulation boxes while allowing for a minimum distance of 10 Å between the complex and the walls of the simulation box. Ions were added to electrostatically neutralize the systems and to attain a salt concentration of 0.15 M NaCl. After initial geometry optimization, the standard Desmond solute relaxation protocol was applied. This protocol involves sequential short MD simulations using restraints of decreasing magnitude on atomic positions, followed by short restraint-free equilibration simulations. In the production runs, each complex was simulated for 10 ns, and system configurations were saved every 5 ps for analysis. Throughout the simulations, the NPT ensemble (i.e., constant number of particles, pressure, and temperature) was employed, using the Nosé-Hoover chain thermostat (300 K) [86] and Martyna-Tobias-Klein barostat (1 atm) [87]. The MD trajectories were analyzed with respect to protein backbone and ligand RMSD, as well as proteinligand contacts. Protein-ligand interaction fraction plots were generated in R using the ggplot2 package [88]. Configurational entropy variations were evaluated from 25 ns long MD simulations of the free compound in solution and bound to tubulin, respectively.
Schlitter's formula was used for the configurational entropy calculation [89]. Prior to entropy calculations using the GROMACS 2020.3 modelling suite [90], all the structures in the MD trajectories were geometrically fitted to the first frame in order to remove the rotational and translational degrees of freedom. The fitting procedure included all heavy atoms of compound 10a. To evaluate the convergence of the entropy profile, the calculations included steeply increasing time intervals along the entire MD trajectories.

In Silico ADME and Toxicity Predictions
The ADME in silico evaluation for the most active compound 10a was performed using the Swiss ADME web tool (http://swissadme.ch/index.php (accessed on 29 March 2023)) in terms of molecular properties, pharmacokinetics, drug-likeness, and medicinal chemistry.
The in silico toxicological evaluation for the most active compound 10a was performed using the web-service Cell-Line Cytotoxicity Predictor (CLC-Pred), which screens for in silico cytotoxicity on a panel of 278 tumor cells and 27 normal human cell lines from different tissues [91].

Conclusions
Five series of potential microtubule-targeting anticancer pyrrolo-fused heterocycles were designed as analogs of phenstatin and synthesized through 1,3-dipolar cycloaddition. Their antiproliferative activity was tested against NCI's panel of 60 cancer cell lines. Pyrrolo[1,2-a]quinoline 10a proved to be the best in terms of growth inhibition properties, with almost all GI 50 values in the submicromolar range. The best anticancer behavior was displayed on the A498 renal cell line (GI 50 = 27 nM), melanoma MDA-MB-435 cell line (GI 50 = 35 nM), and non-small cell lung cancer NCI-H522 line (GI 50 = 50 nM). In vitro, tubulin polymerization inhibitory properties were found for compound 10a, while PRISM analysis indicated a strong profile correlation to other microtubule inhibitors (even if with a totally different structure). QSAR simulation of ADMET properties showed a promising drug-likeness profile for compound 10a with very good selectivity for tumor cells, good physicochemical properties, good GI absorption, moderate solubility, and adequate bioavailability. The binding mode of 10a was investigated through molecular docking and molecular dynamics simulations, as well as configurational entropy calculations. This compound was confirmed to bind to the colchicine site of tubulin through global and local docking experiments, but molecular dynamics simulations revealed that some of the predicted interactions in the docking experiments were not stable. As such, the lowestscoring conformation from local docking, BM I, drifted away from its original orientation during MD simulations, while BM III, which was also identified through global docking experiments, was the most stable throughout the MD simulation. The second-lowest scoring conformation, BM II, which was the most similar in terms of relative orientation and hydrophobic amino acid interactions to other co-crystallized trimethoxy-substituted ring-containing tubulin binders including colchicine, was also stable throughout the MD simulation. However, longer MD simulations are necessary in order to demonstrate the full stability of the investigated systems. Configurational entropy calculations indicated similar values for configurational entropy loss upon binding for all three docked solutions, suggesting that all three binding modes are equally favorable in terms of entropic contribution. Taken together, our results suggest that for compound 10a, docking experiments alone are not sufficient for the adequate description of molecular interaction details in terms of target binding, which makes subsequent scaffold optimization more difficult to achieve. This highlights that a thorough in silico investigation should be performed for novel pyrrole-fused agents in order to have more meaningful insights into the molecular details of anticancer activity. Hopefully, our efforts will guide the discovery of novel potent antiproliferative compounds with pyrrolo-fused heterocyclic cores, especially from an in silico perspective.

Supplementary Materials:
The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/ph16060865/s1; Figure S1: HRMS spectrum of compound 10a; Figures S2-S19: NMR-spectra of active compounds; Figures S20-S45: NCI data for tested compounds; Figure S46: conformation distribution after blind docking of compounds (a) colchicine, (b) phenstatin, and (c) 10a; Figure S47: binding score distribution and RMSD from the lowest-scoring solution from docking experiments. Figure S48: local docking to α,β-tubulin and superimposition with colchicine for (a,d) BM I, (b,e) BM II, and (c,f) BM III; Figure S49: superimposition of BM II of 10a from local docking and D64131 (PDB ID 6K9V); Figure S50: superimposition of BM I of 10a from local docking and the last frame of the MD simulation; Figure S51: ligand root mean square fluctuation (RMSF) throughout the simulations for (A) 10a and (B) colchicine; Figure S52: timeline representation of the number of H-bond interactions throughout the MD simulations; Figure S53: timeline representation of the number of polar interactions throughout the MD; Figure S54. timeline representation of the number of all interactions throughout the MD simulations; Figure S55. configurational entropy for free 10a, BM I, BM II, and BM III when bound to tubulin, and the difference between bound and free state.