Abstract: Angiotensin Converting Enzyme (ACE) inhibitors have been acknowledged as the first-line agents for the treatment of hypertension and a variety of cardiovascular disorders. The present study aims to develop Quantitative Structure-Activity Relationship (QSAR) models for a series of N- (mercaptoalkanoyl)- and [(acylthio)alkanoyl]glycines derivatives for the prediction of the activity of novel compounds as ACE inhibitors. Multiple Linear Regression (MLR) and Partial Least Square (PLS) analyses were used to establish the QSAR between ACE inhibitory activities and molecular descriptors. The models generated by MLR and PLS analyses were comparable and demonstrated good predictive ability with statistical values of R2 = 0.76, F = 33.35 and s = 0.328 for MLR and statistical significance value of 0.536 and fraction of variance explained = 0.733 for PLS analysis. Cross-validated R2 values of 0.706 (for MLR) and 0.712 (for PLS) were obtained which indicates good internal predictive capacity of the models. Validation through external test set provided R2 values of 0.698 and 0.726 for the MLR and PLS analyses, respectively. The results obtained indicate the importance of steric (Verloop B1, K-alpha 3) hydrophobic (lipole Z component) and topological (Wiener index) descriptors in determining the activity of ACE inhibitors. Also, druggability of the compounds was checked by assessing the Lipinski’s rule of five and no violation to this rule was found. The utility of this study is for the effective design of novel, more potent ACE inhibitors as it provides important information about the geometrical and structural requirements for the interaction of glycine derivatives with ACE.
INTRODUCTION
The Renin-Angiotensin System (RAS) is pivotal in the regulation of blood pressure and electrolyte balance. Angiotensin-Converting Enzyme (ACE) plays a crucial role in the RAS by the production of a potent vasoconstrictive octapeptide angiotensin II which affects peripheral resistance, renal function and cardiovascular structure (Jackson, 2001).
ACE is a chloride-dependent zinc metallopeptidase that contains 1277 amino acid residues and has two homologous domains, each with a catalytic site and a region for binding Zn++. It is non-specific and cleaves dipeptide units from substrates with diverse amino acid sequences. Bradykinin is one of the many natural substrates for ACE, whose inactivation by ACE further contributes to hypertension (Kuoppala et al., 2000). The myocardium of failing heart in patients of heart failure shows increased levels of ACE mRNA, ACE protein and ACE activity (Wright et al., 2008).
Over-stimulation of the RAS causes an increase in blood pressure, blood volume and changes in cardiac and vascular structure. If this continues for a long time, may lead to excessive vasoconstriction, fibrosis and cardiac remodeling. If untreated, there may be a chain of events, progressing through a number of physiological steps, including atherosclerosis to Left Ventricular Hypertrophy (LVH) coronary thrombosis, myocardial infarcts to heart failure in the end (Wright et al., 2008).
Since the development of first marketed ACE inhibitor, captopril, these agents have become the first-line agents for the treatment of hypertension and a variety of cardiovascular disorders including heart failure, left ventricular hypertrophy, post myocardial infarction, chronic kidney diseases (including diabetic and non-diabetic nephropathy) and proteinuria (Sutters, 2008). As a summary of evidence from clinical trials, it is reported that treatment with ACE inhibitors has a beneficial role in patients selected for the treatment of left ventricular dysfunction after Acute Myocardial Infarction (AMI) and in relatively unselected patients with AMI (Latini et al., 1995). Several clinical trials have been performed to study the beneficial effects of ACE inhibitors on diabetes mellitus induced AMI and it was found that apart from the beneficial effects in vascular remodeling, they also reduced recurrent ischemic events after myocardial infarction (Nesto and Zarich, 1998). ACE inhibitors are more effective than any other antihypertensive drug in treating chronic renal diseases, even in normotensive patients (Fogo, 1999). ACE inhibitors are also suggested to be used in patients with diabetic nephropathy. It is reported that ACE inhibitors significantly retard the development of end-stage renal failure compared to other antihypertensives, particularly in patients with proteinuria below the median (Yokoyama et al., 1997). A brief report of a patient with congenital nephrotic syndrome (development of nephrotic syndrome in the first three months of life) of unusual etiology suggested responsiveness to an ACE inhibitor alone (captopril) (Shreedharan and Bockenhauer, 2005). A brief review of literature cited above clearly shows the superiority of ACE inhibitors for the treatment of cardiovascular diseases.
Since its advent, the Quantitative Structure Activity Relationship (QSAR) technique (Hansch and Fujita, 1964) has been widely employed in modeling biological activities as well as ADME/Toxicity properties. QSAR models are mathematical equations which try to correlate the structural and chemical characteristics of drug molecules with their biological activities. Once the relationships are established, the information helps in rationally designing more potent compounds and the predictions of biological activities can be done for many new compounds, as suggested by several researchers (Jamloki et al., 2006; Abhilash et al., 2006; Karthikeyan et al., 2006; Jain and Chaturvedi, 2008). QSAR studies are of great value in the study of enzyme inhibition in determining the mechanism of inhibition and identifying the important active sites in the receptor. These aspects of QSAR studies become of paramount importance in drug design (Gupta and Kumaran, 2006). Recently, QSAR studies have been successfully applied to guide the design of potent new compounds as anti-HIV (Karthikeyan et al., 2007; Thakur et al., 2007) and anti-cancer agents (Moorthy and Trivedi, 2006).
Various N-substituted (mercaptoalkanoyl)- and [(acylthio)alkanoyl] amino acids derivatives have been designed, synthesized and evaluated in vitro and in vivo as ACE inhibitors (Panday et al., 2009). One of the active member of the series of compounds used in the present study is (S)-N-cyclopentyl-N-[3-[(2,2-dimethyl-1-oxopropyl)thio]-2-methyl-1-oxopropyl]glycine (pivopril or pivalopril) having potency lower than that of captopril (Kumar et al., 2010). This prompted us to further explore glycine based ACE inhibitors.
We preferred simple and less error prone 2D QSAR approach to 3D QSAR. 2D is more beneficial compared to 3D QSAR methods, as it neither requires conformational search nor structural alignment (Shen et al., 2002; Hoffman et al., 1999). Moreover, 2D methods also have some intrinsic advantages such as, structural key type descriptors which implicitly encode much chemical information that might otherwise be difficult to explicitly calculate (Rogers and Hopfinger, 1994). Several studies support the superiority of 2D over 3D methods (Hoffman et al., 1999; Bajorath, 2005; Zheng and Tropsha, 2000; Golbraikh et al., 2001). Also, 2D QSAR methods are easily automated and adapted to the task of database searching or virtual screening (Tropsha et al., 1999).
MATERIALS AND METHODS
The experiments were done at Banasthali University, Rajasthan by Anubhuti Pandey as part of her M.Sc. project during January 2009 to June 2009.
Generation and 3-dimensional optimization of chemical structures: The chemical structures of 63 derivatives of N-substituted (mercaptoalkanoyl)- and [(acylthio)alkanoyl] glycines (Appendix Table A1) were sketched using standalone module of Discovery Studio 2.0. Compound 23a and 74a, being racemic mixtures, were excluded. Also compound 74b was avoided due to its uncertain IC50 value. The structures were then, imported into the worksheet of TSAR (version 3.3; Accelrys Inc., Oxford, England).
The series consisted of three substitutions that were defined using “define substituent” option in the TSAR worksheet’s toolbar. The understanding of the molecular properties or ligand-receptor interactions requires not only information on how atoms are connected in a molecule (connection Table) but also on their 3D structure. The structures and their substituents were converted from 2D to 3D using CORINA. CORINA is a method for predicting 3D-coordinates directly from the structural formula of the molecule. It generates by default one low energy conformation for each input structure and all structures that can be expressed in a valence bond notation can be expressed. Stereochemical information is also considered and rings and chains are considered separately (Sadowski and Gasteiger, 1993).
Partial atomic charges of the molecules were calculated using “Charge 2 - derive charges” option. Deriving partial charges is a prerequisite for several structural manipulations, such as optimizing a 3D model and for many property calculations.
Further, 3D optimization of the structures was done using COSMIC (Vinter et al., 1987). The components of the COSMIC parameters are valance terms (including bond potentials, bond angle potentials and torsional potentials) and non-bonded terms (including electrostatic interactions and Van der Waals interactions). Total molecular energies were calculated using COSMIC by summing bond length, bond angle, torsion angle, Van der Waals and coulombic terms for all appropriate sets of atoms. The force field supplied by COSMIC for energy calculations during a flexible optimization, ensures that only the more energetically realistic conformations are considered.
| Table 1: | Molecular descriptors employed in present QSAR study |
Calculation of descriptors and collection into QSAR data table: The experimentally determined IC50 values of the 63 compounds taken from the literature (Suh et al., 1985) (Appendix Table A1) were converted into the negative logarithm (log 1/IC50) and were imported into the worksheet to be used as dependent variables. More than 188 molecular descriptors were calculated for use in the QSAR modeling studies, using TSAR software programs. The calculated descriptors included molecular attributes, molecular indices, atom counts and VAMP parameters (Table 1). The TSAR methodology assumes that a suitable sampling of these structural descriptors provides all the information needed for understanding their biological properties (Klocker et al., 2002; Kovatcheva et al., 2003).
Data set preparation and data reduction: The 63 molecules of the series were randomly divided into training set and test set. Training set consisting of 53 molecules was used for Multiple Linear Regression (MLR) and Partial Least Square (PLS) model development and the test set consisting of 10 molecules was used later to check the predictive power of the developed model.
Initially, the descriptors with constant values were discarded. Before model development, the correlation of descriptors with each other was investigated. Among the detected collinear descriptors (i.e., R2>0.5) one of them which had the highest correlation with activity was retained and the rest were removed from the data matrix of descriptors. This was done with the purpose of getting the descriptors that were less correlated to each other, as high co-linearity among descriptors can lead to statistical instability, over-prediction and also makes mechanistic interpretation difficult. The reduced set of molecular descriptors was used as independent variables in the regression analysis to establish the quantitative structure-activity relationship.
Model development and validation: The Multiple Linear Regression (MLR) was carried out to derive best QSAR models. Various MLR models were generated using biological activity data as dependent variable and selected descriptors as the independent variables. These models were used to quantify the relationship between dependent Y variables and independent X variables. The models generated in the form of regression equation explain the activity data and can be used to predict the activity of new compounds. Positive values of the regression coefficients indicate that the given descriptor contributes positively to the biological activity, where as negative values indicate that the increase in the value of the descriptor lead to a decrease in the activity value. Statistical significance of the regression equations are tested on the basis of conventional regression coefficient (R2) Fischer statistic (F) and the standard deviation (s).
The models obtained were validated using both internal and external validation techniques. Internal validation was performed by applying cross-validation analysis using Leave-One-Out(LOO) method in which one molecule is removed from the training set and the model is recalculated. The predicted activity for that molecule is then compared with its actual value. This is repeated until each molecule is omitted once. Further, as a method of external validation, the activities of molecules in the test set were predicted using the model developed by the training set.
Partial Least Square (PLS) is a multivariate analysis based on principal component analysis. This technique carries the generation of principal components and the multiple regression in a single step and in addition gives maximal correlation between the principle components (independent variables) and the dependent variable. The linear equation obtained as a result of PLS analysis Eq. 2 indicates the relationship between a dependent (activity) variable Y and independent variables X which are latent variables or principal components. Positive and negative values of the coefficients of independent variables signify their positive or negative effect to the biological activity. PLS has been recommended as an alternative approach to enlarge the information content in each model and avoid danger of over-fitting (Kubinyi, 1996). As an approach to check the robustness and the predictive ability of the models generated using MLR analysis, PLS analysis was performed on the same training set of compounds. Models generated using PLS were also validated by LOO method (Vengurlekar et al., 2010).
Assessment of druggability: Knowledge of Absorption, Distribution, Metabolism and Excretion (ADME) provides a mechanistic platform to understand pharmacokinetic and pharmacodynamic relationships of a potential drug compound. There is a significant relationship between the molecular properties and membrane permeability of a compound. With the aim of predicting which compounds are likely to display poor in vivo absorption characteristics, we decided to check the violation of the Lipinski’s rule of five (Lipinski et al., 1997). According to which, poor oral absorption is more likely when a compound possesses more than five H-bond donors, more than ten H-bond acceptors, clog P (calculated log P) greater than five and molecular weight over 500 Da. Computational assessment of these properties was done with the help of “ADME check” option of the TSAR software.
RESULTS
As an initial approach the MLR model was developed using ten independent molecular descriptors retained after data reduction which were independent to each other and were least inter-correlated. Further, in order to improve the statistical quality of the model, the number of independent descriptors was reduced employing backward elimination using t-values. So, the t-values of all the descriptors were checked keeping the stepping to zero and the less significant variables, having lower t-value were deleted from the model.
The above process helped in the retrieval of seven independent molecular descriptors, Verloop B1 substituent-2, dipole moment X component substituent-2, dipole moment Y component substituent-2, lipole Z component whole molecule, lipole Z component substituent-2, K-alpha3 index whole molecule and Wiener topological index substituent-1, having good correlation with the biological activity and minimum intercorrelation. The model developed with seven descriptors showed the statistical values of R2 = 0.564, s = 0.501, F = 15.537 and R2 (CV) = 0.437 (Table 2). These descriptors were used to perform the stepwise multiple linear regression analysis. In stepwise regression, the independent variables were correlated in a stepwise manner with the dependent variable. With stepping 4, Verloop B1 substituent-2, lipole Z component whole molecule, K-alpha 3 index whole molecule and Wiener topological index substituent-1 entered into the model while dipole moment X component substituent-2, dipole moment Y component substituent-2 and lipole Z component substituent-2 were excluded from the model. On deleting these three descriptors we arrived to a model with comparatively improved statistical values of R2 = 0.590, s = 0.444, F = 17.266 and R2 (CV) = 0.500 (Table 2). For the model generated with four descriptors, graph was plotted between the actual activity and the predicted activity of the training set of compounds (Fig. 1). In the plotted graph, it was observed that certain molecules are fitted far apart from the regression line. These compounds were identified as the potential outliers. The outliers may point toward an incorrectly measured experimental value, a different binding conformation, a significant difference in the physicochemical properties, or structural uniqueness (Liao et al., 2004). After deleting six outliers (16, 38, 40, 32, 53 and 56) and cross validation with leave out one row, the best model was generated that included four parameters Verloop B1 substituent-2, Lipole Z component whole molecule, K-alpha3 whole molecule and Wiener topological index substituent-1. The descriptors included in the final QSAR model and their statistical significance is given in Table 3. The t-values for the parameters are greater than 2 which indicate their significance in the model. The correlation between biological activity and parameters used is given in Table 4.
The QSAR equation generated from training set using MLR analysis clearly show the importance of descriptors selected in the present study.
| (1) |
n = 47, r = 0.872, R2 = 0.761, R2 (CV) = 0.706, F = 33.35, s = 0.328, Residual sum of square = 4.516, Predictive sum of square = 5.537 where, X1 is Verloop B1 substituent-2, X2 is Lipole Z component whole molecule, X3 is K-alpha3 index whole molecule and X4 is Wiener topological index substituent-1.
| Table 2: | Statistical test values of the models developed by MLR analysis |
(a) s: Standard deviation, F: Fischer statistic, R2:
Regression coefficient, R2 (CV): Cross-validation regression
coefficient |
|
| Fig. 1: | Actual vs predicted activity for the training set of compounds
(including outliers) derived from MLR analysis |
In the above-mentioned equation, the value of R2 is 0.761 which means that the equation accounts for 0.761x100 = 76.1% variance in biological activity which is clearly quite a reasonable fit.
The cross-validation regression coefficient, R2 (CV) is 0.706 (higher than 0.6) and there is small difference between R2 and R2 (CV) indicating good internal predictive ability of the developed model. These statistical values are given in Table 2.
F is the Fischer statistic, a measure of the probability that the correlation has not occurred by chance. The F-test reflects the ratio of variance explained by the model and the variance due to the error in the regression. The value of F-test (= 33.35) was also found to be significant. Standard deviation is measured by the error mean square which expresses the variation of the residuals or the variation about the regression line. Thus, standard deviation is an absolute measure of quality of fit of the model. It should have a low value for the regression to be significant. The value of standard deviation is 0.328 which is low.
To confirm the soundness and predictive ability of the model, PLS analysis was performed using the same data set. For a well defined problem, both MLR and PLS should have comparable results (Paliwal et al., 2009, 2010). The results of the PLS analysis as shown in Eq.2 were also evaluated on the basis of R2 (CV) and statistical significance of the model.
| (2) |
n = 47, Statistical significance = 0.536, R2 (CV) = 0.712, Fraction of variance explained = 0.733, Residual sum of square = 12.284, Predictive sum of square = 13.233 (Table 5).
| Table 3: | Statistical significance of parameters in the best QSAR model
obtained by MLR |
aThe regression coefficient for each variable in
the equation, bAn estimate of the standard error on each regression
coefficient derived from the jack-knife procedure on the final regression
model, cGives an estimate of the standard error on each regression
coefficient derived from the covariance matrix, dMeasures the
significance of each variable included in the final model |
|
| Table 4: | Correlation matrix showing correlation between biological
activity and parameter used |
| Table 5: | Statistical test values of the model developed by PLS analysis |
The experimentally determined log (1/IC50) values and predicted log (1/IC50) values for the training and test set of compounds are given in Appendix Table A2. Graphs plotted between the actual and predicted activity of training (excluding outliers) and test set, obtained as a result of MLR analysis, are shown in Fig. 2 and 3, respectively. Figure 2 shows that a regression coefficient R2 value of 0.761 was obtained indicating 76.1% variance in biological activity. Figure 3 signifies that R2 value of 0.698 was obtained for the test set by MLR, suggesting good external validation. Similarly, the graph of actual versus predicted activity of training set (excluding outliers) obtained by PLS analysis (Fig. 4) shows that R2 value of 0.733 was obtained. This signifies 73.3% variance in biological activity.
| Fig. 2: | Actual vs. predicted activity for the training set of compounds
(excluding outliers) derived from MLR analysis |
| Fig. 3: | Actual vs predicted activity for the test set of compounds
derived from MLR analysis |
| Fig. 4: | Actual vs predicted activity for the training set of compounds
(excluding outliers) derived from PLS analysis |
| Fig. 5: | Actual vs predicted activity for the test set of compounds
derived from PLS analysis |
Figure 5 depicts the plot of actual versus predicted activity of test set of compounds, derived from PLS analysis. It displays the R2 value of 0.726 obtained by PLS analysis, further suggesting good external validation.
DISCUSSION
As the model suggests, ACE inhibitor activity increases with an increase in Verloop B1 parameter of substituent 2 and lipole Z component of the whole molecule and decreases with an increase in the values of K-alpha3 index for whole molecule and Wiener topological index for substituent 1.
The Verloop parameters (Verloop, 1961) are a set of multi-dimensional steric descriptors that define a box that can be used to characterize the shape and volume of the substituent. These parameters are very important in explaining the steric influence of substituents in the interaction of the drug molecule with the target macromolecule. Verloop B1 is a width parameter defined as the smallest width of the substituent in any direction perpendicular to the length of the substituent. The study suggests that the ACE inhibitory activity is strongly correlated with variations in substitutions at nitrogen, by various alkyl, cycloalkyl, bicycloalkyl, aryl, alkynyl and heterocyclic groups.
The model suggests that there is an absolute requirement of optimum value for lipophilic descriptor, lipole Z component of the whole molecule. Lipole Z component is a directional component of lipophilicity. Any decrease in the value of lipole Z component of whole molecule leads to decrease in biological activity, explaining the importance of lipophilic core and substituents in interaction of glycine derivatives and ACE.
Steric parameter K-alpha 3 index of whole molecule is the third order K-alpha index which can be used to encode the shape of the molecule important for biological activity (Kier and Hall, 1999). As the value for K-alpha 3 index of whole molecules increases, their biological activity decreases.
The fourth descriptor in the regression model is the Wiener topological index (Wiener, 1947). The Wiener index characterizes the “compactness” of a molecule, being larger for extended chains. Here, Wiener index for substituent 1 is negatively correlated to biological activity which signifies that presence of a bulky group at R1 position is detrimental to the activity, as with compounds 73 and 74c.
All the compounds of the dataset showed no violation of the Lipinski’s rule of five (Appendix Table A3) which indicates their potential druggability and good absorption characteristics.
The knowledge on structural requirements for ACE inhibition has been derived from a number of SAR studies (Wyvratt and Patchett, 1985; Mayer et al., 1987). In particular, the following three features required for ACE inhibitor activity has been proposed: (1) a terminal carboxyl group (believed to interact with an arginine residue in the receptor) (2) an amido carbonyl group (with hydrogen bonds to an amide carbonyl in the receptor) and (3) a zinc-binding group (which can be a sulfhydryl, carboxylic acid or a phosphinic acid). This structural information defining the minimal set of active site groups necessary for ACE inhibition has formed the basis for several 3D QSAR studies. These studies have analysed databases of diverse structural classes of ACE inhibitors to determine a common three-dimensional geometry for the active site consistent with their activity (DePriest et al., 1993; Bohacek et al., 1996; Bersuker et al., 2000; San Juan and Cho, 2005).
The model of ACE inhibitor binding suggested by DePriest et al. (1993) indicates the importance of hydrophobic and steric interactions. This is in support to our results. Bohacek et al. (1996) constructed 3D models to design dual ACE/NEP inhibitors. Their study reveals that even in different zinc metalloproteases, the inhibitor atoms between the zinc binding group and P1' carbonyl side chain have a common geometry. Present study is in agreement to this postulate as it also reveals the importance of optimum shape and topology especially in the zinc binding domain of the glycine based inhibitors.
| Fig. 6: | The chemical structure of captopril |
Bersuker et al. (2000) proposed an improved electron-conformational method of pharmacophore identification and bioactivity prediction. According to these researchers, in case of many conformations of the same compound only the one having lowest energy, should be parameterized. In this work, we optimized the structures using COSMIC which ensures that during a flexible optimization, only the more energetically realistic conformations are considered. According to Bersuker et al. (2000), the active fragment (pharmacophore) responsible for ACE inhibitor activity contains three oxygen atoms and one atom X (X = S, N, O). In addition, their study reveals the presence of hydrophobic auxiliary groups (AG) surrounding the entire inhibitor structures which may attribute to an enhanced inhibitor activity. The model developed in our study is consistent with this study as it also suggests that hydrophobic character of entire molecules is important for an enhanced ACE inhibitor activity.
San Juan and Cho (2005) suggested that steric groups are favoured along the proline ring of the 2-methylpropionyl-Lproline and the sulfonyl group of 3-mercapto moieties of the captopril chemical structure. Furthermore, the presence of electropositive group on the propionyl moiety of captopril as well as hydrophobic group at the sulfonyl moiety is expected to increase activity. Present 2D QSAR study, conducted specifically on glycine based ACE inhibitors is in line to the above findings indicating the importance of steric influence of N- substitutions in interaction of N- (mercaptoalkanoyl)- and [(acylthio)alkanoyl]glycine derivatives and ACE. In addition, the developed model also signifies the importance of optimum shape and hydrophobic character of the entire molecules for an enhanced ACE inhibitor activity. Furthermore, larger group or extended chain substitutions at the sulfhydryl moiety in mercaptopropanoyl glycines are unfavourable to the inhibitor activity. This postulate is also supported by the fact that captopril having a small sulfhydryl moiety as the zinc binding group (Fig. 6) has high ACE inhibitor activity (IC50 = 0.017 μM).
In conclusion, the validated QSAR models generated in the present study provide important information concerning the geometrical and structural requirements for the desired biological activity and can be used for the effective design of new, more potent angiotensin-converting enzyme inhibitors.
ACKNOWLEDGMENT
Computational resources were provided by Banasthali University and the authors thank the Vice Chancellor for extending all the necessary facilities.
APPENDIX
| Table A1: | Structure and biological activity data of ACE inhibitors used
in QSAR analysis |
| Table A2: | Actual activity, predicted activity and corresponding residuals
for the training and test set of compounds |
| Table A3: | Assessment of the Lipinski’s rule of five |