An approach to identify new antihypertensive agents using Thermolysin as model: In silico study based on QSARINS and docking

2.1 Thermolysin inhibitors database

The database was collected from reports of the international literature (see Fig. S-1 and Table S-1 of supplementary material). These compounds are reversible competitive inhibitors of Thermolysin enzyme and have reported enzyme inhibition constants (K_i) covering a range of 6.8 × 10⁻⁸ to 430 nM. Due to the inhibition wide range of the compounds, we decided to work with pK_i (−log K_i), which varied from −2.63 to 7.16. The database of 176 experimental compounds was divided into training (133 compounds) and prediction (43 compounds).

2.2 Software

For the calculation of molecular descriptors, DRAGON Software (Talete, 2007) was used. All the families of descriptors (0D - 3D) were calculated excluding constant and near-to-constant variables. Program DRAGON is an application for the calculation of the molecular descriptors, originally developed by Milano Chemometrics and QSAR Research Group, and provides more than 1600 molecular descriptors that are divided into 20 logical blocks.

In order to obtain the MLR model, we used the QSARINS (QSAR-Insubria) software (Gramatica et al., 2013a,b) developed at the University of Insubria. It would be ideal to use all combinations of available descriptors for the calculation of the models; however, the number of combinations is usually so great that it is impossible to calculate all models. Therefore, to reduce the model computation time, all combinations were preliminarily used with only a small number of descriptors per model, in order to explore all combinations of low dimension, using the all subset technique, and, later, a genetic algorithm (GA) method was applied to develop models based on a larger number of descriptors.

2.3 Computational techniques

2.4 Model validation according to OECD principles

Algorithm QSARINS is friendly software that can validate the obtained model using OECD principles, in order to increase the confidence of the predictions. According to these principles, the model should have the following: (1) a defined endpoint; (2) an unambiguous algorithm; (3) a defined domain of applicability; (4) appropriate measures of goodness-of-fit, robustness and predictability; and (5) a mechanistic interpretation, if possible (Worth et al., 2007).

Principle 1 is associated with the endpoint definition where it refers to a physicochemical, biological or pharmacological property that can be measured and, therefore, modeled. The aim of this principle is to ensure transparency in the measurement point predicted by a given model. The endpoint in this work was the Thermolysin inhibition constant (K_i), which characterizes substrate inhibition by a competitive inhibitor (Kakkar et al., 1999), expressed as pK_i.

According to Principle 2, QSAR models can be expressed in the form of unambiguous algorithms, considering that the algorithm model is the way of relating the descriptors of the chemical structure and activity (endpoint of the model) through mathematical models or rules based on the knowledge developed by one or more experts. The used algorithm was a mathematical model of Multiple Linear Regression detailed in the previous section.

William graph was used to define the model application domain (DA) (Principle 3), so that reliable predictions of the model have leverage values below the critical leverage with ± 3 standard deviations. Compounds that are outside these ranges may be considered outliers. In addition, Insubria graph was used for DA. This was very useful to assess the reliability of the predictions of experimental compounds, lacking response, and to compare with the predictions of the database.

Following the Principle 4, the goodness of the fit was evaluated using the coefficient of determination R² and a modified form R²_adj that also evaluated the effectiveness by adding a new descriptor to the model.

In order to control the robustness of the model, we carried out internal validation strategies. First, we used the Q²_LOO criterion applying the technique of cross-validation (CV), i.e., iteratively excluding a compound from the dataset (Leaving One Out, LOO), and then calculating a model with the remaining compounds and making the prediction for the excluded molecule. If the internal predictions were good, the value of Q²_LOO was high and comparable with R², so the model was considered internally stable and robust. However, a disturbance of a single compound in a series of 133 cases was too weak to show the real model robustness. Code QSARINS further includes a stronger technique: Leaving Many Out (LMO), which allowed studying the behavior of the model when a larger number of compounds (in this case 30% of the compounds) were excluded. The model was considered stable if the values of R² and Q² calculated for each iteration of the LMO experiment and their averages (R²_LMO and Q²_LMO) were comparable with the values of R² and Q²_LOO of the model.

In order to demonstrate that the model was not the result of a chance correlation, we applied the Y-scrambling procedure. In this procedure, responses are placed randomly, so there was no correlation with the descriptors and, as a result, the performance of the models should decay drastically. In this case, if the model in the validation process was good, the values of R² and Q² of each iteration, and their averages (R²_YS and Q²_YS) must be lesser and lesser with respect to the values of the model.

In order to assess the model’s ability to predict new compounds, an external validation was carried out. The procedure was performed by applying the model equation obtained with the training set, a predictive dataset, i.e., excluding compounds that were never been used for the calculation model. Model performance was evaluated by different criteria such as RMSE_ext, Q²-F1, Q²-F2, Q²-F3, r²M, Δr²_m, CCC, and Golbraikh and Tropsha methods.

The mechanistic interpretation of the model (Principle 5), although it is not mandatory, it is very desirable. In our work, we developed a general interpretation of the model based on the descriptors selected to describe the Thermolysin inhibition. This could be useful to understand the different factors driving the process.

2.5 Virtual screening

With the aim of performing identification/selection of new compounds as Thermolysin inhibitors, the obtained model was used for the virtual screening of the international database Drug Bank (www.drugbank.ca), which is a database available on the University of Alberta (Canada). These data contain more than 1431 molecules approved by the FDA (Food and Drug Agency of the government of the United States), 133 bioactive drugs approved by the FDA and more than 5212 experimental drugs.

3.1 Getting QSAR-model RLM

Using the DRAGON software (Talete, 2007) we calculate 3224 descriptors (0D-3D). To perform the calculation of the 0-2D descriptors we use molecules representation in the ChemDraw software (saved as ∗.mol files). However, to carry out the calculation of 3D descriptors, geometry optimization calculations should be done; several approaches can be used (González and Helguera, 2003; Goodarzi et al., 2013a,b; Perez Gonzalez et al., 2006). In our work, we performed geometry optimization calculations for each compound considered in this study using the semiempirical method AM1 (Dewar et al., 1985) implemented in Mopac. Given the great number of variables calculated, even when 1503 descriptors were excluded, the number of remaining descriptors was 1721. Thus, the data were subjected to a new filtration process using the QSARINS software (Gramatica et al., 2013a,b) and the descriptors that possessed more than 95% of constant values and had more than 90% correlation were eliminated. After that, only 399 variables were considered in the selection process to develop the model, as shown in Fig. 1.

Close

Figure 1

Diagram of variables filtering process.

Export to PPT

Then, we developed several models with low multicollinearity between descriptors and good correlation with the modeled response. The generated population of models were plotted the mean values of R² and Q²_LOO (with their standard deviations) versus the number of variables of the model (Fig. 2), in order to evaluate the performance of the models against their size and to avoid overfitted models. As can be seen, the values of R² and Q²_LOO rose as new descriptors were added to the model, except for six and seven variables, so we can infer that the added descriptor does not give new information. The highest values of R² and Q²_LOO were achieved by the models of eight variables; however, the models with five variables showed very similar values with three less predictive variables. Based on the principle of parsimony, we selected five variables as optimal dimension to develop the final model, which were used to predict the Thermolysin inhibition.

Close

Figure 2

Model performance according to number of variables.

Export to PPT

The best MLR model obtained with its statistical parameters is shown below:

• pK_i = −23.62 + 0.52xMor07u − 20.77xMor12v − 3.12xR5v₊ + 2.87xR5p₊ + 23.58xB01[N–O]

• N = 133 R² = 0.714 R²_adj = 0.702 s = 1.230 F = 63.248

• R²−R²_adj = 0.011 LOF = 1.688 K_xx = 0.231 Delta K = 0.100

• RMSE_tr = 1.202 MAE_tr = 0.975 RSS_tr = 192.024 CCC_tr = 0.833

The model presented an R² of 0.714, so it manifests a proper fit for modeling Thermolysin inhibition. In addition, an R²_adj of 0.702, which is indicative of the convenience to add a new descriptor to the model and, together with the low value of the LOF parameter of 1.688, we can say that no existing overfitting is in the model, as it presents a good fit with minimum number of descriptors. The correlation between the model descriptors is low because the value of K_xx is small (0.231); this allows us to assume that the model has very little redundant information in the selected descriptors. Besides, according to the Delta K parameter (0.1), the model has an appropriate correlation between descriptors and the modeled response, with a small error on the training calculations and estimated parameters (RMSE_tr = 1.202; MAE_tr = 0.975; s = 1.230). Yellow circles in Fig. 3 represent the values predicted by the model equation versus the experimental pK_i values for the training set. This scatter plot permitted detection of systematic trends or grouping of the data and potential outliers therein.

Close

Figure 3

Scatter plot of experimental pK_i versus predicted by model equation.

Export to PPT

As can be seen, the compound NFpLA (compound 43) of the training set has a residual of −4.38 and can be considered an outlier compound. Keep in mind that this is the compound that has the highest pK_i in the database (pK_i = 7.17), and the pK_i predicted by the model is 2.786, which leads to obtain a very high residual value. However, the model predictions were not affected, since a compound with predicted pK_i greater than 2.5 is considered very active. The experimental pK_i for Thermolysin and the results predicted by the MLR model for the training set are shown in Table 1.

3.2 Model validation according OECD principles

In order to validate our model, we followed the OECD regulatory principles to ensure their validity, checking the model performance, i.e., the fitting, stability in the cross-validation and the ability to predict new compounds (Worth et al., 2007). The fitting and stability of the model were evaluated using internal validation procedures; first, we take into account the parameters of the cross-validation Leave-One-Out (LOO). According to the obtained results, it is possible to affirm that the internal predictions are good since the variance explained in the prediction by LOO (Q²_Loo = 0.6868) has a comparable value with R² = 0.7135. Then, the model is considered internally stable or robust with a small error in the predictions (RMSE_cv = 1.2563 and MAE_cv = 1.0190). Fig. 4 shows the values predicted by LOO versus the experimental values of pK_i (clear circles for the training set and dark circles for the prediction set). As in the model predictions, the NFpLA compound is the farthest from the line and is supported by the hypothesis of which it can be considered an outlier compound.

Close

Figure 4

Scatter plot of experimental pK_i versus predicted by LOO.

Export to PPT

A stronger technique included in the QSARINS is Leaving-Many-Out (LMO), which was developed leaving out the 30% of the dataset to study the behavior of our model. This technique was very useful in our database because in opposite to LOO, where for each single disturbance of the data the influence on the predictions was very weak; in LMO technique, each perturbation of data is significant. The model is considered stable because the R² = 0.7135 and Q²_LMO = 0.6810 values are comparable, and calculation in each iteration of LMO and their averages are comparable to the values of R² and Q²_LOO of the model. Fig. 5 shows the Q²_LMO versus the correlation between descriptors and Thermolysin inhibition (K_xy) in which it is observed that the values of Q²_LMO are very similar to each other, which confirms that the model has a good fit, robustness and stability.

Close

Figure 5

Plot of LMO validations compared with the original model.

Export to PPT

In order to demonstrate that the model was not the result of a casual correlation, the Y-scrambling procedure was used, placing the answers at random, so that there was no correlation with descriptors and, as a consequence, the model performance decayed dramatically. The values of R² and Q² of each iteration and its averages (R²_Y-scr and Q²_Y-scr) give us the criteria that the model is good, since these parameters are ever lower with respect to the values of the model (R²_Y-scr = 0.039 and Q²_Y-scr = −0.055). The R²_Y-scr and Q²_Y-scr values against R² and Q² of the model are shown in Fig. 6. Notice that the values of R² and Q² of the model are very far from the values obtained for those parameters in the Y-scrambling experiment, which indicates that the model was not obtained as a result of a random correlation.

Close

Figure 6

Plot of Y-scrambled model compared with the original model.

Export to PPT

The predictive ability of the model was tested using a series of external predictions, which were never used to obtain the model (external validation). Using this procedure, we checked the ability of the model to predict new compounds and their parameters showing results equivalent to the value of model R² (R²_ext = 0.723), RMSE_ext = 1.182, MAE_ext = 0.91, PRESS_ext = 60.091, Q²-F1 = 0.713, Q²-F2 = 0.713, Q²-F3 = 0.723, CCC_ext = 0.817, r²_m_aver = 0.619, r²_m_delta = 0.126.

Where R²_ext is the external determination coefficient (Golbraikh and Tropsha, 2002), RMSE_ext is the Root Mean Square Error in external prediction; MAE_ext is the Mean Absolute Error in external prediction; PRESS_ext is the Predictive Residual Sum of Squares (external validation); Q²-F1 (Shi et al., 2001), Q²-F2 (Schüürmann et al., 2008), and Q²-F3 (Consonni et al., 2009) are the variances explained in external prediction; CCC_ext is the Concordance Correlation Coefficient (Chirico and Gramatica, 2011, 2012; Lin, 1989) and r²_m_aver and r²_m_delta are the Roy criteria: average and delta (Ojha et al., 2011). The predictions of compounds in the external set can be seen in Fig. 3 (blue circles).

A crucial problem in chemometrics and QSAR studies is the definition of the applicability domain (AD) of a classification or regression model. Even with a robust, significant and validated QSAR model, unreliable predictions of the modeled property for the entire universe of chemical compounds can be expected. Indeed, only predictions for compounds that fall within the AD can be considered reliable and not the extrapolation of the models (Dimitrov et al., 2005).

In this study, we used the approach of leverage (h) and standardized residuals described in the technical literature (Atkinson, 1985). Fig. 7 shows the William graph of the model, where open circles represent the compounds of the training and dark circles represent the prediction set.

Close

Figure 7

Williams plot. Hat diagonal values versus standardized residuals.

Export to PPT

As shown in this figure, most of the compounds are within the applicability domain of the model. There is only one compound (compound 34) with the value greater than the critical leverage (h = 0.192), but showing residual within the limits, and one compound (compound 43), having residual off-limits (−3.62) though this is within critical leverage. Both cases must be taken into account as potential outlier compounds.

Besides, the QSARINS provided a new approach for determining the model applicability domain based on leverage and model predictions. The Insubria graphic is very useful to evaluate the position of the molecules lacking experimental response, compared with the structural applicability domain, and to compare their predictions with the compounds having the experimental values. According to the graph of Insubria (Fig. 8), only one compound is outside the AD by leverage (compound 34) and coincides with the William graph.

Close

Figure 8

Insubria graph. Hat diagonal values versus predicted data.

Export to PPT

As shown in the model equation, using only five predictors (molecular descriptors), more than 70% of the inhibitory activity of Thermolysin can be described. Most of the variables are three-dimensional (Mor07u, Mor12v, R5v+ and R5p+) of 3D MoRSE and GETAWAY families, respectively. This result is logical, since there are reports that the enzyme is stereospecific (Jin and Kim, 1997), so the three-dimensional structure of the inhibitors is essential for access and interaction with the active site and subsite pockets of the enzyme. Note that descriptors weighted by the van der Waals volume have a negative coefficient that could be related to the steric hindrance to access the active site inside. In addition, the descriptor weighted by polarizability has a positive coefficient, which reinforces the idea that electrostatic interactions with residues near the active site of the enzyme are crucial to their catalytic activity. Furthermore, the model has a two-dimensional (2D) descriptor (B01 [NO]), based on the search for similar patterns in the database compounds. These suggest that the presence of very close nitrogen atoms and oxygen atoms permits the interaction with the active site of the enzyme.

3.3 Molecular docking

First, the method of the docking to reproduce the crystal structure of the Thermolysin-inhibitor complex (Fig. 9) was made. This test allowed measuring the capacity of the docking to reproduce known structures, which were contained in a crystal.

Close

Figure 9

Conformational comparison of compounds 20, 21 and 22 (in yellow) from the crystal structures in Thermolysin–5H9^a inhibitor complex (inhibitor in red). ^a: N-[(S)-({[(benzyloxy)carbonyl]amino}methyl)(hydroxy)phosphoryl]-L-leucyl-4-methyl-L-leucine-(6S,9S)-1-(benzyloxy)-6-isobutyl-9-neopentyl-4-oxido-1,4,7-trioxo-2,5,8-triaza-4-phosphadecan-10-oate (C22 H36 N3 O7 P).

Export to PPT

As it is shown in Fig. 9, the docked structures had an adequate orientation in the pocket of Thermolysin. It reproduces in an acceptable way the X-ray crystal structure of the Thermolysin-inhibitor complex. The root mean square deviation (RMSD) of the analyzed molecules regarding the crystal structure of the inhibitor was <2.0 Å in the 73.3% of all the analyzed cases (Table 1). In this analysis, it was considered that 2.0 Å is the reference value that identifies an either correct or incorrect resolution of the docking (Gohlke et al., 2000). The RMSD has often been used to measure the quality of reproduction of a known (i.e. crystallographic) binding pose by a docking method (Gohlke et al., 2000). Several authors have suggested that a low RMSD with respect to the true binding pose (in our case the crystallographic structure of Thermolysin-inhibitor (5H9) complex), is good (ideally less than 1.5 Å, or even better, less than 1 Å). This represents good reproduction of the correct pose. Beware reading much significance into relative values of RMSD’s greater than 3 Å. A binding pose with RMSD 4 Å is not better than one of 6 Å. They are both equally poor (i.e. wrong) (Gohlke et al., 2000; Liebeschuetz et al., 2012; Maiorov and Crippen, 1994). Only the compounds 20, 29, 31 and 66 had values higher than 2.0 Å.

Analyzing the docking results, it is worth mentioning that the compound 172 presented the lower union energy (−8.7 kcal/mol) of all the analyzed molecules. As Fig. 10B shows, the orientation of one of the OH groups makes that it coordinates with the metal (Zn²⁺) in the active center of the protein. Besides, one of the terminal benzenes is located in a hydrophobic pocket, formed by Leu133, Val139, Phe130 and Leu202, whereas the other benzene terminal is located in a hydrophobic area that has been formed by the amino acids Trp115 and Tyr157, which stabilize the molecule. The terminal OH groups are oriented toward a polar area that has been formed by the Asn227. The Asp226 is also located in this area. It is anionic and it provides density of charge in the stabilization of these groups.

Close

Figure 10

Representation of the interactions of the compound 43 (A) and 172 (B) in the active center of the Thermolysin. The interaction surface with the solvent has been shaded. The polar area of interaction has been represented with blue ribbons; the hydrophobic area, with green ones, and the density of negative charge, with red ones. The interaction with the Zn²⁺ has been represented with a gray line.

Export to PPT

Of all analyzed compounds, two of them were identified as possible outliers (Compounds 34 and 43) from QSAR method results. In this sense, we made experiments of molecular docking of these ‘outlier’ compounds and a group of compounds with similar structural characteristics, since the model provided three-dimensional information of molecules, and it would be interesting to observe their behavior in the active site of Thermolysin. In order to analyze the actions of compound 43 in the active site of the enzyme, besides, seven compounds were considered by the similarity therewith; they are shown in Fig. 11.

Close

Figure 11

Structurally similar compounds to compound 43 proved with docking.

Export to PPT

Of all compounds in the dataset, the compound 43 is the most active (pK_i = 7.17); in other words, this compound is an excellent inhibitor of Thermolysin. This molecule behaved in a similar way in the docking experiments, with low energy of union (−8.0 kcal/mol), compared with other evaluated compounds (see Table 1) that have similar structural characteristics. This molecule is oriented in the active center of Thermolysin, in such a way that the terminal benzenes are located in a hydrophobic pocket, formed by the amino acids Phe130, Leu202, on the one hand, and Phe114, Trp115, Tyr157 and the Ala113, on the other, which gives stability to the complex. Moreover, the Asp213, when being anionic, stabilizes that area and attracts both terminal OHs, which present hydrogens with a lack of charge, stabilizing even more the complex (Fig. 10A).

According to these results, it is confirmed that this compound is very active. However, the pK_i predicted by our model is similar to those values of structurally related compounds, instead of their real value, so the residual is large. This suggests that the high experimental activity of compound 43 could be a consequence of the difficulty in hydrolyzing the molecule. Note that compounds with lower E_vina (66, 70 and 172) have a hydrolyzable bond near to the phosphorus atom in the chain (P–O– or P–NH–). This enables that the molecule stay longer time in the active site of Thermolysin and, thus, their large value of activity with respect to other structurally similar molecules; for this reason, we decided to keep this compound in the training set.

In the case of compound 34, it has a very close-to-zero (−0.05) residual, so the prediction of this compound by the model is very good. However, this compound is outside the critical leverage, indicating that it may be a structural outlier. As shown in Fig. 12, we can see that it exists in the database other thiols nearby from the structural point of view (compounds 28–33); nevertheless, this compound has the lowest RMSD of all tested molecules (see Table 1) and low energy of union (−7.0 kcal/mol).

Close

Figure 12

Structurally similar compounds to compound 34 proved with docking.

Export to PPT

The low energy of union above is due to its stabilization in the active center of the Thermolysin in a polar area, formed by the histidines 142 and 231, besides the Asn112. This area stabilizes the OH terminal groups (Fig. 13A). Besides, the terminal benzene is stabilized in a hydrophobic area that is formed by Trp115 and Tyr153. Note that the interactions of the compound 34 (Fig. 13A) and compound 32 (Fig. 13B), in the active center of Thermolysin, involve the same residues (His 142, 146 and 231; Glu 143 and 166; Ala113; Phe114; Trp115 and Tyr157). From the foregoing, we decided to maintain also the compound 34 in the training set of the model.

Close

Figure 13

Representation of the interactions of the compound 34 (A) and compound 32 (B) in the active center of the Thermolysin. The interaction with the solvent has been shaded. The polar area of interaction has been represented with blue ribbons; the hydrophobic area, with green ones, and the area with density of negative charge, with red ones.

Export to PPT

3.4 Virtual screening

In this report, the main objective was to identify new compounds with inhibitory activity of Thermolysin, considering that they were predicted by the model as active, were within the AD and possessed suitable pharmacokinetic characteristics for oral administration. For the division into virtual screening, it was considered different pK_i values, which were separated into three groups according to their activity, reported by Klopman and Bendale (1989): inactive (pK < −0522), marginally active (−0.522 ⩽ pK_i < −0022) and active (pK_i ⩾ −0.022). According to that, we decide to carry out a new division for the active compounds: Group 1 (−0.022 ⩽ pK_i < 2.5), Group 2 (2.5 ⩽ pK_i < 3.5) and Group 3 (pK_i ⩾ 3.5). Fig. 14 shows the filtering process carried out in the database Drug Bank. As a result, 141 compounds were predicted by the model as belonging to Group 2 (Table S-2 of supplementary material) and 69 were predicted in Group 3; these compounds and their predicted values can be seen in Table 2.

Close

Figure 14

Scheme of virtual screening.

Export to PPT