Modeling the Antileukemia Activity of Ellipticine-Related Compounds: QSAR and Molecular Docking Study

Abstract

The antileukemia cancer activity of organic compounds analogous to ellipticine representes a critical endpoint in the understanding of this dramatic disease. A molecular modeling simulation on a dataset of 23 compounds, all of which comply with Lipinski’s rules and have a structure analogous to ellipticine, was performed using the quantitative structure activity relationship (QSAR) technique, followed by a detailed docking study on three different proteins significantly involved in this disease (PDB IDs: SYK, PI3K and BTK). As a result, a model with only four descriptors (HOMO, softness, AC1RABAMBID, and TS1KFABMID) was found to be robust enough for prediction of the antileukemia activity of the compounds studied in this work, with an R² of 0.899 and Q² of 0.730. A favorable interaction between the compounds and their target proteins was found in all cases; in particular, compounds 9 and 22 showed high activity and binding free energy values of around −10 kcal/mol. Theses compounds were evaluated in detail based on their molecular structure, and some modifications are suggested herein to enhance their biological activity. In particular, compounds 22_1, 22_2, 9_1, and 9_2 are indicated as possible new, potent ellipticine derivatives to be synthesized and biologically tested.

1. Introduction

Cancer disease represents one of the most significant health problems in the world, being the second most common cause of death around the planet [1,2]. It is estimated that by 2030, the number of cancer cases in the world will have increased by approximately 23.6 million [2]. The situation in the Latin American population is even more dramatic; according to the World Health Organization (WHO) data, 37% of cancer cases are reported in this region, and are associated with the limited development of technology and industry related to the treatment of this disease in these countries [3]. Childhood cancer is considered a rare disease, with between 0.5% and 3% of the malign neoplasms reported around the world found in children. Nevertheless, this cancer type represents a significant public health problem due to the high probability of death in children, which generates a tremendous social impact on the families [4].

In this respect, leukemia is the principal childhood cancer. For children under 15, this type of cancer has a 30% rate of morbidity. Moreover, it seems that this kind of cancer affects more children living in poorer countries. In Latin America, the statistics is alarming: in Colombia, leukemia is responsible for around 49% of cancer deaths, while in Mexico, this kind of cancer is responsible for 22% of childhood deaths [5]. In contrast, in industrially and technologically developed countries, the picture that emerges is encouraging: survival in leukemia patients has increased notably in the last years. Around 90% survival rates are reported in countries such as the United States and the United Kingdom. According to these statistics, 10% of patients were non-responsive to primary chemotherapy and consequently died [4,5,6].

Nowadays, a major concern in public health is the quality of life of childhood leukemia patients, especially related the side effects and toxicity of the drugs used to treat this kind of disease. In this regard, studies of leukemia biology have increased in the past decade [7,8]. Some recent articles have reported the discovery of protein targets that play an essential function in the life cycle of leukemia cells. This knowledge represents a subtle but real thread of hope in the search for new and novel compounds with a high affinity for this target; additionally, understanding the action mechanisms of thse targets could help researchers to improve the activity and toxicity of cancer drugs [9,10,11,12].

In recent decades, quantitative structure–activity relationships (QSAR) are among the techniques most often used to generate new, promising compounds against diseases such as malaria [13], diabetes [14], and cancer [15], among others. By definition, a QSAR model is an equation that involves molecular descriptors with a remarkable influence on a particular biological activity. The knowledge of these chemical characteristics will allow the development of new compounds with better activity and better therapeutic indexes. Moreover, it has been shown that the QSAR-based synthesis has a high probability of success in any disease [16,17].

At the same time, another technique widely used in the development of new drugs is molecular docking [18,19]. This tool has become more and more critical in drug design because of the extraordinary advances in protein purification, nuclear magnetic resonance in proteins, and protein crystallography, which have contributed to elucidating the structural details in ligand–protein complexes. This technique can be used to model the interactions between small molecules and any protein at the atomic scale, allowing researchers to describe the behaviors of these small molecules at the binding site (active site).

The knowledge of both a leader structure from QSAR and its interaction type at the active site of the protein target from docking techniques represents a promising path in the search for and the development of new series of molecules active against a particular disease [20,21,22,23,24,25,26].

Regarding leukemia, the literature reveals extensive effort paid to the search for new molecules with potential applications in the treatment this disease. QSAR techniques have been applied to several families of compounds to find the structural characteristics that influence antileukemia activity [26,27,28]. More recently, some works have reported relationships between the structures of a compound family and affinity with a particular protein that plays a fundamental role in leukemia cell growth [29,30,31,32].

Ellipticine is a natural metabolite from the plant Ochrosia elliptica (Apocynaceae), and it has a potent anticancer activity. Its mechanism of action is related to DNA intercalation or via the inhibition of topoisomerase II protein [33,34,35]. Besides the significant activity of this compound, it is not yet available in the pharmaceutical market because of several side effects like nausea, vomiting, hypertension, and fatigue. We believe that strategies to eliminate or minimize the adverse effects of ellipticin can be derived based on functionalization, structural modification, or, more drastically, the search for new compounds analogous to this nucleus.

Among the compounds similar to ellipticine are the benzodioxinic analogues. Some studies have reported that the presence of oxygen atoms as a cyclic peroxide (dioxygen) or as an ester group on the ellipticine moiety seems to have a significant influence on the compound’s biological activity [36,37,38,39]. The structural similarity with ellipticine is, in theory, a fundamental key to designing new anticancer compounds. Moreover, the vast quantity of reported compounds with structures close to ellipticine represents a unique opportunity to identify the common structural characteristic that any molecular structure must have to be active against leukemia. In this respect, this work aimed to find a quantitative relationship between several molecular descriptors (topological, thermodynamics, and electronics) and the antileukemia activity of compounds related to ellipticine, in order to guide the synthesis of new promising antileukemia compounds. Additionally, to offer more insight into the interaction of ellipticine derivatives with leukemia cells, a docking calculation on the selected molecular target of the L1210 leukemia line cell is presented.

2. Results

Pharmacological data in vitro of several ellipticine analogues with antileukemia activity against L1210 cells were collected from the literature [40,41]. After the application of the Lipinski [42] rule filters, only 23 ellipticine analogues (ellipticine include) were selected. Figure 1 shows the chemical structure of the compounds studied herein.

2.1. Molecular Modeling

The minimum-energy 3D geometries for the compounds shown in Figure 1 were obtained using density functional theory with WB97XD/6-311G(d,p) as a theory level [43], using Gaussian 16 software [44] for Linux available in the high performing computer of the San Francisco de Quito University, Quito, Ecuador. The DFT level and interchange correlation functional was chosen because of its good correlation with experimental results based on the energetics and structure of organic molecules [45,46,47,48].

The minimum geometry structure was verified using the second derivative criteria [49]. In this regard, the vibrational frequency calculations were performed for the entire dataset and displayed no imaginary frequency, indicating that all of the geometries were minimum-geometry structures. Both minimum structures and frequency calculations were used to find electronic and molecular descriptors such as dipolar momentum (μ), HOMO (High Occupied Molecular Orbital) and LUMO (Low Unoccupied Molecular Orbital) energies, polarizability (α), enthalpy (H), entropy (S), free energy (G), ionization potential (PI), electronic affinity energy (EAE), hardness (η), softness (s), electrophilic index (ω), lipophilia (ClogP), polar surface area (PSA), topological index (TI), Balaban index (BI), hydrogen bond acceptor (HA), hydrogen bond donator (HD), AC[1]_K_F_AB_nCi_2_M1_NS0_C_LGL[8-9]_a_MID (AC1RABABMID), and TS[1]_K_F_AB_nCi_2_M1_SS0_T_LGL[2-3]_a_MID (TS1KFABMID). Except for AC1RABABMID and TS1KFABMID, these topological indexes were computed using Chemaxon [50], while AC1RABABMID and TS1KFABMID were calculated using QuBiLs-MIDAS as reported previously [51].

2.2. Statistical Analysis

The linear regression (LR) method was employed to find the most relevant molecular descriptors—the parameters related to the biological activity. Thus, the biological activity as pIC₅₀ (−ln(1/IC₅₀), dependent variable) was plotted against each molecular descriptor. A regression coefficient > 0.5 was considered to indicate an “important” molecular descriptor.

According to the linear regression analysis (LRA), Table 1 reveals that four electronic molecular descriptors (HOMO, softness, LUMO, and dipolar momentum (μ)) and three topological molecular descriptors (AC1RABABMID, TS1KFABMID, and PSA) have a strong influence over biological activity.

QSAR Model

Using the descriptors shown in Table 1 as independent variables, we obtained the mathematical models described by Equations (1)–(3), with three–four descriptors.

pIC₅₀ = 6.3074 + 1.2938 s + 19.794 AC1RABABMID + 1.5161 TS1KFABMIDR² = 0.798, F = 25.0, s = 1.023, Q² = 0.650, a(R²) = 0.064, a(Q²) = −0.453(1)

pIC₅₀ = 12.533 − 0.0189 PSA + 18.828 AC1RABABMID + 1.6302 TS1KFABMIDR² = 0.652, F = 11.9, s = 1.342, Q² = 0.506, a(R²) = 0.036, a(Q²) = −0.287(2)

pIC₅₀ = 20.530 + 51.942 HOMO + 20.763 AC1RABABMID + 1.5447 TS1KFABMID + 1.4174 sR² = 0.836, F = 23.0, s = 0.946, Q² = 0.729, a(R²) = 0.126, a(Q²) = −0.505(3)

The best model, selected according to the statistical robustness, was Model 3, with an R² of 0.836, a high Fisher ratio value of 27.02, and a great correlation prediction index (Q²). A Y-scrambling analysis was also performed on Model 3, and found values of a(R²) = 0.126 and a(Q²) = −0.505, which suggest that the models predictability cannot be explained by chance. The correlation prediction index was estimated by using the leave-one-out crossvalidation. Figure 2 graphically represents the linear relationship between the experimental pIC₅₀ values and those predicted using Equation (3).

The correlation matrix of Model 3 is shown in Table 2. Note that the Pearson’s correlation coefficient in each descriptor was <0.6, indicating that the model did not over fit. The descriptors AC1RABABMID and TS1KFABMID seemed to have the most influence over the activity.

According to Equation (3), the highest occupied molecular orbital, HOMO, has a large influence over the antileukemia activity. Thus, the more negative the HOMO energy, the greater the impact on biological activity. This descriptor is related to the molecule’s capability to participate in dipole/dipole interactions like hydrogen bonds. Additionally, the HOMO energy describes the ionization potential and the molecule’s vulnerability to electrophiles attack. HOMO energy also plays a pivotal role in free radical reactions and redox potential [52,53,54,55] therefore, it is very common to find this descriptor in works related to anticancer activity.

The global softness (s) related to the antileukemia activity shown for the compounds studied herein. Global softness is defined as the reciprocal of global hardness and describes the extent to which the electronic environment surrounding the nucleus/nuclei of an atomic/molecular species tends to loosen itself [56,57,58,59,60]. In this sense, it could be associated with a compound’s ability to deform its electronic cloud via variation of the number of electrons. Therefore, a soft compound necessarily has a high potential to transfer electrons in a redox process. As reported by Table 1, ellipticine had the highest softness value, while the rest of the studied compounds had values between 3 and 5, with the most actives ones values close to 5.

The topographic descriptors AC1RABABMID and TS1KFABMID, which were estimated using QuBiL-MIDAS software, have been widely used for the construction of QSAR models for several biological activities [51,61]. These descriptors presented a good correlation with the pIC₅₀ values, as shown in Table 2, and they were obtained by the application of some algebraic linear indexes, considering an autocorrelation and a total sum invariant with a lag value of 1 (AC1 and TS1) on the 3D optimized structure. These two descriptors were atomic weighted by the physicochemical property logP, denoted with “a” at the end of the descriptor, and this property is related to the water solubility as well as to the lipophilicity of the compounds, which both play an important role in almost any biological activity of an organic compound.

2.3. Molecular Docking

A molecular docking tool was used to gain more insight into the binding modes of the ellipticine-related compounds studied herein with the selected targets. It is necessary to mention that with the exception of ellipticine, the mechanisms of antileukemial cell proliferation inhibition by the compounds studied herein are not well known. One author proved that, in contrast to ellipticine, these kind of compounds are not able to inhibit topoisomerases I and II proteins [40,41]. In this regard, it seems that the antileukemial activity might be associated with other mechanisms that involve different targets than topoisomerases.

According to the literature, spleen tyrosine kinase (SYK), phosphoinositide 3’kinase (PI3K), and Bruton’s tyrosine kinase (BTK) were selected because they have been related to decrease of leukemia cells [62,63,64] in preclinical models. Morover, some compounds structurally close to those studied herein have been reported as SYK [65], BTK [66], and PI3K inhibitors [67]. In this regard, the entire dataset shown in Table 3 was docked against the three targets mentioned above. Due to the inactivity against topoisomerase II, this target was ruled out.

Fostamatinib (R788) is one of the latest reported inhibitors of both SYK and PI3K leukemia proteins; its action mechanism relates to the blocking of antigen-dependent B-cell-receptor signaling [66,67]; thus, based on its already known antileukemial activy, it was used as the control ligand in the docking study.

Table 3 reveals some interesting patterns: the most active compounds (pIC₅₀ > 13) showed scoring values smaller than −8 kcal/mol against any target. It is also of note that the less active compounds (pIC₅₀ < 11) had at least one value scored higher than −8 kcal/mol. A close inspection of Table 3 reveals that compounds 9, 20, and 22 can be highlighted because of their scoring values close to those of compound R788.

According to Table 1, the entire dataset could be divided into two structural groups; compound 13 derivatives and compound 22 derivatives. Compounds 9 and 22 had the highest antileukemia activity and the highest scoring values into their group. Thus, to gain more details about these two compounds, Figure 3 shows these molecules docked with the proteins SYK, PI3K, and BTK. However, because only the SYK protein has a reported inhibitor, this protein was selected to show the interactions with 9 and 22 compared with that inhibitor (fostamatinib). Thus, the molecular interaction between these ligands and SYK protein is shown in Figure 3; this was generated using Pymol [68] and ligplot software [69].

Figure 3 shows the most active compounds (9 and 22) and the structure of fostamitinib with the target protein, as well as a 2D diagram of interactions with the terminal residues of the protein. Further, the figure shows the pose of compounds 9 and 22 bound to the SYK protein; the inhibitor R788 (fostamatinib) is also displayed. It can be noted that compounds 9, 22, and fostamatinib contact with the protein via the same pocket. This result is fascinating because it suggests that these compounds might act via a similar mechanism. However, some differences were found in the molecular interaction types.

As shown in Figure 3, the molecular interactions between SYK and fostamatinib were mostly, dipolar in nature. The figure highlights the hydrogen bonds between Glu48 and Leu45 with the OH fragment, Lys40 and Gly18 with NH, and Met80 with the P–OH fragment in fostamatinib. In contrast, compound 9 did not present hydrogen bond formation with the target; however, other dipolar interactions, such as LEU133, Met82, and Glu81 with the CN aliphatic chain, and hydrophobic interactions such as cation–pi (Lys40) or pi–pi (PHE20), did take place.

Compound 22 showed a similar behavior to fostamatinib. This compound displayed several hydrogen bonds formed with the residues ASN131 and ARG130. Additionally, other dipolar interactions were noted, such as with PRO87, Asp144, and Val23. The rest of the interactions were in the hydrophobic range.

Compounds 9 and 22 were docked against phosphoinositide 3’kinase (PI3K). The 3D and 2D docking interaction diagrams are shown in Figure 4. As shown in Figure 4, both compounds 9 and 22 hit the PI3K protein at the same pocket as fostamitinib. However, the nature of the interaction was different. While fostamitinib interacts mostly via dipolar interactions (several hydrogen bonds), in general, compound 9’s activity was related to several hydrophobic interactions, mostly pi interactions; however, a few dipolar interactions could be observed (Hist650, Glu792, and Met788 with CO group). The same behavior was observed for compound 22. Figure 4 shows that just a dipolar interaction (Lys642 with OCH₃ group) took place for this compound.

According to the result above, it seems that the selected compounds might hit the proteins SYK and PI3K in the same pocket as the leader compound fostamitinib. This outcome represents a good starting point for the design and testing of these compounds and their derivatives experimentally.

As shown by the QSAR result, antileukemia compound activity demanded a low HOMO energy and high softness. Moreover, a close inspection of the molecular docking results revealed the relationship of scoring value and QSAR descriptors in Model 3: side chains with electronegative substitution drive to more polar molecules with smaller HOMO energy; thus, a better interaction with the molecular target could be.

In agreement with the above, i.e., the combination of both molecular behavior and QSAR results, we propose four molecules with the necessary chemical characteristic previously outlined. Figure 5 shows the predicted pIC₅₀ and scoring values for these compounds against SYK protein.

As shown in Figure 5, the QSAR model predicted both activity and scoring values higher than its precursor, encouraging the synthesis and experimental evaluation of the proposed compounds. However, as the scoring values for SYK were higher than for PI3K protein, only the docking against SYK was analyzed.

It was noetd that in the compound 9 derivatives, the enhancing of the side-chain with a hydrophobic substituent like a benzene ring increased the antileukemia activity. The benzene ring with an NH terminal not only improved the softness and diminished the HOMO energy, but also enhanced both hydrophobic and dipolar interactions with the selected target. This behavior is highlighted in the SYK proteins when interactions like hydrogen bonds (Arg130 and Asp126) take place.

On the other hand, of the compound 22 derivatives, 22_2 presented the best results against the SYK protein. It was noted that an aromatic ring bound to the ester moiety enhanced the biological activity. Once again, an aromatic ring meets the electronic requirements demands of the QSAR models, i.e., smaller HOMO energy and a higher softness. Additionally, the docking result from SYK target revealed some new interactions driven by further substitutions. For example, as in compound 22, the hydrogen bonds with ASn131 and Asp130 remain; however, a new hydrogen bond between the ester carbonyl group and Asp144 residue formed. Nevertheless, it seems that the hydrophobic interactions were more critical than for the precursor. Thus, several hydrophobic interactions appeared between the ester moiety aromatic ring and PHE20, Asn19, Val23, Asn19, Lys165, and Ser17 (Figure 6).

The combination of QSAR and docking results in this work seem to support Romero et al.’s idea [41] that the enlargement of the side-chain in ellipticine-related compounds enhance their antileukemia activity. When this enlargement of the side-chain is accompanied by an aromatic ring and electronegative atoms like N, the antileukemia activity is enhanced. Thus, the results reported herein support compounds 9 and 22. Likewise, their derivatives have a reasonable probability of being active against leukemia cells; therefore, their synthesis and posterior biological testing is encouraged.

The literature revealed that the synthesis procedure for compound 9_1 has been reported before [66]. However, it was not tested again cancer cells. The rest of the compounds have not been reported to date; however, some similar structure nuclei have been reported. In this regard, based on earlier reports, Figure 7 shows a proposing synthetic route for the four potential antileukemia compounds predicted in this work.

3. Materials and Methods

3.1. Pharmacological Data Collected

The data of in vitro ability to inhibit 50% L1210 leukemia cell proliferation for 23 ellipticine-related compounds were collected from the literature, mainly from Pujol et al. [40] and Romero et al. [41]. All data collection ensured homogeneity; all of the IC₅₀ data were obtained through the same experimental methods. The selected compounds’ IC₅₀ values < 100 μM.

3.2. Molecular Modeling

The minimum energy geometries for the 23 compounds studied herein were obtained using density functional theory with WB97XD/6-311G(2d,p) as the theory level, using Gaussian 16 software for Linux [44]. Both DFT theory level and interchange correlation functional were chosen because of their good correlation with experimental results based on the energetics and structure of organic molecules [45,46,47].

The minimum geometry structures were verified using the second derivative criteria [50]; the vibrational frequency calculations performed for the entire dataset showed no imaginary frequencies; therefore, all of the geometries were confirmed to be minimum-energy structures. Both the minimum structures and frequency calculations were used to find electronic and molecular descriptors such as dipolar momentum (μ), HOMO and LUMO energies, polarizability (α), enthalpy (H), entropy (S), free energy (G), ionization potential (PI), electronic affinity energy (EAE), hardness (η), softness (s), and electrophilic index (ω) using Equations (4)–(7).

(4)$\mu =-\chi =\frac{-\left(IP+EA\right)}{2}$

(5)$\eta =\frac{\left(IP-EA\right)}{2}$

(6)$S=\frac{1}{2\eta }$

(7)$\omega =\frac{{\mu }^2}{2\eta }$

where (µ) is electronic chemical potential and (X) is electronegativity.

The partition coefficient (ClogP), a measure of lipophilicity (a highly influenced descriptor), was obtained using MarvinSketch software for Windows 2017 [50], whereas Gibbs energy, enthalpy, and entropy were obtained by combining frequency calculation and statistical mechanics [49]. The 3D-optimized structures for the all data studied herein were used for the calculation of MAM descriptors, as described previously [14,15]. MAM descriptors are molecular attributes based on n-linear transformation matrix representations weighted over atomic properties like lipophilicity, van der Waals volume (ν), refractivity (r), polar surface area (PSA), polarizability (α), and electronegativity (χ), and these calculations were performed with QuBiLs-MIDAS software [51].

3.3. Statistical Analysis

The molecular descriptors were plotted against IC₅₀ (as pIC₅₀) to find the significance over the biological activity. Descriptors with higher association with antileukemia activity (r > 0.5) were considered statistically relevant and used in the construction of the mathematical model (QSAR).

3.3.1. Quantitative Structure–Activity Relationship (QSAR) and Statistical Validation

The relationships of antileukemia activity and the most relevant molecular descriptors were studied using multiple linear regressions. The method used has been described [13,70] before. Briefly, attempts were made to map the relationship between two or more independent variables with a dependent variable by fitting a linear equation involving the observed data. The independent variables of the model, selected according to forward selection and backward elimination [71,72,73] methods, were three statistical variables: the correlation coefficient (R), the Fisher ratio values (F), and the standard deviation (s).

After the QSAR models were formulated, they were validated statistically before their application in the designs of the new antileukemia molecules. To this end, the predictive power of the best equation was verified via leave-one-out cross-validation methods [74] and quantified by Equation (8). This method has been explained previously in several articles, and is well known for its extensive use in QSAR studies [75,76,77,78].

(8)$q_{cv}^2=1-\frac{{{\displaystyle \sum }}_{i=1}^n{\left(Y_{exp}-Y_{pred}\right)}^2}{{{\displaystyle \sum }}_{i=1}^n{\left(Y_{exp}-\overline{Y}\right)}^2}$

Likewise, the standard error of prediction (SEP) is calculated as:

(9)$SEP=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n{\left(yi-\ddot{y}i\right)}^2}{n}}$

where y is the experimental value of ln (1/IC₅₀), $\ddot{y}$ is the predicted value, and n is the number of samples used for model building. A y-scrambling analysis was also performed in order to confirm that the predictability of the models was not due to chance. The a(R²) and a(Q²) values were then determined using a total of 300 iterations on the response variable, and small values of these correlation coefficients were associated with a correlation not due to chance.

3.3.2. Molecular Docking

The protocol used herein to perform the molecular docking has been reported previously [68,72], and has been used to model the interactions of drugs with the active sites in different diseases, such as Pin1 (peptidyl-prolyl cis-trans isomerase NIMA-interacting 1) [15] inhibition and antimalarial activity [13]. Autodock4 software was used to this end [79]. Tridimensional structures of spleen tyrosine kinase (SYK, PDB-ID:4F4P), phosphoinositide 3´kinase (PI3K, PDB-ID:6DGT), and Buton´s tyrosine kinase (BTK, PDB-ID:3PIY) were obtained from the reported RX in the protein data bank (PDB) as.pdb format [80]. Each protein was refined before use; water molecules and any ligand associated with the protein were removed. Additionally, both polar hydrogen atoms and Kollman-type charge were added. The non-ligand proteins were then saved in .pdb format. Next, the optimized ligands from the WB97XD/6-311G(d,p) theory level were converted to .pdb format and added to the protein. Next, both protein and ligand saved in PDBQT format were used to perform the docking calculation with autodock4 software. The interaction points between ligands and proteins were analyzed using autogrid software, building the grid box into the active center of each protein considering 50 points through the x, y, and z directions, taking into account at least 2.5 million interactions through genetic algorithm, and taking into account a binding site size of 22 Angstrom. The best scoring factors obtained were compared with standard ligands, i.e., the reported inhibitors of the proteins used. These poses were saved in .pdb format and visualized using pymol for Linux. Next, the best pose interaction images were generated using Ligplot; this software allowed us to gain information about the molecular interactions (dipole–dipole and hydrophobic) around the protein active site and the ligand. The identified interactions were validated via the QSAR model.

4. Conclusions

The minima energy structures for 23 compounds—including ellipticine, a prominent anticancer compound—were obtained by means of density functional theory, combining the WB97XD method with a 6-311G(d,p) basis set. Multiple linear regression methods were applied over the entire dataset, using both electronic and topological molecular descriptors as independent variables and the pIC₅₀ as the dependent variable. After this stepwise approach was used, three mathematical models were generated; according to statistical descriptors, Model 3 had the best results. Moreover, this model was statistically validated through leave-one-out cross-validation methods. According to the results herein, antileukemia activity can be attributed to four molecular descriptors: HOMO, softness, AC1RABAMBID, and TS1KFABMID. These manageable descriptors allowed association of the chemical characteristics with antileukemia activity. The entire dataset was used to perform docking studies with three fundamental protein targets: spleen tyrosine kinase (SYK), phosphoinositide 3’kinase (PI3K), and Bruton’s tyrosine kinase (BTK). Compounds 9 and 22 were highlighted based on their scoring values; moreover, these compounds hit the SYK protein into the same active site as compound R788 (fostamatinib), a reported inhibitor of SYK protein. Furthermore, the scoring values of the four compounds were close to the R788 values. Based on QSAR analysis and the analysis of docking molecular interactions with SYK protein, four compounds are proposed as possible SYK inhibitors; the results achieved herein encourage the synthesis of these four compounds and subsequent biological tests.

Author Contributions

Conceptualization, J.R.M. and E.M.; Methodology, J.R.M., L.C., and E.M.; Software, J.R.M.; L.C. and E.M.; Validation, V.F.-M.; J.R.M. and E.M.; Formal analysis, V.F.-M., J.R.M. and E.M.; Investigation, V.F.-M., J.R.M., D.I., L.C. and E.M.; Writing—original draft preparation, J.R.M., E.M. and D.I.; Writing—review and editing, J.R.M. and E.M.; Project administration, J.R.M. and E.M.; Funding acquisition, V.F.-M., J.R.M. and E.M. All authors have read and agreed to the published version of the manuscript.

Funding

The authors are grateful with the Instituto of Biomedicina and SINDE-1831-2019 of Universidad Católica Santiago de Guayaquil (UCSG) for the financial support to this research. The High-Performance Computing systems available in the USFQ and UNINORTE were used for this project.

Conflicts of Interest

The authors declare no conflict of interest.

Figures and Tables

Figure 1. Chemical structures of ellipticine-analogous compounds studied in this work.

Figure 2. Experimental pIC50 vs. pIC50 predicted using Equation (3).

Figure 3. On the left, the 3D representation is shown, and on the right, the 2D interactions of compound 9 (a), fostamitinib (b), and compound 22 (c) at the active site of spleen tyrosine kinase (SYK).

Figure 3. On the left, the 3D representation is shown, and on the right, the 2D interactions of compound 9 (a), fostamitinib (b), and compound 22 (c) at the active site of spleen tyrosine kinase (SYK).

Figure 4. On the left, the 3D representation is shown, and on the right the 2D interactions of compounds 9 and 22 at the active site of PI3K.

Figure 5. Structures, pIC50, and scoring values against SYK protein for new proposed compounds based on QSAR results.

Figure 6. On the left, the 3D representation is shown, and on the right, the 2D interactions of the proposed compound at the spleen tyrosine kinase (SYK) active site.

Figure 7. Possible synthetic schema for proposed anticancer compounds. (1) Compounds 9_1 and 9_2. (2) Compounds 22_1 and 22_2.

Table 1

Values of the seven most important properties calculated and their respective values of antileukemial activity. HOMO, softness, and LUMO are expressed in hartree/particle, while μ is expressed in Debye.

Compound	pIC50	HOMO	AC1RABABMID	TS1KFABMID	Softness (s)	LUMO	μ	PSA (A)
1	12.94	−0.288	0.015	0.762	3.615	0.011	5.636	87.38
2	13.23	−0.288	0.007	0.187	4.613	0.011	5.646	91.98
3	13.12	−0.288	0.038	0.142	3.614	0.011	5.610	96.58
4	12.45	−0.291	0.011	0.197	3.576	0.011	5.575	101.18
5	11.21	−0.294	0.010	0.767	3.535	0.011	5.475	105.55
6	10.87	−0.288	0.015	0.342	3.502	0.003	5.099	102.75
7	9.47	−0.295	0.010	−0.029	3.427	0.003	9.575	73.94
8	14.43	−0.291	0.214	−0.199	4.528	0.007	6.154	73.94
9	9.33	−0.289	0.022	−0.303	3.612	0.012	6.567	88.38
10	11.27	−0.289	0.061	0.036	3.685	0.011	4.510	73.94
11	12.94	−0.281	0.022	0.169	4.011	0.015	6.952	69.18
ELLIPTICINE	15.48	−0.281	0.095	−1.020	6.100	0.000	3.545	24.06
12	13.63	−0.286	0.169	−1.020	3.495	13.630	−0.286	23.59
13	13.63	−0.286	0.169	−1.020	3.495	0.007	3.274	30.49
14	10.59	−0.274	0.007	−0.193	3.914	0.016	5.912	59.95
15	11.01	−0.262	0.055	−0.093	2.046	0.016	6.150	69.18
16	10.82	−0.270	0.014	−0.453	3.835	0.007	5.121	41.93
17	11.62	−0.287	0.027	0.382	3.369	−0.010	3.638	52.93
18	11.17	−0.285	0.046	−0.018	3.410	−0.009	2.842	41.93
19	14.51	−0.274	0.019	1.636	3.571	−0.008	7.458	100.85
20	15.42	−0.274	0.022	2.150	3.567	−0.008	7.597	100.85
21	10.11	−0.278	0.026	−0.660	3.560	−0.003	6.429	82.39
22	17.73	−0.288	0.201	2.337	3.292	−0.015	7.236	100.85

Table 2

Pearson’s correlation coefficient for the molecular descriptors used in Equation (3).

	pIC50	HOMO	AC1RABABMID	TS1KFABMID	Softness
pIC50	1	0.012	0.497	0.512	0.210
HOMO		1	−0.151	0.412	0.285
AC1RABABMID			1	−0.010	0.043
TS1KFABMID				1	−0.400
Softness					1

Table 3

Scoring values for the entire data studied against the three selected target proteins. SYK: spleen tyrosine kinase; PI3K: phosphoinositide 3’kinase; BTK: Bruton’s tyrosine kinase; R788: fostamatinib.

Compound	pIC50	SYK	PI3K	BTK
1	12.94	−8.1	−8.3	−9.0
2	13.23	−8.9	−9.5	−8.9
3	11.84	−8.3	−8.2	−8.9
4	13.12	−8.8	−9.1	−9.1
5	12.45	−7.6	−8.1	−9.0
6	11.21	−7.6	−8.6	−9.7
7	10.87	−7.3	−7.5	−7.9
8	13.00	−9.5	−9.2	−9.1
9	14.43	−9.8	−9.0	−9.8
10	9.32	−7.7	−8.4	−8.9
11	11.27	−8.8	−8.9	−8.9
12	12.94	−8.4	−8.4	−8.8
Ellipticine	15.48	−9.5	−9.9	−9.5
13	13.63	−9.1	−9.0	−8.4
14	10.58	−7.5	−7.4	−10.1
15	11.06	−9.0	−10.0	−10.5
16	10.83	−7.7	−6.9	−6.8
17	11.61	−7.7	−7.9	−7.4
18	11.16	−7.5	−9.8	−10.1
19	14.51	−8.9	−9.3	−10.1
20	15.42	−8.9	−8.7	−8.9
21	10.11	−7.7	−7.3	−6.7
22	17.73	−9.7	−8.7	−9.7
R788	13.23	−9.1	−9.5	−9.4