{'fraction_non_alphanumeric': 0.03193585337915235, 'fraction_numerical': 0.01809851088201604, 'mean_word_length': 5.086447295036252, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 4, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Development of Predictive in Silico Cytotoxic Activity Model to Predict the Cytotoxicity of a Diverse Set of Colchicine Binding Site Inhibitors

DOI: 10.14744/ejmo.2022.44123 EJMO 2022;6(2):172–181

kaloid, taxane, and colchicine binding site (CBS). [2] In comparison to other tubulin inhibitors, microtubule-targeting agents usually have limitations in developing resistance. CBS agents, on the other hand, have various advantages. The CBS is located between the αand β-tubulin dimers. [4] Colchicine itself is an FDA-approved drug for gout, but it also has implications for different inflammatory diseases such as familial Mediterranean fever pericarditis and Behcet's disease. [5] Many drugs that bind to CBS are less sensitive to multidrug resistance (MDR), which reduces the effectiveness of other inhibitors. [6] The overexpression of class III tubulin, which affects the conformation of the taxane binding site but does not confer resistance to drugs that bind to the CBS, is one of the major causes of drug resistance, suggesting that the CBS may circumvent resistance. [2] Microtubules are important regulators of endothelial cells and vasculature network formation, and it has been shown that CBS is an attractive target for a tumor's established vasculature. [7] Publications related to CBS inhibitors have increased over the last decades. [2,[8][9][10][11][12][13][14][15][16][17][18] Most of the inhibitors were developed to address conventional MDR in the last few years, which gives enough data to build a cytotoxic activity model purely based on known tubulin inhibitors. [1,4] QSAR is the most used chemometrics approach for obtaining data that are useful for drug discovery and medicinal chemistry. It relates a molecule's bioactivity to its physical or chemical properties. Descriptors are numerical values that represent all of the physical and chemical attributes stored with the chemical structure. The IC 50 value represents the response or dependent variable in the QSAR model, whereas these descriptors represent independent variables. In the literature, there are numerous examples of QSAR models being successfully used to screen chemicals for their biological activities. [19][20][21][22] Previously, many QSAR and pharmacophore models were built to predict the biological activity values of colchicine binding inhibitors and to analyze the mechanism of some particular scaffolds. [18,23,24] Some models were developed with noncolchicine binding inhibitors and are based on a set of compounds with few scaffolds. [18] Thus, there is a need for models based on known colchicine binding inhibitors to predict and identify novel colchicine binding inhibitor molecules, which can give a more accurate result.
Tian et al. previously built a 3D-QSAR pharmacophore model using 26 tubulin inhibitors using the HypoGen algorithm; however, there were insufficient data to build a better QSAR model. [25] Similarly, Zhang et al. built a model employing 116 compounds (e.g., arylthioindoles, thiazoles, benzo[b]furan, colchicine, and CA4) as tubulin inhibitors, but did not specify which binding site they were targeting, although the binding site plays an essential role. [18] We built a single QSAR model for predicting the biological activity of heterogeneous compounds that bind to the CBS because earlier QSAR models used only a specific scaffold.
Recent QSAR studies have relied on vast, heterogeneous data sets containing a wide range of structurally varied molecules. A total of 213 chemicals that bind to tubulin's colchicine binding domain were obtained from PubMed literature, which was tested against the A549 cancer cell line. Regression models were developed using multiple linear regression (MLR) based on significant molecular descriptors. Finally, the contribution of several descriptors was studied, and the best prediction model was chosen by comparing model measures, model performance, and accuracy as per the Organization for Economic Cooperation and Development criteria (OECD). Figure 1 depicts a flow diagram of the study's workflow. This includes a large-scale QSAR model for predicting and analyzing the inhibitors of the CBS, which was carried out according to OECD guidelines: (i) a defined endpoint for the data set; (ii) an unambiguous learning algorithm; (iii) a defined QSAR model applicability domain (AD); (iv) acceptable goodness of fit, robustness, and productivity measurements; and (v) a mechanistic interpretation of the QSAR model. Workflow of QSAR modeling, in which is as follows: (1) molecular structures were collected and curated; (2) 2D and 3D optimization and descriptors generation and assigning pIC50 value to the respective structure; (3) the data set was divided into training and test set; (4) built model on training set data by deleting highly correlated data and selecting signification descriptors; (5) the performance of the model is evaluated by R2, Q2, RMSE, and other statistical parameters; (6)an applicability domain (AD) is defined.

Data Collection and Activity Evaluation
The publication related to inhibitors of CBS of tubulin increased in the last few years. This gives us the required data to achieve a predictive QSAR model. A total of 213 colchicine binding inhibitors against the A549 cancer cell line were collected from the PubMed literature of the last 10 years (Fig. 2). All of these compounds were divided into three groups, Sublib-I, Sublib-II, and Sublib-III, [24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40] based on the types of molecules. All of the molecules with their IC 50 values are presented in Table S1-S17, which are available in supplementary material.
Sublib-I: Combretastatin A4 (CA-4) connecting cis double bond bridge and two hydrophobic rings. Several derivatives were reported by structural modification of CA-4. All of the structures with different structural modifications are represented in the following libraries.

Building of Molecular Structures and Descriptor Calculation
All the collected molecules are 2D represented using Mar-vinSketch (https://chemaxon.com/products/marvin) software and are converted to their 3D representation and optimized geometrically using ChemSketch 3D viewer by acdlab (https://www.acdlabs.com/). Each of these compounds had associated in vitro cytotoxic activity values (IC 50 values) against the A549 cell line.
Here, we generate all the 2D and 3D descriptors using PaDEL-Descriptor software including atom-type electrotopological state descriptors, 2D autocorrelations, WHIM, Petitjean shape index, and count of chemical substructures identified by Laggner. File formats were converted to SD format and saved and exported to PaDEL-Descriptor software V2.20, which is an open Source Software available for the calculation of molecular descriptors and fingerprints. A set of 1400 molecular descriptors was calculated for all the 213 molecules.

Data Preprocessing and Selection of Significant Descriptors
The data set is preprocessed, which is the deletion of missing values. Here, our response value is IC50, but as the unit of all IC 50 values is in the micrometer range, we convert IC- 50 to pIC 50 (-log IC50) because higher values indicate exponentially greater potency, so pIC 50 becomes our Y (dependent variable). We use all the other features as our inputs x (independent variables). The calculated molecular descriptors (x) were normalized by a technique that preserves range (maximum and minimum) which scales the data, so every feature sits in the range 0-1.
The selection of relevant descriptors is most important because the process reduces the number of input descriptors when developing a predictive model which reduces the computational cost of modeling and improves the performance of the model. In the first phase, any parameter that was not calculated (missing value) for any number of compounds in the data set was rejected. Some of the descriptors were discarded because all of the compounds had a zero value (zero tests). To reduce the influence of collinearity and eliminate redundancy, a correlation matrix was created with a cutoff value of 0.9, and variables that showed exact linear dependencies between subsets of the variables and multicollinearity were removed from the analysis (high multiple correlations between subsets of the variables). A systematic search was performed in the order of missing value test, zero tests, multi co-linearity, and descriptors were discarded by applying a genetic algorithm (GA) to select the best significant descriptors toward the biological activity value. [41] There are numerous examples in the literature where MLR-GA methods were used successfully for the selection of a significant set of descriptors as powerful search techniques based on the evolution of biological systems for QSAR modeling. [42][43][44] GA is a natural variable selection method that belongs to the class of evolutionary algorithm techniques inspired by natural evolution. Here, we set the parameters: number of iterations to 100, equation length (number of descriptors) to 15, mutation probability to 0.3, the initial number of equations generated to 500, and the number of top equations selected in each equation to 30 based on mean absolute error-based criteria.

Building and Evaluating the QSAR Regression Model
The development of a model MLR is an approach that demonstrates a direct link between a dependent variable Y (IC50) and an independent variable x (descriptors). The model is fitted to minimize the sum of square discrepancies between experimental and predicted values of s biological activities. The dependent variable (IC50) Y in regression analysis is dependent on x (descriptors); however, MLR examination expands this concept to include more than one variable, and the regression equation comes in the form: where Y is the dependent variable, bs are standard errors (SEs) of constant, cs are regression coefficients for corresponding x values (descriptors), and a is a regression constant. In the current work, the models were built with the simple MLR method with the selected variables from GA using DTC-QSAR_v1.0.5 (https://dtclab.webs.com/software-tools).
The established QSAR models are judged by the statistical parameters N (number of substances in regression), K (number of descriptors), and R 2 (squared correlation coefficient); Q 2 (cross-validated correlation coefficient), pred R 2 (R 2 for external test set), and F-test (Fischer's value) for statistical significance. The difference in the experimental activity of the data set determined by the regression equation is measured by the regression coefficient R2. A QSAR model is thought to be predictive if the accompanying conditions are fulfilled: R 2 > 0.6, Q 2 > 0.6, and pred R 2 > 0.5. The F-test reflects the proportion of variation clarified by the model, which varies due to the regression error. The F-test has high estimations, indicating that the model is statistically significant. The model's overall strength is demonstrated by the low standard error of projected R 2 , Q 2 , and pred R 2 .

Internal and External Validation of the QSAR Model
The cross-validation technique was used to internally validate the modeled QSAR equation. This method provides more information regarding the predicted reliability of the QSAR equation. In this study, the leave-one-out cross-validation technique was used, and the cross-validated Q 2 cv was evaluated using the expression: where Y tr is the average observed concentration of the training set, Y is the observed concentration, and Y pred is the predicted concentration. The squared correlation coefficient (R2) was calculated to compare the predicted value obtained by the QSAR equation with the observed concentrations from the experiment, but the issue with R 2 is that it always increases as we add more features to the model, even if they are unrelated to the response. Therefore, choosing the model with the highest R 2 is not a reliable approach for choosing the best model. R 2 adjusted for this is also defined, and unlike R 2 , the error predicted by adjusted R 2 will begin to increase as model complexity becomes very high. R 2 adjusted is defined as Because of its diagnostic means of evaluating the model's predictive power, the cross-validated Q 2 cv value is often smaller than the R 2 value of the QSAR model. [45] The best combination of training and test data (i.e., Compd ID, descriptor matrix, and response) was fed into the MLR-plusValidation 1.3 (https://dtclab.webs.com/software-tools) program, which uses the cross-validation method (leaveone-out) and test set validation based on model acceptable criteria to validate the model internally and externally. Golbraikh and Tropsha provided the following statistical properties of the test set for a better QSAR model with strong predictive power [46] : Here, R 2 represents the squared correlation coefficient between observed and predicted activities, r 20 represents the squared correlation coefficient between predicted and observed activities, and k and k' represent the regression slopes passing through the origin.
Here Y pred_test and Y obs_test are the predicted and observed activities of test set compounds and Y traning is the average values of the training set observed activity compounds.

Development of AD
The AD of the proposed model is defined as the area of chemical space compound structure and response where the model can predict accurately. One of the most significant goals of QSAR modeling is to anticipate the activity of new chemical compounds that fall inside the AD of the created model. The accuracy of any QSAR model is dependent on the accuracy of these novel compounds' predictions. The AD of a QSAR model specified the chemical structure space defined by the properties of the molecules in the training set. In this case, the AD is defined by Roy and Kar's basic standardization approach. [47] The descriptors of the training set should ideally follow a normal distribution pattern. Therefore, according to this, 99.7% of the population set remains within 3 standard deviations (SD), and if the associated standardized descriptors of the compound are larger than ±3 SD, the compound may be an outlier in the training set or outside AD in the test set.

Mechanistic Interpretation of Feature Importance
Feature importance analysis can assist in identifying features that are crucial for bioactivity. The dependency between the descriptors and the pIC 50 values is defined by the mutual information (MI) technique using the scikit-learn python mutual information library. [48,49] It is equal to zero if the two variables are independent and a high value means high dependency.

Results
In silico QSAR analysis selecting descriptors based on the GA, a multilinear regression model was developed containing 15 optimum descriptors. The final selected MLR-GA model is:  50 values, as well as the residual between experimental and predicted values. All descriptors with their decreasing order of ranking for dependency toward pIC 50 are shown in Figure 7, which indicates that descriptor WTPT-4 shows a higher dependency for pIC50.   Smallest absolute eigenvalue of Burden Burden modified eigenvalues Negative modified matrix -n 3 / weighted by relative polarizabilities SpMin2_Bhm Smallest absolute eigenvalue of Burden modified matrix -n 2 / weighted by relative mass Burden modified eigenvalues Negative

Discussion
The important step in building a QSAR model is the selection of significant descriptors. The highest-performing model was not necessarily the best. The chemical space is diverse, and the goal of a good machine learning model is to generalize well from the training data. Here, the performance of the model may be less compared with previously built models but it covers all of the chemical spaces without any loss of data. This allows to make predictions accurately the model has never seen. A total 1400 number of descriptors were generated and only 15 optimum descriptors were selected. The whole data set was split randomly into training (70%) and test set (30%). Further, all the calculations are done with training set compounds only. A systematic search performed in the order of missing value test, zero tests, multi co-linearity, and descriptors was discarded by applying a GA to select the best significant descriptors toward the biological activity value.
Previously, a QSAR model was built for CBS using structural diversity of the selected 116 compounds as tubulin inhibitors by Zhang and group using Discovery Studio, who concluded that a new molecule should have strong van der Waals attraction, which affects the activity value. Similarly, the addition of polar groups to the model increases the activity of colchicine. Also, variations of the steric and electrostatic nature of compounds can lead to an increase or decrease in the activity. [18] Here, the model was built using freely available tools, built using 213 molecules of heterocyclic scaffolds taken from PubMed literature from the last 10 years reported as CBS inhibitors. Feature importance analysis by mutual information was used to evaluate the relative importance and contribution of each descriptor to the model. The first descriptor in the model, which has the highest contribution is WTPT-4 defined as the sum of path lengths starting from oxygen. It indicates that the presence of oxygen with different path lengths plays a major role in anticancer activity. [50] SpMin3_Bhp, SpMin2_Bhm, SpMin5_Bhi, and SpMin8_Bhi have a negative contribution to the model and are the descriptors of the absolute eigenvalue of Burden modified matrix of n = 1-8, weighted relative polarizabilities, relative mass, and relative first ionization potential. The first ionization potential is an atomic property that represents the outermost electronic state. It is an estimation of the energy needed to remove one valence electron from a neutral atom. It indicates that the presence of a less electronegative atom or group increases the anticancer activity, as a less electronegative atom or group requires less ionization potential. [51] SpMax4_Bhv is the largest absolute eigenvalue of Burden modified matrix n = 4, weighted by relative van der Waals volumes. SpMin2_Bhs is the topological distance matrix descriptor. Increasing the relative van der Waals volumes affect the anticancer activity. [52] Similarly, VR1_DzZ, ATS6v, and MATS2c have a negative contribution to the model and are the descriptors of the absolute eigenvalue of Burden modified matrix of n = 1-8, weighted by relative polarizabilities, relative mass, and relative first ionization potential. The descriptors are the index of molecular branching with the smallest value corresponding to the chain graph and the highest value to the most branch graph. [52] As the descriptors are negatively related to the anticancer activity, it may be concluded that decreasing the complexity of the molecule may increase the anticancer activity. [53,54] SPC-5 is defined as a simple path cluster of order 5 that has a positive impact on the model. [55] VE3_D is described as the logarithmic coefficient sum of the last eigenvector from the topological distance matrix. This descriptor encodes electronic, topological, and other geometrical aspects of the compound. The presence of the descriptor in the QSAR model indicates the role of steric and electronic interaction influencing anticancer activity. It utilizes the molecular linear free energy relation. [56] AATS3e is negatively correlated to the model defined as average Broto-Moreau autocorrelation-lag3 weighted by Sanderson electronegativity. It is distributed along with the topological structure and calculated at a given spatial lag usually ranging from 1 to 8. It may be concluded that surrounding electronegativity decreases the anticancer property. [53] MLFER_L is defined as the overall or summation of solute hydrogen bond acidity that has a positive contribution to the model. [57] The model with selected descriptors passed the acceptability criteria proposed by Golbraikh and Tropsha (Table 2). Regression statistics, that is the p-value of the model, suggests that the coefficient of the descriptors used in this model are statistically significant at a confidence interval of 95%. A standardized approach of the AD defined all the compounds of the training set present within the AD except 6 compounds (3, 56, 27, 36, 69, and 75). However, the difference between observed and predicted values is low. As a result, these compounds can be considered influential in fitting the model performance not necessarily outliers to be deleted from the training data set. Similarly, compounds 21, 39, and 5 show outside the AD, but maximum compounds of the test set present within the AD show the confidence within the defined AD.

Conclusion
In this study, we built a QSAR model for the prediction of IC 50 of the unknown chemical compound for the CBS of tubulin on the training set compounds by using the MLR-GA approach. A set of 213 compounds collected from PubMed literature that binds to CBS of tubulin is used for the model building. From the set of 1400 descriptors, 15 optimum descriptors were selected as highly contributing to the biological activity value. Cross-validation of the model (LOO), Troposha's metrics, and rm 2 metrics validate the internal and external predictabilities of the model developed using training and test sets. The selected MLR-GA model has R 2 train = 0.77372, R 2 adjusted = 0.7895, and R 2 pred = 0.7419. Also, the evaluation of AD shows that the model is reliable to make predictions within the chemical space for which it is developed.

Disclosures
Ethics Committee Approval: The study was approved by the Local Ethics Committee.