Molecular Modeling and In Vitro Evaluation of Piplartine Analogs against Oral Squamous Cell Carcinoma

Cancer is a principal cause of death in the world, and providing a better quality of life and reducing mortality through effective pharmacological treatment remains a challenge. Among malignant tumor types, squamous cell carcinoma-esophageal cancer (EC) is usually located in the mouth, with approximately 90% located mainly on the tongue and floor of the mouth. Piplartine is an alkamide found in certain species of the genus Piper and presents many pharmacological properties including antitumor activity. In the present study, the cytotoxic potential of a collection of piplartine analogs against human oral SCC9 carcinoma cells was evaluated. The analogs were prepared via Fischer esterification reactions, alkyl and aryl halide esterification, and a coupling reaction with PyBOP using the natural compound 3,4,5-trimethoxybenzoic acid as a starting material. The products were structurally characterized using 1H and 13C nuclear magnetic resonance, infrared spectroscopy, and high-resolution mass spectrometry for the unpublished compounds. The compound 4-methoxy-benzyl 3,4,5-trimethoxybenzoate (9) presented an IC50 of 46.21 µM, high selectively (SI > 16), and caused apoptosis in SCC9 cancer cells. The molecular modeling study suggested a multi-target mechanism of action for the antitumor activity of compound 9 with CRM1 as the main target receptor.

Molecules 2023, 28, 1675 2 of 25 usually located in the mouth, and it represents about 90% of cancers in the tongue and floor of the mouth. The pathology is considered a principal cause of morbidity and mortality worldwide since it is a highly complex and very aggressive cancer that presents a high probability of evolution towards both locoregional and distant metastasis [6][7][8][9].
In research involving medicinal plants, certain species of the Piperaceae family have become the subject of chemical-pharmacological studies due to their potential as a source of bioactive secondary metabolites, especially against tumor cells [10][11][12]. According to Costalotufo et al., the substance piplartine, found in some species of the genus Piper, has a trimethoxylated phenylpropanoid structure and can be classified as an alkamide. Piplartine is also known as piperlongumine [13].
The present study aimed to investigate the cytotoxic potential of a collection of piplartine analogs against SCC9 tumor cells derived from oral squamous cell carcinoma. For structure-activity relationships and evaluation of the results obtained in the biological tests, parameters such as lipophilicity, steric effects, and electronic effects were considered for each compound.

Chemistry
Preparation of the compounds is demonstrated in Scheme 1. The compounds were prepared by Fischer esterification (a); esterification using alkyl or aryl halides (b); and amidation reactions via the PyBOP coupling agent (c).
Molecules 2023, 28, x FOR PEER REVIEW 2 of 25 in women, breast, lung, cervix, colorectal, thyroid, and stomach cancers are more common [3][4][5]. Squamous cell carcinoma-esophageal cancer (EC) is a type of malignant tumor usually located in the mouth, and it represents about 90% of cancers in the tongue and floor of the mouth. The pathology is considered a principal cause of morbidity and mortality worldwide since it is a highly complex and very aggressive cancer that presents a high probability of evolution towards both locoregional and distant metastasis [6][7][8][9].
In research involving medicinal plants, certain species of the Piperaceae family have become the subject of chemical-pharmacological studies due to their potential as a source of bioactive secondary metabolites, especially against tumor cells [10][11][12]. According to Costalotufo et al., the substance piplartine, found in some species of the genus Piper, has a trimethoxylated phenylpropanoid structure and can be classified as an alkamide. Piplartine is also known as piperlongumine [13].
The present study aimed to investigate the cytotoxic potential of a collection of piplartine analogs against SCC9 tumor cells derived from oral squamous cell carcinoma. For structure-activity relationships and evaluation of the results obtained in the biological tests, parameters such as lipophilicity, steric effects, and electronic effects were considered for each compound.
The 1 H and 13 C NMR spectra of the products synthesized displayed comparable patterns whose signal multiplicity was in accordance with the proposed structures. For the hydrogen signals obtained, two aromatic hydrogens were observed as a singlet ranging from δ H 7.20 to 7.36 ppm and methyl shifts from methoxyls appeared as a singlet ranging from δ H 3.89 to 4.06 ppm. For the 13 C NMR spectral data, the common signals obtained

Cell Viability (Cytotoxicity) Assay
To determine the new compounds' cytotoxicity, the ester and amide collections derived from 3,4,5-trimethoxybenzoic acid were first tested on SCC9 cells and fibroblasts by the MTT cell viability assay. Initial screening was performed exclusively in SCC9 cancer cells, as they exhibit a lower stringent phenotype (higher sensitivity) to cytotoxic agents than other SCC cell lines, as previously demonstrated [14]. In addition, DMSO was used as a negative control, and Carboplatin, which is a chemotherapy drug widely used in the clinic for the treatment of this type of cancer, was used as a positive control.
Nonlinear regression curves of viability were generated in relation to the concentration of treatment used to calculate the IC 50 of each compound (Table 1). Only compounds with IC 50 lower than 500 µM were tested in fibroblasts to select compounds cytotoxic enough to be useful in future pre-clinical studies and to avoid crystallization at higher concentrations needed to calculate the IC 50 and selectivity in normal fibroblasts. Compounds that did not reach IC 50 below 1000 µM were considered non-determined (N.D.; Table 1) for the same reason.              Determination of IC 50  Different species of the genus Piper have the alkamide piplartine, which has several biological activities, such as antitumor, cytotoxic, antinociceptive, antiplatelet, and antimetastatic activity [15]. Although piplartine was not tested in this study, the compounds tested are analogous to this molecule. The most active compounds that fit the requisites in SCC9 cells were 3 (IC 50 = 299 ± 0.06), 4 (373.1 ± 0.05), 5 (204.2 ± 0.07), 6 (160.2 ± 0.07), 7 (257.6 ± 0.08), 9 (46.21 ± 0.16), and 13 (468.8 ± 0.07). Thus, normal fibroblasts were treated with the selected compounds and with Carboplatin to determine selectivity (Table 1). Of these compounds selected for testing on non-tumor cells (normal primary oral fibroblasts), only 7 was not selective (I.S. < 2); others showed intermediary selectivity, as in 5 and 6 (I.S. between 2 and 3), some were highly selective (I.S. > 4), and 9 was the most selective with an S.I. of 16.02. Compounds that present an S.I. ≥ 3 are considered highly selective, because they are three times more cytotoxic to the neoplastic lineage compared to the normal cell and therefore have a large therapeutic window [16]. It is worth noting that 9 was one of the novel compounds and had a 6.5-fold higher selectivity than Carboplatin (S.I. = 2.46), which is the drug of choice in the clinic.

Selectivity
To further characterize the cytotoxicity of the most selective compounds (S.I. > 4), we used another OSCC cell line, SCC4 (Table 2). Of these compounds, only 4 (IC 50   The IC 50 (µM) of two different OSCC cell lines (SCC9 and SCC4) and normal fibroblast cells followed by the average of the selectivity index of compounds 3, 4, 9, and 13 were calculated in the same way as in Table 1. All values are results of at least three independent experiments.

Hemolytic Activity
Since compound 9 had the highest cytotoxicity and selectivity against all cancer cells tested, we further analyzed its tolerability for future in vivo assays by discarding any surfactant activity of the compound which could lead to unspecific cytotoxicity through cellular membrane damage by hemolytic assays. Compound 9 displayed no significant hemolysis ( Figure 1A; lower than 1%) in human erythrocytes even in a very high concentration (1000 µM); hemolysis levels were comparable to the controls (carboplatin and DMSO) [17].

Hemolytic Activity
Since compound 9 had the highest cytotoxicity and selectivity against all cancer cells tested, we further analyzed its tolerability for future in vivo assays by discarding any surfactant activity of the compound which could lead to unspecific cytotoxicity through cellular membrane damage by hemolytic assays. Compound 9 displayed no significant hemolysis ( Figure 1A; lower than 1%) in human erythrocytes even in a very high concentration (1000 µ M); hemolysis levels were comparable to the controls (carboplatin and DMSO) [17].

Investigation of Cell Death Pathway
In order to elucidate the compound-9-induced cell death pathway, flow cytometry assays were performed to analyze apoptotic phenotypes (phosphatidylserine exposition, DNA fragmentation, and caspase 3/7 activation) and cell cycle distribution ( Figure 1B-E). Compound 9 showed a four-fold increase (~17% versus~4% in control) in annexin V singleand/or double-positive (PI + ) labeling at 48 h ( Figure 1B), which is indicative of greater apoptotic processes compared to the control. Compound 9 was also able to induce DNA fragmentation (~21%), as shown by increased SubG1 DNA content ( Figure 1C), and a massive executer caspases 3/7 activation (~66%; Figure 1D) when compared to the control (~4% SubG1 and~15% active caspases), supporting a classic apoptotic phenotype. Thus, our results suggest that 9 induces a classical apoptotic cell death pathway.
Furthermore, treatment of the SCC9 tumor cell with 9 led to a substantial increase in the number of cells exhibiting a G2/M DNA amount compared to the DMSO control ( Figure 1E;~59% ×~23%, respectively).
Together these results demonstrate that the piplartine analogue 9 is highly cytotoxic and selective against different types of cancers possibly through induction of cell cycle arrest at the G2/M phase and an intense cell death induction through apoptosis.
Driver genes can have negative, positive, or ambiguous regulation regarding tumorigenesis or cell proliferation, and individual analyses are required. PTPN1 inhibition has been ambiguous in tumorigenesis (acting as positive or negative regulator) [18][19][20]. A similar situation was found for PPARA [21]. Moreover, THRA (thyroid hormone receptor alpha) has been studied in thyroid, breast, and ovarian cancer [22,23]. However, the resulting protective effect from its direct inhibition is not as clear as for the beta type receptor [24]. On the other hand, ACHE and CASP9 have a protective effect against tumor progression [25,26]. Therefore, these proteins were not selected for further modelling. Table 4 summarizes the potential targets of compound 9 that were selected from the computational target fishing analyses and from the scientific literature that could explain its antitumoral activity. The table includes the UniProt accession of each target, the IDs employed along the manuscript, a brief functional description, and the source of the structure employed for modeling. The structures of the receptors were mainly obtained from the Protein Data Bank (PDB) database with the exceptions of STAT3 and AHR. In the case of AHR, no structure was available, and a homology model was generated. On the other hand, there were X-ray structures deposited in the PDB for STAT3, but these lack an important segment in the C-terminal region important for ligand binding. For this reason, the AlphaFold model was selected for modeling STAT3.   Table S1. In total, the 13 investigated targets lead to 33 ligand-receptor complexes to be further investigated. In all cases, these 33 complexes showed the compound within the binding sites. The weakest points of docking algorithms are the large approximations required in scoring functions and the limited exploration of the complex's conformational space. One of the alternatives to overcome these limitations is the post-processing of docking solutions with more accurate and computationally intensive methods based on Molecular Dynamics (MD) simulations. MD-based methods have aided in the identification of the molecular targets of bioactive compounds [27,28], while the estimation of free energies of binding from MD trajectories are effective at discriminating between correct and incorrect ligand-binding poses [29]. According to this evidence, the final selection of the most probable targets of compound 9 was carried out using the MD simulations protocol presented in the Materials and Methods section.
The length of the MD simulations required for predicting the free energy of binding of a ligand to a protein is a topic highly discussed in the scientific literature. Usually, the length of MD simulations is directly related to the available hardware resources. Additionally, there is no consensus on the optimal length of MD simulations required to obtain a conformational ensemble for MM-PBSA calculations. Different authors have concluded that short (fewer than 5 ns) simulations would be sufficient for MM-PBSA calculations [30,31]. Another aspect considered in the literature is that when short MD simulations are used for MM-PBSA calculations, a multiple trajectories approach is preferred over a single trajectory one. Thus, we performed five short MD simulations per complex, each one lasting for 4 ns, and the free energies of binding were estimated from an ensemble of 100 MD snapshots evenly extracted from all five trajectories. With this approach, the total MD simulation time performed in the current research was 660 ns.
The conformational stability of the complexes during MD simulations was measured as the RMSD of the ligand relative to the starting docking pose along each of the 5 ns MD simulations. These plots are provided as Supplementary Materials in Figures S1-S33. Overall, in most of the simulations, the ligand deviated less than or around 2 Å from the starting docking pose during the 4 ns MD simulation time. In no case did the observed RMSD values increase beyond 3 Å. This indicates that the obtained MD trajectories were conformationally stable. Furthermore, the differences in the RMSD plots obtained for a target suggest that a complex explores different conformational spaces during each MD replica.
The predicted free energies of binding obtained from the MD simulations are given as Supplementary Materials in Table S2 and summarized in Figure 2. Only the top-scored ligand conformer per target is represented in the figure. According to these results, it can be concluded that the most probable target of the compound is CRM1, followed by AHR, MCL1, and RIPK2. This suggests a multi-target mechanism of action for the antitumoral activity of the compound, with CRM1 as the principal receptor for compound 9. Consequently, the predicted binding mode of the ligand to CRM1 was analyzed in detail.
snapshots evenly extracted from all five trajectories. With this approach, the total M simulation time performed in the current research was 660 ns.
The conformational stability of the complexes during MD simulations was measur as the RMSD of the ligand relative to the starting docking pose along each of the 5 ns M simulations. These plots are provided as Supplementary Materials in Figures S1 to S Overall, in most of the simulations, the ligand deviated less than or around 2 Å from t starting docking pose during the 4 ns MD simulation time. In no case did the observ RMSD values increase beyond 3 Å . This indicates that the obtained MD trajectories w conformationally stable. Furthermore, the differences in the RMSD plots obtained fo target suggest that a complex explores different conformational spaces during each M replica.
The predicted free energies of binding obtained from the MD simulations are giv as Supplementary Materials in Table S2 and summarized in Figure 2. Only the top-scor ligand conformer per target is represented in the figure. According to these results, it c be concluded that the most probable target of the compound is CRM1, followed by AH MCL1, and RIPK2. This suggests a multi-target mechanism of action for the antitumo activity of the compound, with CRM1 as the principal receptor for compound 9. Con quently, the predicted binding mode of the ligand to CRM1 was analyzed in detail.  Figure 3 represents the predicted pose of compound 9 in the binding cavity of CRM as well as the ligand-receptor interactions. The complex structure selected for represen tion was the centroid of the most populated cluster obtained from clustering the MD sna shots used for MM-PBSA calculations. The ligand was predicted to bind in the know mainly hydrophobic groove of CRM1 responsible for the recognition of the nuclear exp signals in cargo proteins [32]. Among the 21 HEAT repeats that form CRM1, the inhibito binding cleft was at the interface of repeats 11 and 12.
Specifically, the methoxy of the para-methoxyphenyl group was predicted to ori tate toward C528. This is a key residue that is involved in the covalent binding of ma known CRM1 inhibitors [32,33]; however, non-covalent inhibitors of this protein have a been reported [34]. The rest of this ring is predicted to interact with L525, A541, I544, K5 and F572. The central linker of the compound interacts with L525, T564, and V565. On other hand, the trimethoxyphenyl moiety is predicted to occupy a cavity lined by L5 V518, I521, V548, F554, L555, and F561, positioning favorably to form π-π stacking int actions with F561.  Figure 3 represents the predicted pose of compound 9 in the binding cavity of CRM1 as well as the ligand-receptor interactions. The complex structure selected for representation was the centroid of the most populated cluster obtained from clustering the MD snapshots used for MM-PBSA calculations. The ligand was predicted to bind in the known mainly hydrophobic groove of CRM1 responsible for the recognition of the nuclear export signals in cargo proteins [32]. Among the 21 HEAT repeats that form CRM1, the inhibitors' binding cleft was at the interface of repeats 11 and 12.
Specifically, the methoxy of the para-methoxyphenyl group was predicted to orientate toward C528. This is a key residue that is involved in the covalent binding of many known CRM1 inhibitors [32,33]; however, non-covalent inhibitors of this protein have also been reported [34]. The rest of this ring is predicted to interact with L525, A541, I544, K568, and F572. The central linker of the compound interacts with L525, T564, and V565. On the other hand, the trimethoxyphenyl moiety is predicted to occupy a cavity lined by L517, V518, I521, V548, F554, L555, and F561, positioning favorably to form π-π stacking interactions with F561. Figure 3. Predicted binding mode and ligand-receptor interactions for the compound 9-CRM1 complex. The ligand is represented as orange balls and sticks. The receptor and its surface are colored as carbon: tan, nitrogen: blue, oxygen: red, and sulfur: yellow. The figure was produced with UCSF Chimera [35] and LigPlot+ [36].
Our modeling results were consistent with the observed antitumoral activity of compound 9. CRM1 is a validated molecular target for anticancer compounds [37][38][39]. This protein is usually overexpressed in cancer and is classified as an oncogenic driver. Specifically, piperlongumine has been shown to be an inhibitor of CRM1 with antitumoral properties [40]. The predicted interaction of compound 9 with C528 is in agreement with the previously observed direct interaction of piperlongumine with this amino acid in CRM1 [40]. The models described here will aid in directing future research focusing on improvement of the antitumoral and CRM1 inhibitory properties of piperlongumine derivatives.
Finally, the ADMET properties of compound 9 were predicted as described in the Materials and Methods section. The results of these predictions are given in Table 5, and they show that compound 9 contains no PAINS alert. In addition, it has favorable physicochemical properties for oral bioavailability. Compound 9 was also predicted to be soluble in water with high gastrointestinal absorption, to be skin permeable, and to not be a substrate of the P-gp. It was predicted to poorly distribute to the brain and to be an inhibitor of three isoforms of the Cytochrome P450. This last aspect, along with its predicted mutagenic potential, should be fully clarified with additional experiments during the hit optimization phase of this compound.  Our modeling results were consistent with the observed antitumoral activity of compound 9. CRM1 is a validated molecular target for anticancer compounds [37][38][39]. This protein is usually overexpressed in cancer and is classified as an oncogenic driver. Specifically, piperlongumine has been shown to be an inhibitor of CRM1 with antitumoral properties [40]. The predicted interaction of compound 9 with C528 is in agreement with the previously observed direct interaction of piperlongumine with this amino acid in CRM1 [40]. The models described here will aid in directing future research focusing on improvement of the antitumoral and CRM1 inhibitory properties of piperlongumine derivatives.
Finally, the ADMET properties of compound 9 were predicted as described in the Materials and Methods section. The results of these predictions are given in Table 5, and they show that compound 9 contains no PAINS alert. In addition, it has favorable physicochemical properties for oral bioavailability. Compound 9 was also predicted to be soluble in water with high gastrointestinal absorption, to be skin permeable, and to not be a substrate of the P-gp. It was predicted to poorly distribute to the brain and to be an inhibitor of three isoforms of the Cytochrome P450. This last aspect, along with its predicted mutagenic potential, should be fully clarified with additional experiments during the hit optimization phase of this compound.

Discussion
Two ester derivatives (compounds 2 and 10) presented no cytotoxic activity against SCC9 cells. The others presented low to satisfactory cytotoxicity. In general, the esters were more promising than the amides. However, the starting material, 3,4,5-trimethoxybenzoic acid was inactive.
The alkyl derivatives with medium carbon chains were more potent than the aryl derivatives, except for 4-methoxy-benzyl trimethoxybenzoate (9), which was the most potent analog. In fact, compound 9 was 4.5 times more potent and 6.5 times more selective than Carboplatin, the positive control. Carboplatin is used clinically against various types of cancer, such as uterine-cervix cancer, lung carcinoma, gynecological cancers, breast cancer, and others [41][42][43].
It was evident that reducing the size of the carbonic chain reduced the cytotoxic activity for alkyl derivative compounds 4, 5, and 6. However, when comparing straight-chain alkyl derivative propyl 3,4,5-trimethoxybenzoate (3) with isopropyl 3,4,5-trimethoxybenzoate (4), which is a branched-chain alkyl derivative, greater potency was noted for compound (3). N-butyl-3,4,5-trimethoxybenzoamide (13) was the only bioactive derivative, yet comparison with butyl 3,4,5-trimethoxybenzoate (5) revealed that the ester was 2.3 times more potent. Derivative compound 5 was equipotent to the reference drug. Isopentyl 3,4,5-trimethoxybenzoate (compound 6) was slightly more potent than compound 5, and 2-methylnaphthalene 3,4,5-trimethoxybenzoate (11) presented greater potency and selectivity than diphenyl 3,4,5-trimethoxybenzoate (12), revealing that the fused (naphthalene) ring system potentiates cytotoxic activity. The results obtained in the present study corroborate the results obtained by Chen and collaborators (2018) [46], who evaluated piplartine against SAS and CGHNC8 cells, which are oral cancer cell lines. It was shown that piplartine suppressed tumor sphere formation in the oral cancer cells tested, suppressed expression of the SOX2, Oct-4, and NANOG genes associated with stem cells, inhibited cell migration and invasion, and regulated expression of molecules associated with epithelial-mesenchymal transition (EMT). It also increased chemo-and radiosensitivity while even inhibiting tumorigenesis in vitro and in vivo.
Compared to other chemotherapies and to new putative compounds for treatment of OSCC cells, the selectivity results for the tested derivatives were significantly high [47][48][49][50].
Comparing the results obtained in the present study with data from the literature reveals similar cytotoxicities (IC 50 ) for piplartine and other analogs against various tumor cells lines, but none demonstrated a higher selectivity than the results displayed by compound 9 [51][52][53].
Test derivative (9) results (as well as those of piplartine) did not reveal hemolytic activity in mouse erythrocytes, suggesting that the cytotoxicity of this compound class does not involve membrane damage [17].
Based on the cell death pathway trials, compound 9 is likely involved in apoptotic processes. The derivative was able to induce DNA fragmentation. A similar apoptotic phenotype was observed in PC-3 cells treated with piplartine, with accumulation of fragmented DNA (SubG1) and activation of caspase 3 [54]. Cell cycle arrest in the G2/M phase causes a delay in cell proliferation [55], and this arrest in G2/M has been observed in PC-3 and V79 cells treated with piplartine [54,56].

Materials and Methods
All reagents necessary for chemical synthesis were obtained from Sigma-Aldrich. Adsorption column chromatography (CC) was performed using silica gel 60, ART 7734 from MERCK as the stationary phase (particle sizes between 0.063 and 0.200 mm). The nuclear magnetic resonance ( 1 HNMR and 13 CNMR) spectra were recorded in Bruker Avance III HD spectrometers operating at 400 MHz ( 1 H) and 100 MHz ( 13 C) and with a Varian NMR operating at 500 MHz ( 1 H) and 125 MHz ( 13 C). Chemical shifts were expressed in parts per million (ppm), and the multiplicities of the 1 HNMR bands were indicated according to the conventions: s (singlet), d (doublet), t (triplet), quint (quintet), sext (sext), sept (septet), and m (multiplet). The Fourier Transform Infrared Spectroscopy (FTIR) spectra were obtained using an Agilent technologies Cary 630 FTIR instrument, in the 4000-400 cm −1 spectral range. High-resolution mass spectra were obtained at the Biomolecules Mass Spectrometry Center (CEMBIO) of the Health Sciences Center (CCS/UFRJ). The spectrometer used was Solarix XR-FT-ICR configuration (Bruker), with an external temperature with 100 µM sodium formate in water/acetonitrile 1:1 and electrospray ionization source (ESI-nebulizer pressure: 1 bar, 4, 0 L/min, temperature 200 • C, flow rate 4 µL/min), pos-itive mode, detection window 75 to 1200 m/z. The theoretical resolution was 99,000 for m/z 400. The samples were solubilized in 1 mL of methanol and then diluted 1000× in water:acetonitrile 1:1 and injected directly into the mass spectrometer. Melting points were determined in a Microquímica apparatus (Microquímica equipamentos LTDA, Model MQAPF 302, Serial No.: 403/18, Palhoça, Brazil) with temperature measurements in the 10 • C to 350 • C range. All reaction processes were monitored by analytical thin layer chromatography using silica gel sheet chromatographs (Merck 60 F 254).

Preparation of Compounds 1-6
Sulfuric acid (0.2 mL) was added slowly to a solution containing 3,4,5-trimethoxybenzoic acid (0.1 g; 0.47 mmol) solubilized in 20 mL of alcohol. The reaction mixture was refluxed with magnetic stirring for 5-24 h. After alcohol rotary evaporation, distilled water (10 mL) was added and extraction was performed with dichloromethane (3 × 10 mL). The organic phase was neutralized with 5% sodium bicarbonate (3 × 10 mL), treated with 10 mL of water, and dried with anhydrous sodium sulfate. After purification using column chromatography with silica gel 60 and a mixture of hexane and ethyl acetate in an increasing polarity gradient, yields ranging from 45.3 to 99.6% were obtained.

Preparation of Compounds 7-12
In a solution containing 3,4,5-trimethoxybenzoic acid 1 Eq. solubilized in 15 mL of anhydrous acetone, 4 Eq. of triethylamine, and 1 Eq. of the halide were added. The reaction mixture was refluxed with magnetic stirring for 72 h. After acetone rotary evaporation, distilled water (10 mL) was added and extraction was performed with dichloromethane (3 × 10 mL). The organic phase was treated with 10 mL of distilled water and dried with anhydrous sodium sulfate. After purification, using column chromatography with silica gel 60 and a mixture of hexane and ethyl acetate in an increasing gradient of polarity, yields ranging from 40.6 to 75.2% were obtained.

Preparation of Compounds 13-19
A solution containing 3,4,5-trimethoxybenzoic acid (0.1 g; 0.47 mmol) and 0.47 mmol of the corresponding amine solubilized in dimethylformamide (1.0 mL; 0.94 mmol) in the presence of triethylamine (0.06 mL; 0.47 mmol) at 0 • C was prepared. At 30 min, PyBOP coupling agent 0.24 g (0.47 mmol) was added, which had previously been dissolved in dichloromethane (1.0 mL; 0.94 mmol) with magnetic stirring. The reaction time was 6-24 h. After alcohol rotary evaporation, distilled water (10 mL) was added and extraction was performed with dichloromethane (3 × 10 mL). The organic phase was treated with 10% hydrochloric acid (3 × 10 mL) and then neutralized with 5% sodium bicarbonate (3 × 10 mL), treated with 10 mL of distilled water, and dried with anhydrous sodium sulfate. Yields ranging from 44.4 to 91.3% were obtained and purified using column chromatography with silica gel 60 and a mixture of hexane and ethyl acetate in the proportions of 6:4 or 3:7.  The viability of cancer cells and primary fibroblasts was tested using the MTT assay as described in MACHADO et al., 2021 [49]. Briefly, cancer cells (5 × 10 3 cells/well) were plated in duplicate in a 96-well plate, as indicated. DMSO was diluted in culture medium at the same concentrations and was used as a negative control, representing 100% cell viability. After 48 h of treatment, cells were incubated for 3.5 h with 5 mg/mL of MTT reagent (3,4,5-dimethiazol-2,5-diphenyltetrazoliumbromide) (Sigma-Aldrich Co., St. Louis, MO, USA) diluted in culture medium. Subsequently, the formazan crystals were dissolved by MTT solvent (50% Methanol and 50% DMSO), and the absorbance was read at 560 nm using the EPOCH microplate spectrophotometer (BioTekInstruments, Winooski, VT, USA).

Hemolysis Assay
Human blood was used as approved by the Research Ethics Committee of the Fluminense Federal University-Nova Friburgo, RJ (CAAE: 43134721.4.0000.5626). The erythrocytes were collected by centrifugation at 1500 rpm for 15 min, washed with PBS (phosphate buffer saline) and 10 mM glucose, and counted in an automatic cells counter (Thermo Fisher, Waltham, MA, USA). The erythrocytes were then plated in 96-well plates at a concentration of 4 × 10 8 /well in triplicate, and 10 µL of the compounds were added at a final concentration of 1000 µM in PBS with glucose (final volume of 100 µL). Ten microliters of PBS were used as negative control and 10 µL of PBS with 0.1% Triton ×100 as a positive control. Data reading was performed with EPOCH (BioTek Instruments, Winooski, VT, USA) at 540 nm absorbance, and statistical data were generated with the GRAPHPAD Prism 5.0 program (Intuitive Software for Science, San Diego, CA, USA).

Cell Cycle and SubG1 Analysis
To evaluate the action of the compound 9 on the cell cycle, cells of the SCC9 lineage were plated in a six-well plate (5 × 10 5 cells/well) as described in LUCENA et al., 2016. After 48h of treatment, cells were trypsinized and stained with propidium iodide (75 µM;Sigma) in the presence of NP-40 (Sigma). DNA content was analyzed by collecting 10,000 events using a FACScalibur flow cytometer. Data were analyzed using CellQuest (BD Biosciences) and FlowJo (Tree Star Inc., Ashland, OR, USA) software. (Version 9.9.4).

Caspase Analysis
For active caspase detection, 5 × 10 4 SCC9 cells were plated in a 24-well plate containing 1 mL of DMEM/F12 with 10% FBS per well. CellEvent™ Caspase-3/7 reagent (#R37111, Invitrogen) was diluted in culture medium according to the manufacturer's instructions. Twenty-four hours after plating, cells were treated with Caspase-3/7 reagent and 2× IC50 of 9 or DMSO as a control. Cells were analyzed by flow cytometry after 48 h of treatment.
4.4.6. Statistical Analysis, Calculation of IC 50 , and Selectivity Index (SI) IC 50 values for the MTT assay were obtained by non-linear regression using Graph-Pad Prism 5.0 software (GraphPad, San Diego, CA, USA) from at least three independent experiments. Data are presented as means, ± standard deviation (SD). A log dose response (inhibitor) vs. response curve using the least square method was used to determine the IC 50 and SD of the data. The selectivity index was calculated as: S.I. = IC 50 of the compound in fibroblasts/IC 50 of the same compound in tumor cells.
The computational methodology used for target prediction consisted in the consensus of different models of target-compound interaction predictions previously applied in other investigations [66,67]. Briefly, all the possible target proteins of compound 9 in Homo sapiens were predicted using the following algorithms: MolTarPred [68], Swisstarget Prediction [69], Targetnet Scbdd [70], Targetnet Scbdd Ensemble [70], RF QSAR [71], and PPB2 [72]. The following algorithms from PPB2 were employed for predictions: PPB2-Extended Connectivity fingerprint ECfp4 NN, PPB2-Shape and Pharmacophore fingerprint Xfp NN, PPB2-Molecular Quantum Numbers MQN NN, PPB2-Extended Connectivity fingerprint ECfp4 NNNB, PPB2-Shape and Pharmacophore fingerprint Xfp NNNB, PPB2-Molecular Quantum Numbers MQN NNNB, PPB2-Extended Connectivity fingerprint ECfp4 NB, and PPB2-Extended Connectivity fingerprint DNN. Therefore, a total of 13 algorithms were used to predict the possible target proteins of compound 9. For each target, a final score was computed as: where: F i is the average normalized score of target i across all methods that predict it, N is the number of methods identifying the i-th target, and M is the number of target fishing algorithms [73] employed. Usually, many proteins comprising different metabolic pathways or biological processes are found with this consensus strategy. However, only those involved in cancer are relevant for studying the potential antitumoral mechanism of action of compound 9. To detect possible proteins related to cancer we used two different databases: (1) DriverDBv3 [74] and (2) Cancer Target Discovery and Development (CTD 2 ) [75]. The first database is an in-tegrative multi-omics database comprising driver genes organized across different types of cancer. All the 1457 driver genes reported in the database were considered. It is important to notice that a driver gene does not necessarily imply that targeting it will produce an antitumoral activity. The second database is a team effort organized by the National Cancer Institute and the NIH with the aim of validating discoveries from large-scale adult and pediatric cancer-genome research and their use in precision medicine. From this database, 89 genes classified as molecular targets with different evidence levels (Tier 1 to 3) were selected.
It is possible that some of the potential targets identified by the consensus analysis are not necessarily targets for antitumoral compounds (as defined in the CTD 2 database) or a driver gene (as defined in the DriverDBv3 database). Therefore, we constructed a proteinprotein interaction network with the top 89 targets from CTD 2 restricted to only direct interactions using the BioGrid [76] database and Cytoscape [77]. The top 100 candidate targets obtained from the consensus target fishing approach were filtered to keep only those also present in the list provided by the analysis of the DriverDBv3 and CTD 2 databases.

Molecular Docking
Molecular docking was performed following the previously described methodology [78,79]. One initial three-dimensional (3D) structure of compound 9 was generated with OpenEye´s Omega [80,81], and partial atomic charges of type am1bcc were added to it with Molcharge [82]. Docking calculations were carried out with Gold [83].
The explored binding sites for compound 9 were defined from the ligands and inhibitors co-crystallized with the target proteins when these were present in their PDB files. This was the case for NFKB (PDB 1SVC), CRM1 (PDB 6TVO), PI3K (PDB 4JPS), mTOR (PDB 4JSP), GSTP1 (PDB 5J41), MTNR1B (PDB 6ME6), MMP2 (PDB 1HOV), RIPK2 (PDB 5J79), DUSP3 (PDB 3F81), and MCL1 (PDB 6UDV). For STAT3, despite using the model available in the AlphaFold repository, the ligand-binding cavity was defined after superimposing this model with an experimental structure of the protein in complex with an inhibitor (PDB 6NJS). In the case of AHR, the antagonist present in one of the model templates identified by the SwissModel server (PDB 4ZQD) was used as a reference for defining the binding cavity after structural superimposition with the homology model. On the other hand, for RELA the ligand-binding region was defined from the segment of DNA interacting with it on the X-ray structure (PDB 3GUT). In all cases, the detect cavity option of Gold was switched on, and the binding region was defined as any receptor residue at a distance lower than 6 Å from the reference ligand.
The ChemPLP scoring function was selected for primary docking, and 30 diverse binding hypotheses were generated for each target. The search efficiency parameter of Gold was set to 200%. The binding modes predicted were rescored with the scoring functions GoldScore, ChemScore, and ASP. Scoring values were converted to Z-scores and aggregated as in our previous publications [78,79]. All ligand conformations with a Z-score value larger than one were selected for further analyses.

Molecular Dynamics Simulations and Free Energies of Binding
Molecular dynamics (MD) simulations were performed with Amber 22 [84] as in our previous reports [79,85]. The ff19SB, gaff2, and lipid17 force fields were used for the parametrization of amino acids, compound 9, and lipids, respectively. The Zn 2+ ion and its coordinating residues in MMP2 were prepared according to the Cationic Dummy Atom (CaDA) model [86,87].
Except for the MTNR1B membrane protein, all systems were prepared following the same procedure. All these complexes were enclosed in truncated octahedron boxes that were solvated with OPC water molecules. Excess charges were neutralized by the addition of Na + and Cl − counterions at a concentration of 150 mM following the methodology described in [86]. The solvated systems were energy minimized in two stages with all atoms except those solvent restrained during the first of these. All constraints were released for the second energy minimization step. Energy minimized systems were heated from 0 K to 300 K and subsequently equilibrated at constant temperature and constant pressure.
The complexes with MTNR1B were prepared with the CHARM-GUI server [28,29]. Bilayer membranes containing 50 1-Palmotoyl-2-oleoylglycero-3-phosphocholine (POPC), 50 1-Palmitoyl-2-oleoylphosphatidylethanolamine (POPE), and 25 Cholesterol (CHL) molecules on each side were constructed to embed these systems. Excess charges were neutralized as in the soluble proteins (see above). The energy minimization, heating, and equilibration steps for these complexes were performed using the configuration files provided by CHARM-GUI.
Five short (4 ns) MD simulations were performed for each, including the membrane protein. The atomic initial velocities were randomly initialized for every production run to ensure a diverse exploration of the complexes' conformational space. Parameters for the production runs were the same for all the systems.
The free energies of binding were estimated using the MM-PBSA methods as implemented in the MMPBSA.py script of Amber 22. These calculations took place over 100 MD snapshots evenly selected from the 1 ns-4 ns time interval each of the five production runs. The ionic strength for MM-PBSA calculations was set to 150 mM. Default implicit solvent parameters were employed for the complexes of compound 9 with soluble proteins. For MTNR1B, a heterogeneous dielectric implicit membrane model (memopt = 3) was employed. The thickness and center of the implicit membrane was defined from average positions of the N31 atoms in the explicit bilayer used for MD simulations.

ADMET Predictions
The absorption, distribution, metabolism, excretion, and toxicology (ADMET) properties of compound 9 were predicted with the SwissADME [88] and pkCSM [89] servers. The prediction of the physicochemical properties as well as the PAINS analysis were performed with SwissADME, while the rest of the predictions were made with the pkCSM server.

Conclusions
This study presents preparation, cytotoxic evaluation, and molecular modeling of a series of piplartine analogs against tumor cells. SAR revealed the influence of carbon chain size on pharmacological activity. Compound 4-methoxy-benzyl 3,4,5-trimethoxybenzoate (9) presented the best cytotoxic profile (IC 50 = 46.21 µM with a selectivity index of 16.02). Its potency and selectivity were superior to the reference drug (carboplatin). The in silico study suggested that the most likely compound 9 target was CRM1, followed by AHR, MCL1, and RIPK2, suggesting a potential multi-target mechanism of action for its antitumor activity and corroborating the results obtained in the cytotoxicity assays. In accordance with the results obtained in the present study, compound 9 may also be useful in other studies to synthesize analogs in an attempt to optimize its activity. Other compounds derived from 3,4,5-trimethoxybenzoic acid may also be tested against differing cancer cell lines to assess their cytotoxicity. The results of our study may contribute to the development of new antitumor agents.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/molecules28041675/s1, Table S1: Results of docking compound 9 to its potential targets; Table S2: Predicted free energies of binding according to the MM-PBSA method; Figures S1-S33: RMSD of compound 9 relative to its docking pose along the MD simulations.