2.1

Pearson correlation coefficient

To calculate the correlation between variables, our first idea comes from Pearson correlation coefficients (PCCs), which are mainly used to calculate the linear correlation between operational variables¹. In particular, suppose there are two sets of data, A: {A₁, A₂, ‧ ‧ ‧ A_n} and B: {B₁, B₂, ‧ ‧ ‧ B_n}, and the mean of the sample is:

The covariance of the sample is:

So the Pearson correlation coefficients between A and В is:

in which S_A means the standard deviation of sample A, and in the same way, S_B means the standard deviation of sample B. In this paper, we take the bioactive pIC50 as В and 729 molecular descriptors as A, and separately calculate the Pearson correlation coefficients of molecular descriptors on bioactivity.

2.2

Cosine similarity

Cosine similarity (cos) was initially used to calculate the similarity between two vectors by calculating the cosine similarity between two vectors to digitize the similarity between vectors, which related expression is:

Suppose the vector X is {X₁, X₂, ‧ ‧ ‧ X_n} and the vector Y is {Y₁, Y₂, ‧ ‧ ‧ Y_n}, then the cosine of the angle θ between the vectors X and Y is:

Then, the cosine values are normalized for the similarity:

In this case, the closer the similarity is to 0, the more similar the vectors are to each other. Similar to the Pearson correlation coefficient, we took the bioactive pIC50 as separate vectors. Next, select the main variables according to the cosine similarity separately calculated between the 729 molecular descriptors and pIC50.

2.3

Grey relation analysis

Grey Relation Analysis (GRA) is widely used and is very suitable for calculating correlations between data with unknown relationships between them. It refers to the process of unifying the data with a uniform magnitude and then using one set of sequences as the parent and the rest of the sequences as the subsequence. Then determine the degree of similarity in geometry between the parent and subsequence to find whether the two sets of sequences are closely related². The relevant steps are as follows:

2.4

Solutions

In this paper, we first extracted the top 20 main variables using each of the three feature extraction methods mentioned above. Considering that only the correlation calculated by the Pearson correlation coefficient is meaningful whether linear relationship exists between two variables, the scatter plot matrix shown in Figure 1 expressed the screened main molecular descriptors after obtaining the results calculated by Pearson correlation coefficient. The figure which is used to verify the usability of the results clearly shows that variables screened are basically linearly correlated with bioactivity. Then the molecular descriptors that make ADMET have better properties were screened on this basis, which means the desired molecular identifiers that ultimately have a greater effect on the biological activity and ADMET properties were obtained.

Figure 1.

Correlation scatter plot of major molecular descriptors.

We assigned scores to the molecular descriptors obtained by each method according to their ranking, with the first ranking given 20 points and decrease when the ranking goes ahead. Table 1 shows the 20 molecular descriptors which are most significant effective on bioactivity and their corresponding weights judging by the score, from which we can find that the most significant effect on bioactivity is MDEC-23, following by MLogP, LipoaffinityIndex, nC, LipoaffinityIndex, maxsOH, and minsOH. It is known that MDEC-23 represents the edge of molecular distance between secondary carbons and tertiary carbons, which is important in determining the bioactive properties. Also, variables such as MLogP, CrippenLogP are also important for bioactivity. These facts proved that the variables we screened are of strong practical significance.

Table 1.

The 20 variables with the most significant effects on biological activity.

In order to filter out the molecular descriptors that make ADMET have better properties, we supposed a molecular descriptor need at least three good properties to be regard as suitable. In the file “ADMET.xlsx”, the better data of Caco-2, CYP3A4 and HOB are already marked as 1. But for hERG and MN, we can find 0 means the better properties, so we need to normalize the data, turn 0 in hERG and MN to 1. After normalizing, we counted the number of compounds with ADMET properties of 1 for each compound as count, and used count = 3 as the dividing line to divide the data into two parts, and Figure 2 shows the histogram of this result.

Figure 2.

The histogram of count.

We calculate the mean value of each molecular descriptor of each part separately after normalization to get two matrices of size 1×729, and then subtract the two matrices to plot the difference of the two matrices and plot the difference into the difference histogram shown in Figure 3. Figure 3a shows the histogram plotted by the original difference. Considering that we mainly pick out the numerator descriptors with larger difference values as the variables that affect the larger nature of ADMET exist negative values which are not easy to operate. We plotted the graph (b) in Figure 3 with all values taken as absolute values.

Figure 3.

The histogram of differences.

The graph shows obviously that some molecular descriptors have a larger difference between the count ≥ 3 and count < 3 groups. In order to find the better molecular descriptors that affect the properties of ADMET, we plot the values of the differences after taking the absolute values as the scatter plot in Figure 4a. It is obvious to find a lot of discrete points in the high out. In order to pick out all these discrete ones as much as possible without introducing too many molecular descriptors, we screened the difference of 0.11 as the dividing line and drew Figure 4b. All molecular descriptors above the black line in Figure 4b were picked out to obtain a primary sieve of better molecular descriptors affecting the properties of ADMET, and a total of 52 primary molecular descriptors are obtained.

Figure 4.

Primary screening of the main variables affecting ADMET.

Since our first step was only a primary screening of the main molecular descriptors affecting the properties of ADMET, molecular descriptors were sieved out are more than needed. To further reduce the number of variables, the same 20 main variables were extracted using the previously mentioned method. And the molecular descriptors which are most significant effective on biological activity extracted earlier were intersected to obtain the desired final molecular identifiers with a high impact on biological activity and ADMET properties. Theses molecular descriptors consist of Table 2.

Table 2.

Molecular descriptors effect in both ADMET and bioactivity.

3.1

Construct regression model

A BP neural network with 30 hidden layers was first built using the neural network toolbox in MATLAB, the regression model of molecular descriptors with pIC50 is modeled with the Bayesian Regularization training algorithm. This algorithm has good generalization and can deal with small or noisy datasets by taking more time.

Back Propagation Neural Network (BP neural network) is a multilayer feedforward neural network trained by error back propagation using gradient descent algorithm, which belongs to a kind of artificial neural network. It calculates the gradient of the loss function for individual weights by the chain rule, which, unlike the direct calculation, is valid for one layer at a time, and its associated schematic is shown in Figure 5.

Figure 5.

The schematic of BP neural network.

3.2

Construct classification model

Random Forest (RF) is an integrated learning algorithm which consists decision trees and is capable of performing classification and regression. For regression tasks, Random Forest usually takes the mean value of the output of each decision tree as the final result. Compared with decision trees, random forests are more accurate and more stable. The stochastic process of random forest is shown in Figure 6.

Figure 6.

The schematic of random forest.

The decision tree itself can be used to classify and regress the data. The decision tree for regression is mainly divided into cells in the data feature space, and each cell will have a corresponding output. The regression decision tree mainly answers the “yes” and “no” questions, so its structure is a binary tree, in which there are mainly root nodes, decision nodes and leaf nodes, where leaf nodes refer to the output of a specific category.

Figure 7 illustrates the combination strategy in random forest, and the three main methods commonly used are averaging, voting and learning.

Figure 7.

The combination strategy of random forest.

For numerical regression prediction problems, the averaging method is generally used as a combination strategy, and the averaging method has arithmetic averaging and weighting. In this case, we mainly weight the decision tree, and the final prediction is:

in which h_i regards to the weak learning machine, w_i ≥ 0, .

In this paper, the data obtained after feature filtering is fed into a decision tree classifier to train. And we can classify and predict the unknown data by the trained random forest classification algorithm.

3.3

Construct optimization model

Particle swarm optimization’s (PSO) central idea is derived from a group collaboration-based search algorithm developed to simulate the foraging behaviour of a flock of birds¹³. The PSO algorithm designs a massless particle to simulates the birds in a flock, which regards velocity and position as its attributes. Each particle individually records its current individual extreme value by searching optimal solution in the space, and they share their own extreme value to find the global optimal position. By continuously iterating and updating, the global optimal solution is found after the departure stopping condition are met.

Figure 8 illustrates the overall flow of the PSO algorithm. The first step is to initialize the particle swarm. Suppose the swarm size is n, the position of the particle is x_i, the initial velocity of the particle v_i, and the initial adaptability value of the particle P_i. Then set the swarm iterates m times, denotes the best position that the ith particle passes through after d iterations; gbest^d denotes the best position that all particle pass through in d iterations. When , is replaced by . Similarly, when , the current gbest^d is replaced by . The update equations for the particle positions and velocities in the swarm are shown below:

Figure 8.

The schematic of PSO.

where r₁, r₂ are random numbers between [0,1], w represents the inertia weights, and c₁, c₂ represent the learning factors. The above steps are repeated until the maximum number of iterations is reached. The optimal solution of the objective function and the optimal position of the particle x_i are the output.

We use the regression model established by BP neural network, assuming the output of the random forest regression algorithm is f (x). After that, we use the classification model established by the random forest classification algorithm, assuming that the output of the random forest classification algorithm is h(x). So the constrained optimization problem established is as follows:

4.1

The evaluation criteria

The evaluation indexes of the regression model mainly include MAE, MSE and R². The model’s regression effect is better if MAE and MSE is closer to 0 and R² is to 1. MSE and R² are selected as the evaluation indexes to evaluate the three models.

where Pre denotes the predicted value, and Obj denotes the actual value.

where y is the data to be fitted, denotes the mean of the data, denotes the fitted data, SST denotes the sum of total squares, SSR denotes the sum of the regression squares, SSE denotes the sum of squared residuals.

Equations (15)-(18) show the evaluation metrics for the classification problem. In which True Positive (TP) refers to the correct rate of true judgment; True Negative (TN) refers to the correct rate of false judgement; False Positive (FP) refers to the false alarm rate; True Negative (TN) refers to the miss rate. ROC represents the perception curve, and the value closer to 1 means the result is better.

4.2

Experimental results

To facilitate the operation, we copy the attributes pIC50 and the extracted features from the “ERα_activity.xlsx” file into one table, and take the extracted features as independent variables and pIC50 as dependent variables. Then it is input into MATLAB’s neural network toolbox to fit using BP neural network, and the results of the training and test sets obtained at the end of the fitting are shown in Figure 9.

Figure 9.

The result of fitting with BP neural network.

Figure 10 shows the effect of fitting the predicted and true values of the model created by BP neural network, and it proves the satisfactory of the results.

Figure 10.

The effect of fitting the predicted and true values of the BP neural network model.

To visually compare the advantages and disadvantages of the regression prediction effects of the three models, we calculate MSE and R² for the BP neural network and the other two comparison models. The relevant data are shown in Table 3.

Table 3.

The comparison of the three models.

By comparing the three models, it can be concluded that the BP neural network not only has the best fit closest to 1, but also the mean square error can be kept at a low level.

The data obtained from the previous feature extraction is fed into the random forest classifier to obtain its classification results, and the results are plotted in Figure 11. The experimental results of some other methods are also shown in Figure 11 in order to reflect the superiority of the method used in this paper. It can be clearly seen that the random forest algorithm has the best classification results from Figure 11.

Figure 11.

Results of the comparison.

We use the PSO algorithm for modelling analysis. Since the default PSO algorithm in sklearn is to solve the minimum value problem but the optimization problem requires the solution of the maximum value problem, we make slight changes to optimize our constructed constraint. Figure 12 shows the change curve of the fitness of the PSO algorithm during the solution process.

Figure 12.

The change curve of the fitness of the PSO algorithm.

As the result the optimization algorithm get is the optimal value point, not a range that we want, we repeated iterations of the PSO algorithm for many times to get all the optimal results in the fluctuation of each variable maximum and minimum value. These are regarded as the value range fluctuations of the seven molecular descriptors. The seven molecular descriptors optimal value fluctuation graph are as shown in Figure 13.

Figure 13.

The seven molecular descriptors optimal value fluctuation.

The values of the optimal fitness obtained after each run of the particle swarm optimizer are as shown in Figure 13, which is to verify the accuracy of the range obtained. Also, we selected two molecular descriptors, LipoaffinityIndex and CrippenLogP, and plotted the values of fitness and the corresponding molecular descriptors for each particle swarm training in Figure 14, in which red line represents the change in the optimal fitness obtained each train, and the purple and green represent the upper and lower bounds. It is obvious from Figures 13 and 14 that the optimal fitness of the output obtained from each data entry into the particle swarm optimizer varies very little and can be assumed as constant. On the contrary, the corresponding numerator descriptors have relatively large variations. So, after various iterations with the optimal fitness remaining basically unchanged, we can assume that the maximum and minimum values obtained for each numerator descriptor in multiple iterations are an approximate range of values for the variable.

Figure 14.

The values of fitness and the corresponding molecular descriptors.

From the above analysis, it is known that our method is effective in solving for the approximate range of values when the seven molecular descriptors can make the compound biologically more active for inhibiting ERα while having better ADMET properties, and finally we filled in the obtained range of values in Table 4.

Table 4.

The range of the molecular descriptors.