Computational Analysis of Artimisinin Derivatives on the Antitumor Activities

The study on antitumor activities of artemisinin and its derivatives has been closely focused on in recent years. Herein, 2D and 3D QSAR analysis was performed on the basis of a series of artemisinin derivatives with known bioactivities against the non-small-cell lung adenocarcinoma A549 cells. Four QSAR models were successfully established by CoMSIA, CoMFA, topomer CoMFA and HQSAR approaches with respective characteristic values q2 = 0.567, R2 = 0.968, ONC = 5; q2 = 0.547, R2 = 0.980, ONC = 7; q2 = 0.559, R2 = 0.921, ONC = 7 and q2 = 0.527, R2 = 0.921, ONC = 6. The predictive ability of CoMSIA with r2 = 0.991 is the best one compared with the other three approaches, such as CoMFA (r2 = 0.787), topomer CoMFA (r2 = 0.819) and HQSAR (r2 = 0.743). The final QSAR models can provide guidance in structural modification of artemisinin derivatives to improve their anticancer activities.


Introduction
Cancer ranks the top public health threat and is the main cause of death in China [1]. In 2015, more than 4 million new cancer cases and nearly three million deaths occurred in China, in which lung cancer is the most common incident cancer and it is also the highest mortality cancer. Although a series of chemical therapies are available, many of which would accompany with serious side effects due to the chemical toxicity to normal cells. Additionally, drug resistance during the cancer treatment is of great challenge to be overcame. These challenges urged researchers to develop novel antitumor drugs.
Natural products are treasure in the process of drug discovery. Artemisinin (ART), a sesquiterpene lactone natural product, isolated from the traditional Chinese medicinal herb Artemisia annua L., is a powerful antimalarial drug. In 1991, Scholars reported that artemisinin derivatives showed inhibitory activity to leukemia P388 cell [2]. Afterward, Henry's group revealed that dihydroartemisinin combined with holotransferrin exhibited powerful inhibitory activity against leukemia cell line and breast cancer cells with minor effects on normal human lymphocytes [3,4]. In the past decades, increasing publications have been disclosed in the area of the artemisinin and its derivatives as antitumor compounds [5][6][7].
Over the years, Artemisinin was widely used as antimalarial drug and potential anticancer lead compound, but the clear target and mechanism are still unknown. Furthermore, there are some pharmacokinetic limitations of artemisinin, such as low solubility in water or organic media, low bioavailability and short plasma half-life in vivo [8]. To overcome these problems, a series of artemisinin derivatives such as artemether, Artemisia ether, dihydroartemisinin, artesunate were designed to improve both antimalarial and anticancer activities [9].
Quantitative structure activity relationship (QSAR)which encompasses several analysis methods such as comparative molecular field analysis (CoMFA), H  comparative molecular similarity indices analysis (CoM-SIA) [10,11], topomerCoMFA [12] and hologram QSAR (HQSAR) [13] was widely used in drug development process. Numerous researches documenting QSAR studies on the antimalarial artemisinin have been published. Cheng and co-workers established CoMFA and CoMSIA models to study artemisinin and its analogues as antimalarial agents. They concluded that the CoMFA and CoMSIA models had a good predictive ability and well matched the docking results [10]. Avery group reported the CoMFA and HQSAR analysis on a series of 211 artemisinin derivatives with known antimalarial activity [14]. Yadavand coworkers employed QSAR and molecular docking to screen the novel antimalarial artemisinin derivatives [15]. By contrast, there are few QSAR studies on the anticancer property of artemisinin. In this work, traditional QSAR, topomer QSAR and HQSAR based on 46 compounds against A549 cell were performed to get useful information to guide the structural modification of artemisinin for improvement of its antitumor activities.

Data Sets
An array of artemisinin derivatives in Table 1 with reported anticancer activities against the non-small-cell lung adenocarcinoma A549 cells were selected [16,17]. The molecular structures were prepared by D.S 4.0 (Discovery Studio 4.0 Client). Those 3D structures were optimized by energy minimization with SYBYL-9 2.0. The IC 50 values in units were converted in logarithmic unit (pIC 50 ). The relationships between the chemical structures and pIC 50 against A549 were listed in Table 1. During the development of QSAR models, the rational selection of training and test sets have an impact on the predictive ability of the established model [18]. The marked nine compounds by * in Table 1 were selected as the test set.

Structure Preparation and Alignments
The results of CoMFA and CoMSIA analysis depend on the molecular alignment [19]. It is vital to choose a template compound at this stage [20]. The geometry optimizations of the ligands were performed by Powell conjugated gradient algorithm method using the TPIPOS force field until the root-mean-square (RMS) deviation reached 0.005 kcal/mol. The Gasteiger Hückel approach was applied to calculate partial atomic charges. In this study, a ligand-based alignment was adopted since the receptor is unknown. Compound 14, the most potent inhibitor among those selected ligands, was chosen as the template which all molecules were aligned based on. The alignment of structures was shown in Fig. 1.

CoMFA and CoMSIA Methodology
In CoMFA, the energies parameters of steric and electrostatic of the CoMFA were set by the default probe, i.e., a  sp3 carbon atom and ? 1 charge as steric and electrostatic probe respectively. And Tripos force field with a distancedependent dielectric constant at all intersections with a fixed grid spacing of 2 Å was used. The maximum steric and electrostatic energy thresholds were set to be 30 kcal/mol. Firstly, the partial least squares (PLS) analysis was performed based on the CoMFA descriptors as independent variables and pIC 50 values as dependent variables. The leave-one-out (LOO) cross-validation method was used in regression analysis. The column filtering was set to be 3 kcal/mol to improve the signal-noise ratio by omitting those lattices points which energy variation was below this threshold. Then the cross-validation was performed to obtain the optimum number of components (ONC), standard error of predictions (SEP) and crossvalidation squared correlation coefficient (q 2 ). The optimum number of components (ONC) was used for the noncross validation PLS analysis to build the final model with the corresponding conventional correlation coefficient (r 2 ), standard error of estimate (SEE)and the F value. Five CoMSIA fields including steric, electrostatic, hydrophobic, hydrogen bond donor and hydrogen bond acceptor have been evaluated using the default probe, charge, and grid spacing which were similar to the CoMFA model. The column filtering was set to be 3 kcal/mol and the attenuation factor was set to be 0.3. The subsequent parameters were obtained in the same way as described above in CoMFA.

TopomerCoMFA
TopomerCoMFA, an objective QSAR methodology, accelerates lead optimization [12]. Being different from the traditional CoMFA, topomerCoMFA is a No-Alignment CoMFA [21]. In topomerCoMFA, electrostatic and steric fields were calculated for the fragments from 3D structural cleavage.

Hologram Quantitative Structure-Activity Relationship
Compared with the traditional 2D fingerprint, HQSAR encodes more information including branched and cyclic fragments as well as stereochemistry. The hologram was constructed by the HQSAR descriptors which encoded topological and compositional molecular information [22]. In general, the best HQSAR model was obtained by various hologram lengths and the fragment. In order to instruct the predictive ability of QSAR models, the predictive correlation coefficient (r 2 ) was calculated using the following equation based on the test set: In this equation, SD indicates the sum of square of the difference between the true values of the test set and the mean data of the training set; PRESS is the sum of square of the difference between the predicted and the observed activities of the test set.

Results and Discussion
In this work, 37 artemisinin derivatives were randomly selected as a training set to generate QSAR models. The rest of molecules were used to validate the model as a test set. For an ideal model, the following values with regard to statistics parameters should be required. The predictive correlation coefficient r 2 of the test set must be greater than 0.50 and the cross-validation squared correlation coefficient q 2 of the training set must better than 0.5. In addition, Pearson correlation coefficient R 2 values and standard error of estimate (SEE) as well as the former quantity vary in the range 0-1, of which 1 represents a perfect model and 0 means a model without any predictive ability. The optimum results are listed in Table 2. The statistics parameters q 2 , R 2 (training set) and r 2 (test set) are 0.567, 0.968 and 0.991 respectively, which indicated that the CoMSIA represented the best model. PLS analysis was performed using the training and test sets with pIC 50 values. Debugging the column filtering of 3 led to q 2 = 0.547 and optimum number of components = 7. The statistics parameters of CoMFA model were listed in Table 2, in which a high R 2 (0.98) and a low SEE (0.122) were obtained from the training set and an estimated F-ratio value of 199.966. In the model, the steric and electrostatic fields contributed 47.2 and 52.8% respectively. These data showed that the CoMFA model is reasonable.
A series of models were established by random combination of five fields such as electrostatic (E), steric (S), hydrophobic (H), hydrogen bond acceptor (A) and hydrogen bond donor (D). Among these models, four fields (S/E/ H/A) were selected to establish the best CoMSIA model. The statistics results in Table 2 showed that the crossvalidation correlation coefficient q 2 = 0.567, optimum number of components = 5, and the SEP = 0.547. The non-cross-validated coefficient R 2 of 0.968 with a low SEE of 0.148 was obtained. The contributions of S, E, H and A were 9.5, 36.3, 31.2 and 23% respectively. These values indicate that a reliable CoMSIA model was constructed successfully.
In topomerCoMFA, the first step is to split molecules into fragments from the rotatable bonds. The best data was obtained by splitting the C-O bonds at the C10 of training and test sets, as shown in Fig. 2 when compound 14 was chosen as template. Test set was selected in the same way as described in CoMFA and CoMSIA models. The statistics results were listed in Table 2 with q 2 = 0.559 when the optimum number of components was set to be 7 and R 2 = 0.921 with a low SEE of 0.24. The r 2 of 0.819 for test set was obtained by the equation as mentioned above. The corresponding predictive activity of the training and test sets and contribution values of each fragment were exhibited in Table 3.
The training and test sets were selected in the same means for CoMFA, CoMSIA and topomerCoMFA models. QSAR analysis was performed by screening the hologram lengths of 97, 151, 199, 257, 307 and 353, using the fragment size of 5-8 to select the best model based on the least standard error. The best results were summarized in Table 2, in which a cross-validated q 2 = 0.527 with 6 optimal components and a hologram length of 151 were obtained. The conventional R 2 value and standard error of estimate were 0.921, 0.238 respectively.

Validation of the QSAR Models
The statistics results of the best models for these four QSAR methods were collected in Table 2. The experimental data, prediction values and its residues of the training set of the CoMFA, CoMSIA and HQSAR models were summarized in Table 4, and the rest of 9 molecules in the test set were used to verify the predictive capacity of these models as shown in Table 5. The relevant data of the best topomerCoMFA was exhibited in Table 3. From the Table 2, rational r 2 values were obtained to validate the reliability and predictability of the CoMFA, CoMSIA, HQSAR and topomerCoMFA models (0.787, 0.991, 0.743 and 0.819 respectively). The Fig. 3 graphically explained the correlation of predicted and experimental data.
These four QSAR models showed that there were good linear relationships between the true values and predicted values for the group of the training set and the test set (Fig. 3). In comparison, the dispersity of the HQSAR is distinctly higher than the other QSAR models, indicating that 3D QSAR methods (CoMFA, CoMSIA, and topo-merCoMFA) are more suitable to provide crucial structural modificative information than 2D QSAR (HQSAR). Meanwhile, the CoMSIA model has better predictive ability than CoMFA and topomerCoMFA models observed from the Table 2 and the Fig. 3. It is assumed that four fields (A/H/S/E) from the CoMSIA model are more comprehensive than two fields (S/E) from CoMFA and topo-merCoMFA models.

Graphical Analysis of QSAR Models
In this work, contour maps were calculated using the PLS analysis (StDev * Coeff) and the contour plots of 3D QSAR models with compound 14 as the template are   Fig. 4. In this section, the molecular skeleton was labeled as R 1 , R 2 , and R 3 regions for the sake of contrastive analysis.
In the steric contour maps, the green contours (80% contributions) represent that bulk group is a favorable substituent and the yellow contours (20% contributions) show that bulk group is an unfavorable substituent at this position. In the CoMFA model, a big green polyhedron at the R 2 portion means that big groups are useful to improve the activities, which matches the case that when the groups like Br, methoxyl replace the hydrogen atom in this position, compounds 3, 4, 5, 6, 7 have a higher activity than compound 2, and compound 15 is higher than compound 8. Two yellow contours appear in the map, one of which is near the R 3 region and the other is above the C-9 methyl. The steric contour maps of the CoMSIA (Fig. 4e) and the Similarly, in the electrostatic contour maps, the blue contours (80% contributions) means that the positive charged groups are good for the active values while the red contours (20% contributions) indicate that positive charged groups would decrease the activity in this region. In the CoMFA model, a large blue contour near the region of R 2 indicates that the positive charged groups in this position are favorable to the activity. Some red contours appearing at the R 1 and R 3 suggest that the negative charged groups were useful to increase the inhibitory activity in this area. It can be explained that fluorinated and bromine substituent compounds 3 and 9, 12 have a better activity than non-   Fig. 4g, h. The favored and disfavored hydrogen bond acceptor is represented using cyan (80% contributions) and orange (20% contributions) polyhedrons. There are two cyan polyhedrons near R 1 , R 3 and an orange contour at the region of R 2 as shown in Fig. 4g. On the other hand, in the hydrophobic contour maps, the violet regions (80% contributions) means that hydrophobic groups are favorable for improving activity in this place, where as white regions (20% contributions) indicate that hydrophobic groups are disfavored. There is a large violet polyhedron near the R 2 of the compound 14, and it means that the hydrophobic groups in this area are benefic for increasing the bioactivity. It is exampled by Compounds 6, 7 and 15 with alkoxy substituent showing good bioactivity in dataset. There is also a disfavored white contour around the benzene ring which means that hydrophilic groups improve the bioactivity in this region.
The HQSAR model was performed to reach two major objectives listed as following: (a) accurate prediction of the activities of untested compounds; (b) visual display of the contributions of fragments to the activity of the compounds. The atomic contribution maps which display the (a) individual atomic contributions to the molecule's activity by the color of atoms are shown in Fig. 5. Red, red orange and orange atoms represent negative contribution, while the white atoms represent medium contributions to the model. Yellow, green blue and blue exhibit positive contribution. The most and least active molecules as the template were shown as compound 14 and 1 respectively in the contribution maps. The common structural fragments of artemisinin skeleton are displayed by green blue which contribute positively to the model (Fig. 5).

Mode to Optimize Structures
Based on the results of QSAR models, the mode to optimize structures with compound 14 as the template was summarized in Fig. 6. As illustrated in the oval region, the small groups, positive charged groups and hydrogen bond acceptor are favorable for improving activity, while using the bulky hydrophobic groups or hydrogen bond donor at the rectangle place is beneficial to promote the bioactivity of the molecules. In addition, it can improve the molecular bioactivity if using the negative charged groups, small groups or hydrogen bond acceptor at the circular region.

Conclusion
In this study, Four QSAR models (CoMFA, CoMSIA, topomerCoMFA, HQSAR) were successfully established based on a series of artemisinin derivatives with known bioactivity data, which show promising predictive ability validated by the test set. We found that the 3D QSAR models are more suitable than 2D QSAR. Structure-activity relationship information (A/H/S/E) were established from the CoMSIA model with a higher predictive ability (r 2 = 0.991) than CoMFA and topomerCoMFA as well as HQSAR. The topomerCoMFA showed better results than CoMFA and also take much less time. In addition, the 2D atomic contribution maps from HQSAR showed the contribution of individual atoms to the activity of each molecule. Thus, this study will provide reasonable suggestions for design of new artemisinin derivatives with potent antitumor activities.  Fig. 6 The summary of structural modification for artemisinin derivatives from the QSAR models Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.