Exploring differential evolution for inverse QSAR analysis

Methodological concept

Inverse QSAR depends on the derivation of descriptor coordinates for a given model and data set. The goal of the methodology presented herein is finding desirable coordinates in a pre-defined descriptor space (x space) on the basis of a regression function f(x) representing a QSAR model. Confining the search to the applicability domain (AD) of the model translates this task into a constrained optimization problem (COP). The concept of the optimization is illustrated in Figure 1. Newly derived coordinates should be more desirable with respect to pre-defined evaluation criteria than any other data point used to construct the regression model. Accordingly, COP is formulated as follows:

Figure 1. Optimization concept.

A regression function f(x) fits training data to determine new coordinates in descriptor space. Optimized coordinates based on f(x) fall inside the training data range but yield a higher y value than any other data point.

Minimize f(x)

where x ∈ R^d, f: R^d → R is the function to be optimized, g_j: R^d → R is the j-th inequality function, and h_j: R^d → R is the j-th equality function. The i-th component of x: x_i falls into the range [l_i, u_i].

For the purpose of our analysis, the following assignments are made:

x: descriptors;

–f(x): QSAR model;

–g₁(x): AD model;

g_k(x) (k = 2, …, j), h(x): constraints for descriptors;

l_i, u_i: lower and upper bounds of i-th descriptor.

Constraints are applied to descriptors to ensure meaningful value ranges. For example, if the ‘number of heavy atoms’ (x_p) and ‘number of hydrogen bond acceptors’ (x_q) are selected descriptors, a value of five for the former and six for the latter would be impossible for any given data point (compound). Therefore, in order to prevent such settings, an inequality constraint is required and applied: x_q – x_p ≤ 0.

ε Differential evolution

For addressing COP, the differential evolution (DE) algorithm originally introduced by Stone and Price¹⁶ is investigated herein, which has so far not been considered in inverse QSAR. However, given the conceptual simplicity and computational efficiency of DE, the algorithm has been successfully applied to solve optimization problems in other areas of science and engineering, for example, in scheduling of flow shops¹⁷. In addition, for deriving a COP solution efficiently, ε differential evolution (εDE) was introduced by Takahama et al. as an extension of DE¹⁸, illustrated in Figure 2. A candidate vector v for next generation (also called mutant vector) is derived on the basis of three randomly selected vectors:

Figure 2. Evolutionary algorithm.

The steps involved in evolutionary optimization are outlined. First, a candidate v is obtained from three randomly selected individuals by a differential operation. Second, a crossover operation is applied to an individual x_i and the candidate. Third, the evaluation step involves ε level comparison of v and x_i and results in the next individual.

where x_r1, x_r2, x_r3 are different vectors from the current generation and F represents a scale parameter for the difference vector. If the i-th component of v: v_i falls outside the range [l_i, u_i], v is updated as follows:⁴

An exponential crossover operation with probability-based crossover points is applied to v (the probability is called CR). Either x_i or v is selected as x_i+1, the individual for the next generation, following ε level comparison of the corresponding vectors.

ε Level comparison

For prioritizing candidates, given constraints and the optimized function are taken into account. The constraint violation Φ(x) is defined as follows:

where Φ(x) represents the degree of violation, with p set to one. The ε level comparison determines the order between two sets of pair (f(x), Φ(x)):

where t represents a generation in DE. As a decreasing function of t, ε determines the tolerance of constraint violation and ε(t) is determined as follows:¹⁶

where x_θ is the top θ-th individual, T_c determines the generation in which ε(t) becomes zero, and cp the convergence speed. During the optimization, Φ(x) gradually outweighs f(x). In the initial stages of the optimization ε(t) settings enable the selection of diverse candidates but convergence of the algorithm is determined by Φ(x) becoming zero. Accordingly, T_c was set to one herein.

Regression and applicability domain models

For εDE optimization, any regression function can be employed. In this study, support vector regression (SVR)¹⁹ with ν parameter was applied and the AD was defined by one-class support vector machine (OCSVM)²⁰ classification with ν parameter. This parameter ranges from zero to one and defines the upper bound of the fraction of margin error and lower bound of the fraction of support vectors. AD consists of regions where the output of OCSVM is greater than or equal to zero. For both SVR and OCSVM, the radial basis function (RBF) kernel: k(x_i, x_j) = (−γ||x_i − x_j||²) was used. A hyper parameter set {C, ν, γ} for ν-SVR was determined by cross validation of training data on the basis of Q². For OCSVM model construction, γ was set to maximize the variance of the Gram matrix consisting of the kernel function of the training data²¹ and ν was set to 0.01.

Simulation data

Data points on a (x₁, x₂) plane were randomly generated for x₁: [-2 3], x₂: [-1, 4] to yield 50 training and 20 test instances. The corresponding y values were calculated using Mishra’s bird function (https://mpra.ub.uni-muenchen.de/2718/) adding a Gaussian error with a mean of zero and variance of one, defined as:

Three independent trials were carried out with different random number generators. Training and test data sets were plotted on the output domain of the bird function with color-coded true y values (Figure 3).

Figure 3. Simulation data sets.

In three independent trials, simulation data sets were generated. For each trial, training (black dots) and test (blue squares) data are shown with true y values produced by the bird function f(x₁, x₂).

Compound data sets

From ChEMBL²² (version 22), compound data sets were selected using the following criteria: Maximal assay confidence score of ‘9’, interaction relationship type ‘D’(direct), activity standard unit ‘nM’, activity standard type ‘K_i’, and activity standard relation ‘=’. When multiple K_i values were available for a compound, their geometric mean was calculated to yield its final potency value, provided all measurements fell into the same order of magnitude (otherwise, the compound was discarded). In-house implementation of substructure filters for pan assay interference compounds (PAINS)²³ and other reactive molecules were applied to eliminate compounds with potential chemical liabilities. Filtering was not critical for modeling, but active compounds with sound chemical structures were desired. From qualifying data sets, nine activity classes were randomly selected, as summarized on Table 1.

For each compound, 45 descriptors were initially calculated using RDKit (http://www.rdkit.org). These descriptors included constitutional descriptors (e.g., MW, number of rings, number of rotatable bonds), topological descriptors (e.g., Chi and Kier indices²⁴) and partial charge descriptors based on chemical graph’s topology (i.e., maximum of Gasteiger/Marsali partial charges²⁵). From correlated pairs of descriptors exceeding a correlation coefficient of 0.9, only one was chosen. For each activity class, the final number of descriptors (variables) is reported in Table 1. Compounds from each class were randomly divided into equally sized training and test data sets.

Virtual screening

To test the ability of virtual screening (VS) to identify new active compounds from optimized coordinates, ChEMBL (version 22) was used as a screening source. From 1,414,176 unique compounds passing the substructure filters, training molecules used for modeling of each activity class were removed. All remaining ChEMBL compounds provided a large screening source for VS. For screening compounds, descriptors were calculated as described above and the compounds falling inside the AD of each class-specific model were preselected. Active compounds from each activity class not used for training represented true-positive test instances, regardless of their potency values. The calculation of descriptors for more than 1.4 million screening compounds was computationally demanding and exceeded the requirements for coordinate optimization.

For ChEMBL screening compounds including test instances, two VS ranking were generated. First, Euclidean distances to optimized coordinates were calculated. In this case, compound potency was not considered for ranking. Second, pK_i values were predicted for all screening compounds using the class-specific SVR models. The latter calculations were carried out to determine if true positives were highly ranked on the basis of activity predictions. The area under the receiver operating characteristic curves (AUC) was calculated as an evaluation criterion.

Differential evolution for inverse QSAR

Optimization of descriptor coordinates for preferred values of a given model is a central aspect of the two-stage inverse QSAR process, for which currently only approximate solutions exist. Therefore, a more accurate methodology for coordinate optimization is highly desirable, as investigated herein. The evaluation of εDE as an optimization method for this critical task was inspired by previous results obtained for other types of optimization problems where this approach displayed better performance than alternative evolutionary methods, such as genetic algorithms or particle swarm optimization^27–29. Moreover, ε-based lexicological comparison of individual feature vectors makes this algorithm straightforward to apply to problems where several constraints must be balanced, as is the case in inverse QSAR. Studies on simulation and compound data sets were designed to evaluate whether εDE was capable of effectively optimizing coordinates on the basis of a regression function.

Investigating simulation data

For initial proof-of-concept, εDE-based search for optimized coordinates was carried out using simulation data generated as described above.

For the three simulation data sets, SVR models were built yielding optimized parameters {γ, ν, C} of {4, 0.25, 2}, {2, 0.25, 1}, {1, 0.125, 16}, respectively. For all OCSVM models, γ was 1. As reported in Table 2, these SVR models accounted for the output of the bird function in a statistically meaningful way.

Figure 4 shows the different prediction surfaces of the SVR models for the three trial sets. The surfaces of set one and two were overall similar, whereas the surface of set three differed from the others. In each case, however, individual vectors converged at a single point (Table 3) and optimized coordinates were located in regions of highest predicted y values (Figure 4). In set one, for which the SVR model overall best accounted for the bird function, a training data point was found adjacent to the optimized coordinates, which slightly exceeded the largest predicted y value (Table 3).

Figure 4. Optimized coordinates.

For each of three independent trials, optimized coordinates (red squares) and training data (black dots) were mapped on the SVR prediction surface.

For set two, no training data were mapped to local maxima of the bird function, which resulted in a difficult regression scenario. The maximal measured y value in the training data was 29.46 and for optimized coordinates (1.51, 3.16), the predicted y value was 31.14, also slightly exceeding the largest measured y value. However, the predicted value was much smaller than the corresponding true y value of 48.09 (Table 3), due to the inherent regression limitation.

For set three, the maximal true y value within the domain was 56.18 at (-1.59, 0.06). In this case, several training data points were mapped to regions of y values into which optimized coordinates fell (Figure 4), leading to an extrapolative over-prediction of the corresponding y value of 71.94. However, this over-prediction was correctly balanced when the AD of the model was considered instead, leading to a value of 59.40 and a predicted y value for the optimized coordinates of 59.41 (Table 3).

Despite typical regression limitations highlighted by findings for set two and three, the results obtained for simulation data indicated the potential of εDE for predicting optimized coordinates. Importantly, all solutions converged to single vectors representing novel points in simulation data space with large predicted y values falling inside the AD. Taken together, these results indicated the principal potential of εDE for coordinate optimization on the basis of SVR modeling.

Coordinate optimization for compound data sets

Next εDE optimization was applied to different compound activity classes. In each case, SVR models were derived, optimized coordinates determined, and activity values predicted.

For each compound class, optimized coordinates yielded larger predicted pK_i values than any training or test set compound (Table 4), consistent with the methodological strategy. Optimized coordinates fell inside the AD of each model and were proximal to several active compounds. Nearest neighbors of optimized coordinates were mostly predicted to be highly active (Table 4), indicating the presence of smooth prediction surfaces in the vicinity of optimized coordinates.

Prediction surfaces were further characterized graphically by systematically comparing predicted pK_i values of compounds and calculated distances to optimized coordinates. Figure 5 shows the results for two exemplary activity classes, and Figure S1 shows the results for all classes. For set 51 (5-HT1a receptor ligands) in Figure 5, many highly active compounds were located proximal to the optimized coordinates, indicating that these coordinates fell into a well-populated region of activity-relevant chemical space. For set 194 (factor X inhibitors), training and test compounds tended to exhibit higher predicted pK_i with decreasing distance to the optimized coordinates, hence delineating regions of activity progression, which are relevant for compound optimization and exploitation of optimized coordinates.

Figure 5. Activity prediction.

For two exemplary activity classes, predicted pK_i values are related to the scaled distance of the corresponding compounds to the optimized coordinates. Training data (cyan squares), test data (magenta squares), and optimized coordinates (green circle) are shown.

Data set compounds and optimized coordinates were also projected onto PCA plots of descriptor space (Figure 6,Figure S2). These projections revealed that optimized coordinates were central to activity class regions in feature space. Furthermore, Figure 7 shows structures of the three nearest neighbors of the optimized coordinates for sets 51 and 194. In both cases, these compounds were structural analogs. Hence, similarity in feature space corresponded to close structural relationships. Consequently, this would also apply to structure generation from optimized coordinates, which would result in additional analog(s), consistent with the principles of QSAR and inverse QSAR.

Figure 6. Projection of optimized coordinates.

For the two activity classes from Figure 5, training data (cyan squares), test data (magenta squares), and optimized coordinates (green circle) were projected on a principal component (PC) plot derived from training data. For PC1 and PC2, contributions to data set variance are reported in %.

Figure 7. Nearest neighbors of optimized coordinates.

For the two activity classes in Figure 5 and Figure 6, structures of the three nearest neighbors of optimized coordinates are shown and their ChEMBL IDs, scaled distances to the optimized coordinates and predicted pK_i values, are reported.

Virtual screening

ChEMBL compounds were screened relative to optimized coordinates from the nine activity classes and Euclidian distances were determined. Furthermore, pK_i values were predicted for screening compounds falling into the AD of each class-specific SVR model. The VS calculations ultimately led to alternative distance- and potency-based compound rankings. Table 5 reports screening compound statistics and VS results. For each activity class, screening compounds contained large numbers (497–1152) of true positive test instances. Distance-based VS yielded AUC values of at least 0.6 for five of nine activity classes, with a maximum of 0.76. For the four remaining classes, essentially random predictions were observed. Potency-based VS produced AUC values of greater than 0.6 for seven of nine classes, including values above 0.7 for three classes and a maximum of 0.77. Thus, potency-based predictions led to slightly better compound rankings than distance-based VS relative to optimized coordinates, but prediction accuracy was overall moderate at best. Moreover, the true positive ratio among the 30 top-ranked compounds was generally very low for both distance- and potency-based VS. Figure 8 shows exemplary potency prediction landscapes including optimized coordinates as a reference. Figure S3 shows these representations for all activity classes. Highly potent compounds proximal to optimized coordinates were predicted for several activity classes. However, most true positives were not separated from the bulk of ChEMBL screening compounds on the basis of potency predictions. Overall the ability of VS calculations to identify novel active compounds and separate them from false positives was only limited. Thus, although de novo structure generation from optimized coordinates is challenging, it would be difficult to replace the structure generation step in two-stage inverse QSAR with standard VS calculations. However, despite limited prediction accuracy, the VS calculations provided support for the chemical relevance of optimized coordinates because for each activity class, at least few true positives were among top-scoring screening compounds and proximal to optimized coordinates.

Figure 8. Activity prediction for ChEMBL compounds.

For the two activity classes from Figure 5, predicted pK_i values are plotted against the scaled distance of corresponding compounds to the optimized coordinates. ChEMBL compounds (gray squares) and test compounds according to Figure 5 (magenta squares) falling inside the applicability domain are shown. Optimized coordinates are displayed as a green circle.