To mitigate the risk of overfitting, a random partitioning was applied to the set of 34 compounds at a ratio of 1:3, resulting in 8 compounds assigned to the test set and 26 compounds allocated to the training set. The primary objective of the training set is to establish and refine the model, encompassing the construction, calibration, as well as the identification of key variables and algorithms. Meanwhile, the test set serves primarily for model calibration purposes, ensuring the assessment remains unbiased and does not involve parameter modification. Ultimately, decisions regarding algorithm adjustments or model retraining are contingent upon evaluating the overall fit of the model.

The performance of the QSAR model relies heavily on the appropriate selection of molecular descriptors, necessitating the structural optimization of the compound under investigation. In this study, the chemical structure was initially sketched using ChemDraw (Evans, 2014) and subsequently subjected to optimization using HyperChem (Froimowitz, 1993). The optimization process involved employing the molecular mechanics field (MM+) for the initial optimization, followed by the selection of the AM1 or PM3 model based on the presence or absence of S and P atoms. Furthermore, the structure was cyclically optimized using the Polak-Ribiere method until the root mean square gradient reached a threshold of 0.01. Finally, the CODESSA program (Katritzky et al., 2001) was utilized to compute molecular descriptors encompassing quantum chemistry, structure, topology, geometry, and electrostatic properties.

In the process of constructing linear models, the Heuristic Method (HM) was employed to extract all molecular descriptors. Subsequently, feature selection was conducted to determine the optimal number of descriptors that effectively represent the chemical structure while excluding descriptors with minimal impact. Objective measures, such as the F-test, R², R² _CV, and t-test, were used to evaluate the correlation coefficients between two parameters. Additional descriptors were iteratively added until further inclusion of descriptors had little influence on the results. The linear model obtained through this procedure consisted of six descriptors.

Gene Expression Programming (GEP) is a powerful technique rooted in programming and algorithms (Allen, 1992), surpassing the capabilities of both. Unlike coding numbers or analyzing trees, GEP utilizes linear chromosomes as candidates (Kaydani, Mohebbi, and Eftekhari, 2014). The coding of constant-length linear symbols and the derivation of individual phenotypes, similar to coding codes and expression trees, respectively, are employed (Teodorescu and Sherwood, 2008). The candidate chromosomes are generated from the feature set and the end set, and then encoded into an expression tree (ET) format to calculate the equation (Gharagheizi et al., 2012).

Figure 1 illustrates the overall process of GEP. Fitness functions are applied to a random number of chromosomes, with termination conditions being either the achievement of the predicted value or reaching the maximum number of iterations. When the termination requirement is not met, individuals are selected using the elite roulette approach. Genetic operations such as mutation, transposition, and recombination are applied to the selected individuals to form a new generation. This process is repeated iteratively to obtain improved results.

At the operational level, the molecular descriptor values are inputted using automated problem software (APS), and the GEP technique is employed to derive the nonlinear model. The quality of the model is evaluated using R² as an objective assessment indicator, and appropriate parameter selection is critical to achieve this objective.

Thus, two models were designed: one linear and the other nonlinear. Clearly, the nonlinear model generated by the GEP method outperforms the linear model in terms of both predictiveness and stability. However, it should be noted that 2D-QSAR alone fails to fully capture the three-dimensional relationship between molecular structure and activity, highlighting the need for further exploration through 3D-QSAR studies.

The IC₅₀ denotes the concentration of a compound necessary to inhibit a biological process by 50% (Sebaugh, 2011), and it frequently covers a wide range of magnitudes. In light of this, we employed the formula log (IC₅₀) + 9 to substitute the IC₅₀ values for the 34 compounds. This approach was undertaken to facilitate the analysis and processing of the data, ultimately enhancing both accuracy and stability. These modified data were randomly divided into validation and training sets in a 1:4 ratio.

Furthermore, the ChemDraw structures of the compounds were imported into the SYBYL program for additional optimization. Unlike previous optimizations, this step aimed to minimize the energy of the CoMSIA structure as much as possible. To achieve this, the Tripos force field and Powell gradient technique were implemented within the software (Yu et al., 2015).

This optimization process was underpinned by a systematic calibration of salient operational parameters. The Dielectric Function was set to “Distance,” the MB Cut-off was designated at 8, and the Dielectric Constant was adjusted to 1. Further refinement was achieved by determining the Max Displacement at 0.01, establishing the Minimum Energy Change at 0.05, setting the LS Accuracy to 0.001, configuring the RMS Displacement at 0.001, and specifying the Gradient at 0.05. Following this intricate calibration, the software initiated a comprehensive sequence of 1,000 iterations, leading to the generation of an optimal minimal structure. This resultant structure subsequently served as the foundational 3D conformation for ensuing analytical stages.

The selection and comparison of conformations are crucial as the compound’s structure plays a vital role in subsequent modeling processes (Patel et al., 2008; Ai et al., 2011). Among the 34 compounds, 21E, which exhibited the highest activity value, was chosen as the reference or stacking template. The bold segment of 21E is exceptionally designated due to its manifestation of the identical structure that compound 21E shares with the other compounds. Based on this reference, the remaining compounds were aligned and arranged accordingly (as shown in Figure 2).

The Comparative Molecular Similarity Index Analysis (CoMSIA) is based on the concept that changes in molecular bond affinities are strongly correlated with changes in molecular properties (Yu et al., 2015). This technique involves the calculation of molecular fields as numerous contour plots using a Gaussian function that depends on distance (Li et al., 2012). In contrast to CoMFA, CoMSIA utilizes contour plots to depict the five distinct spatial molecular fields. This approach eliminates the need for arduous procedures, such as aligning the grid with the molecules in the dataset, which is required in CoMFA (Bordas, Komives, and Lopata, 2003). As a result, CoMSIA enhances both clarity and precision in the analysis.

The region for molecule stacking is filled with a cubic grid with a pitch of 2 Å and extending 4 Å in all directions. The default probe generates a multi-molecular field grid, and the affinity is associated with molecular properties using partial least squares (PLS). This leads to the construction of a 3D quantitative conformational model.

The leave-one-out (LOO) cross-validation test is employed to determine statistical significance, providing cross-validated correlation coefficients (Q²) and best group scores (ONC) (Hadni and Elhallaoui, 2020). Additionally, the ONC is further analyzed through non-cross-validation to formulate the ultimate regression model. Objective assessment indicators such as the standard error of estimate (SEE), F-value, and non-cross-validation correlation coefficient (R²) are used for the evaluation of non-cross-validation modeling (Yan et al., 2020).

The method of external validation was employed to forecast the activity of the compound set in the test group (Yan et al., 2020), yielding the correlation coefficient R e x t 2 , which was determined by the following equation: R e x t 2 = 1 − ∑ i = 1 n t e s t y i − y ∼ i 2 ∑ i = 1 n t e s t y i − y ¯ t r 2

In the aforementioned equation, ntest represents the total quantity of compounds present within the test set. Symbol y ¯ t r represents the mean value of the compounds’ activity observed within the training set. Additionally, symbols yi and y ∼ i represent the experimental and predicted values, respectively, of the compounds’ pharmacological activity within the test set. It is worth noting that when the correlation coefficient R² exceeds 0.5, the model demonstrates robustness and exhibits excellent statistical predictive capability (Yang et al., 2011; Mouchlis et al., 2012).

Moreover, we performed an supplementary validation of the model utilizing R m 2 to assess its rationality (Roy and Mitra, 2012).

The formula is expressed as follows: R m o v e r a l l 2 = R 2 * 1 − R 2 − R 0 2

Within the formula, R 2 symbolizes the square of the correlation coefficient between the predicted values and experimental values of all compounds in both the test set and validation set, while R 0 2 denotes the square of the correlation coefficient with a zero intercept. When R m 2 exceeds 0.5, it serves as an indicator of the model’s substantial stability (Pratim Roy et al., 2009).

Meanwhile, an internal validation, specifically, a 20-times Y-randomization validation, was performed to ascertain the optimal model. We subjected the dependent variables to randomization and generated new QSAR models. On each occasion, these newly created models exhibited lower q² and R² values than the original model. Subsequently, the parameter R p 2 was introduced to assess the disparity between the randomized and original model R². The formula utilized is as follows: R p 2 = R 2 * R 2 − R r 2

Where R r 2 represents the average R² value over the 20 iterations of the stochastic model, and a favorable outcome for the model R p 2 should yield a value greater than 0.5 (Rücker, Rücker, and Meringer, 2007).

The compounds underwent initial optimization using the SYBYL software, employing the Tripos force field. Following that, the conformation displaying the minimum energy value was chosen for the molecular docking procedure. The PLK1 receptor (ID: 3db6) was obtained from the RCSB Protein Data Bank (PDB) during the second stage. After removing water molecules and hydrogenated atoms, the protein ligand was extracted and discarded, while retaining the binding site for further analysis.

Following this, a flexible docking approach was employed to facilitate the interaction between the ligand and receptor, allowing for the exposure of the active pocket at the binding site. Docking was performed utilizing Sybyl-Dock as the benchmark, employing a docking threshold with a value of 0.5, an expansion factor of 1, and a retention of 20 conformational alterations. The evaluation of ligand-receptor interactions was objectively conducted using a comprehensive scoring function, with higher values indicating a stronger binding impact of the drug. Undoubtedly, the considerations of hydrophobicity, enthalpy, and polarity are crucial in this context.