3D-QSAR studies on CCR2B receptor antagonists: Insight into the structural requirements of (R)-3-aminopyrrolidine series of molecules based on CoMFA/CoMSIA models.

OBJECTIVE: Monocyte chemo attractant protein-1 (MCP-1) is a member of the CC-chemokine family and it selectively recruits leukocytes from the circulation to the site of inflammation through binding with the chemotactic cytokine receptor 2B (CCR2B). The recruitment and activation of selected populations of leukocytes is a key feature in a variety of inflammatory conditions. Thus MCP-1 receptor antagonist represents an attractive target for drug discovery. To understand the structural requirements that will lead to enhanced inhibitory potencies, we have carried out 3D-QSAR (quantitative structure-activity relationship) studies on (R)-3-aminopyrrolidine series of molecules as CCR2B receptor antagonists.
MATERIALS AND METHODS: Comparative molecular field analysis (CoMFA) and comparative molecular similarity indices analysis (CoMSIA) were performed on a series of (R)-3-aminopyrrolidine derivatives as antagonists of CCR2B receptor with Sybyl 6.7v.
RESULTS: We have derived statistically significant model from 37 molecules and validated it against an external test set of 13 compounds. The CoMFA model yielded a leave one out r(2) (r(2) (loo)) of 0.847, non-cross-validated r(2) (r(2) (ncv)) of 0.977, F value of 267.930, and bootstrapped r(2) (r(2) (bs)) of 0.988. We have derived the standard error of prediction value of 0.367, standard error of estimate 0.141, and a reliable external predictivity, with a predictive r(2) (r(2) (pred)) of 0.673. While the CoMSIA model yielded an r(2) (loo) of 0.719, r(2) (ncv) of 0.964,F value of 135.666, r(2) (bs) of 0.975, standard error of prediction of 0.512, standard error of estimate of 0.180, and an external predictivity with an r(2) (pred) of 0.611. These validation tests not only revealed the robustness of the models but also demonstrated that for our models r(2) (pred), based on the mean activity of test set compounds can accurately estimate external predictivity.
CONCLUSION: The QSAR model gave satisfactory statistical results in terms of q(2) and r(2) values. We analyzed the contour maps obtained, to study the activity trends of the molecules. We have tried to demonstrate structural features of compounds to account for the activity in terms of positively contributing physicochemical properties such as steric, electrostatic, hydrophobic, hydrogen bond donor, and acceptor fields. These contour plots identified several key features, which explain the wide range of activities. The results obtained from models offer important structural insight into designing novel CCR2B antagonists before their synthesis.

Chemokines or chemotactic cytokines are small molecular weight (6–15 KD) proteins that modulate inflammatory and immune responses through promotion of cell migration, cell adhesion, and transmigration.[1-7] They are divided into four classes dependent on the arrangement of conserved cysteines in the N-terminal region. Monocyte chemo attractant protein-1 (MCP-1) is a member of the CC-chemokine family and it selectively recruits leukocytes from the circulation to the site of inflammation through binding with the chemotactic cytokine receptor 2B (CCR2b).[89] This receptor is a member of the seven transmembrane receptor families (7-TM) and is expressed on the surface of monocytes and macrophages.[10] The recruitment and activation of selected populations of leukocytes is a key feature of a variety of inflammatory conditions. This response is crucial for host defense during inflammation, but the secretory products of white blood cells may increase injury by damaging surrounding healthy tissue.[11-14] These effects are mediated principally through activation of intracellular signaling pathways following binding of MCP-1 to the CCR2b. MCP-1 is a potent chemotactic and activating factor for monocytes and memory T-cells. It regulates adhesion molecule expression and cytokine production.[15] It also induces superoxide anion and lysosomal enzyme release from human monocytes.[16] The role of MCP-1 is seen in pathophysiology of a wide range of acute and chronic inflammatory diseases such as rheumatoid arthritis,[17-19] atherosclerosis,[20-23] asthma,[24-26] psoriasis,[2728] and transplant rejection.[29-32] This involvement is evidenced by studies showing elevated MCP-1 expression correlated with leukocyte infiltration in vivo,[33-35] the use of neutralizing antibodies,[3637] and through both animal receptor[38] and ligand[39] knockout studies.

An MCP-1 receptor antagonist thus represents an attractive target for drug discovery, and this has prompted an intense period of pharmaceutical research. Several companies have reported the discovery of potent small molecule antagonists of the CCR2b, showing varying degrees of selectivity over closely homologous receptors.[40-45]

Nowadays, the use of three-dimensional quantitative structure–activity relationship (3D-QSAR) techniques, such as comparative molecular field analysis (CoMFA) and comparative molecular similarity indices analysis (CoMSIA),[46-48] is routine in modern drug design to help understand drug–receptor interaction. It has been shown in the literature that these computational techniques can strongly support and help the design of novel, more potent inhibitors by revealing the mechanism of drug–receptor interaction.[49-51] In this study, we have developed predictive 3D-QSAR models using (R)-3-aminopyrrolidine derivatives as antagonists of the CCR2b.[10]

Materials and Methods

Data sets

We took the in vitro inhibitory activity data (IC50, nM) of a series of (R)-3-aminopyrrolidine derivatives, reported by Moree et al,[10] for the study. We have omitted 12 molecules whose IC50 values were not quantitatively reported from the given 62 antagonists.

We have considered the remaining 50 compounds for QSAR study. The IC50 values used in this study span a range of three orders of magnitude ranging from 3.2 to 4330 nM. Although major portion of compounds weighted toward the high-potency end of the spectrum, activities are acceptably distributed across the range of values. Thus, these molecules provided a broad and homogenous dataset for the 3D-QSAR study. The generation of reliable models is dependent on the creation of appropriate training and test sets. We have carefully chosen the training set molecules considering the following rules: composition of the QSAR training and test sets individually is necessarily representative of the whole data set in terms of structural complexity. Thus, these sets should show the appropriate chemical diversity and distribution of biological property across the range of IC50 values. The most active and the least active molecules should always be a part of training set, as these are the representative compounds for diversity in biological property. Thus, data set was divided into training (37 compounds) and test (13 compounds) sets [Table 1] in the ratio of 2.85:1.

Table 1

Structures and biological activities of the training and test set molecules

We have used the partial least squares (PLS) method for all QSAR analyses. We have converted CCR2b inhibitory activities into the corresponding PIC50 (-log IC50) values and used them as dependent variables, whereas CoMFA and CoMSIA descriptors are used as independent variables in the PLS regression analyses to derive QSAR models. The predictive ability of the models was assessed by their q2 values.

Molecular modeling

We performed the three-dimensional structure building and all modeling using the Sybyl program package, version 6.7 on a Silicon Graphics Fuel workstation.[52] CoMFA and CoMSIA studies require that the 3D structures of the molecules to be analyzed should be aligned according to a suitable conformational template, which is assumed a bioactive conformation.[53] As the bioactive conformations of these inhibitors are not known, the lowest energy conformations were reasonable initial structures to perform 3D-QSAR calculations. We performed energy minimizations using the Tripos force field[53] and the Gasteiger–Hückel[54] charge with a distance-dependent dielectric and Powell conjugate gradient algorithm. The criterion of convergence was 0.05 kcal/mol. Subsequently, the lowest energy conformation of each structure was submitted to optimization with the semi-empirical program MOPAC 6.0 and applying the AM1 Hamiltonian.[55]

Compound alignment

One of the most important adjustable parameters in 3D-QSAR is the relative alignment of all the molecules to one another so that they have a comparable conformation and a similar orientation in space.[56] Usually the bound conformation of a ligand from X-ray crystallographic or NMR studies makes a good starting point for the alignment. However, due to the unavailability of experimentally determined structures, we have opted to take the most potent inhibitor of the data set, compound 71[3] as a template for superimposition, assuming that its conformation represents the most bioactive conformation of the (R)-3-aminopyrrolidine derivatives. Each analog was aligned to the template by rotation and translation to minimize the RMSD between atoms in the template and the corresponding atoms in the analog using the DATABASE-ALIGN option in SYBYL. This method involves the alignment of molecules in a structurally and pharmacologically reasonable manner based on the assumption that each compound acts via a common macromolecular target-binding site. The most active compound is shown in Figure 1, and the aligned compounds are shown in Figure 2.

Figure 1

Most active molecule with labeled scaffold

Figure 2

Alignment of all molecules

CoMFA studies

To derive the CoMFA descriptor fields, a 3D cubic lattice with grid spacing of 2Å in x, y, and z directions was created to encompass the aligned molecules. In CoMFA, Lenard-Jones (6–12 interactions), the steric interaction field, and Columbic electrostatic potentials (1/r) were calculated at each lattice intersection. The grid box dimensions were determined automatically in such a way that region boundaries were extended beyond four Å in each direction from coordinates of each molecule. The Vander Waals potentials and Columbic terms, which represent steric and electrostatic fields, respectively, were calculated using Tripos force field. A sp3 hybridized carbon atom with radius 1.52 Å bearing +1 charge served as probe atom to calculate steric and electrostatic fields.

The CoMFA steric and electrostatic fields generated were scaled by the CoMFA standard option available in SYBYL. A PLS approach, an extension of multiple regression analysis, was used to derive 3D-QSAR, in which the CoMFA descriptors were used as independent variables and PIC50 values as dependent variables.[57-59] We performed cross-validation analysis for selecting the model that is most likely to have predictive values. The intensity of the cross-validation process is controlled by selecting the number of 10 groups. After the optimum number of components was determined, we performed a non-cross-validated analysis. We have computed the r2, PRESS (the root mean predictive error sum of squares), r2, F value, and standard error of estimate values according to the definition in the SYBYL. The cross-validated coefficient was calculated using the following equation

where γ, γ, and γ are predicted, actual, and mean values of the target property (PIC50), respectively. We have used the following formula to calculate lowest standard error of prediction

The non-cross-validated PLS analyses were performed with column filtering value of 2.0, to reduce analysis time with small effect on the q2 values. We further assessed the robustness and statistical confidence of the derived models, through bootstrapping analysis for 100 runs. We have examined the predictive power of 3D-QSAR models, derived by using the training set with an external test set of 13 molecules. The predictive ability of the models is expressed by the predictive r2 (r2) value, which is analogous to cross-validated r2 (r2) and is calculated using the following formula:

where SD (standard deviation) is the sum of the squared deviations between the biological activities of the test set and the mean activity of the training set molecules and PRESS is the sum of the squared deviations between predicted and actual activities for every compound in the test set. The activity of the test set was predicted by the CoMFA model using the predict Command. CoMFA coefficient maps were generated by interpolation of the pair wise products between the PLS coefficients and the standard deviations of the corresponding CoMFA descriptor values.

CoMSIA studies

CoMSIA was performed to evaluate steric, electrostatic, hydrophobic, hydrogen bond donor, and hydrogen bond acceptor properties of molecules by using the standard options in SYBYL. The steric, electrostatic, hydrophobic, H-bond donor, and H-bond acceptor fields were calculated separately using the sp3 carbon atom probe with a charge of 1 provided in SYBYL7.0. Similar to CoMFA, a data table has been constructed from similarity indices calculated at the intersections of regularly spaced lattice (2Å spacing). We have calculated similarity indices AF, K between the compounds of interest, and a probe atom according to the following equation:

where q is the grid point for molecule j; ω is the actual value of the physic chemical property k of atom i; ω indicates probe atom with charge 1, radius1Å, hydrophobicity 1, H-bond donor, and acceptor property 1; α is an attenuation factor; and r2 is the mutual distance between the probe atom and grid point q and atom i of the test molecule. The default value of α is 0.3.

Model validation

The predictive power of CoMFA and CoMSIA models was further validated by using an external test set (inhibitors marked with “d” in Table 1). The inhibitors in the test set were given exactly the same pretreatment as the inhibitors in the corresponding training set. The correlation between the experimental and predicted activity for models was calculated as r2 value. We have also performed a cross-validation that is based on Fischer randomization test method.

Results and Discussion

We have used CoMFA and CoMSIA techniques to derive 3D-QSAR models on novel series of (R)-3-aminopyrrolidine-based compounds acting as CCR2b antagonists. The biological activity of negative logarithm PIC50 was used as a dependent variable.

We have used the low-energy conformer obtained from the AM1 optimization as template and aligned all compounds using DATABASE ALIGNMENT method. We generated various 3D-QSAR models and selected the best one based on statistically significant parameters obtained. We obtained the final model with 37 and 13 molecules in the training and test sets, respectively. The predictive power of the 3D-QSAR models, derived using the training set, was assessed by predicting biological activities of the test set molecules. In 3D-QSAR studies r2 of 0.3 is considered statistically significant.[60] In view of it, we can consider the models having r2 > 0.5 much better and statistically significant. Final predicted versus experimental PIC50 values for both CoMFA and CoMSIA models and their residuals (for training and test set compounds) are given in Table 2.

Table 2

Predicted and residual activities of the training and test set molecules from CoMFA and CoMSIA analyses

CoMFA analysis

The CoMFA analysis with contribution of the steric and electrostatic fields [Table 2] yielded a cross-validated, r2 = 0.847 with five components, non-cross-validated r2 of 0.977, a conventional r2 (leave one out), r2 of 0.856, an F value 267.930, and a predictive r2, r2 of 0.673. The results of CoMFA study are given in Table 3. The graphs of actual versus predicted activities for the training and test sets of molecules are depicted in Figure 3. CoMFA contours were generated using this model. To further assess the robustness of the model, bootstrapping analysis (100 runs) was performed and an r2 of 0.988 (S.Dbs 0.005) was obtained, further establishing the strength of the model. Figure 4 shows the histogram of residual values obtained from CoMFA analysis. The steric and electrostatic contributions were found to be 54.6% and 45.4%, respectively. We have further used data set and alignment of CoMFA for CoMSIA analysis.

Table 3

Summary of CoMFA results

Figure 3

Graph of actual versus predicted activity of training and test set molecules from CoMFA analysis

Figure 4

Histogram of residual values obtained from CoMFA analysis

CoMSIA analysis

CoMSIA is similar to CoMFA but uses a Gaussian function rather than Columbic and Lennard–Jones potentials to assess the contribution from different fields. CoMSIA was performed using steric, electrostatic, hydrophobic, hydrogen bond donor, and hydrogen bond acceptor fields. 3D-QSAR models were generated using all the above fields, and the results of study are summarized in Table 4.

Table 4

Summary of CoMSIA results

The CoMSIA model yielded the cross-validated q2 = 0.715 with six components, non-cross-validated R2 of 0.964, F value of 135.666, bootstrapped R2 of 0.975, and a predictive r2 of 0.611. The steric, electrostatic, hydrophobic, donor, and acceptor field's contributions were 15.8%, 27.2%, 8.2%, 29.3%, and 19.5%, respectively. The graph of the actual versus predicted biological activities for the molecules is shown in Figure 5. Histogram of residual values obtained from CoMSIA analysis is depicted in Figure 6.

Figure 5

Graph of actual versus predicted activities from CoMSIA model

Figure 6

Histogram of residual values from CoMSIA analysis

The field values generated at each grid point were calculated as the scalar product of the associated QSAR coefficient and the standard deviation of all values in the corresponding column of the data table (STDDEV*COEFF), plotted as the percentage contributions to QSAR equation. The CoMFA and CoMSIA contour maps developed are shown in Figures 7 and 8, respectively.

Figure 7a

CoMFA STDEV COEFF contour maps: steric fields; green contours indicate regions where bulky groups increase activity, while yellow contours indicate regions where bulky groups decrease activity. Active compound 71 is displayed in the background for reference

Figure 8a

CoMSIA STDEV*COEFF contour maps: steric fields; green contours indicate regions where bulky groups increase activity, while yellow contours indicate regions where bulky groups decrease activity Active compound 71 is displayed in the background for reference

Internal validation

In model validation, we examined the internal predictive power of the models and their ability to reproduce biological activities of the compounds for the training set. The computed CCR2b inhibitory activity from the CoMFA and CoMSIA showed a good correlation with experimental inhibitory activity.

External validation

We have performed external validation of QSAR models to verify the excellent statistical parameters that were obtained and to investigate whether the activity of (R)-3-aminopyrrolidine derivatives from external data could be predicted well with this model. We have selected a set of 13 compounds, as a test set, from the 50 compounds for the validation experiments. The ultimate test for predictability of QSAR analysis in the drug design process is to predict the biological activity of compounds that have not been included in the training set. The r2 was calculated and we have obtained values of 0.673 and 0.611 for CoMFA and CoMSIA, respectively. Thus, the CoMFA model displays higher predictivity both in regular cross-validation and in the prediction of the test compounds.

Fischer statistics (F test)

Fischer statistics (F) is the ratio between explained and unexplained variance for a given number of degrees of freedom. F test is a variance-related statistics that compares two models differing by one or more variables to see if the more complex model is more reliable than the less complex one. The model is supposed to be good if the F test is above a threshold value, i.e., tabulated value. The larger the value of F, the greater is the probability that the QSAR equation is significant. The F values for the CoMFA and CoMSIA models were 267.930 and 135.666, respectively.

Visualization and analysis of contour maps

CoMFA analysis

Steric analysis: The CoMFA contour maps have permitted an understanding of the steric and electrostatic requirements that represent the QSAR result. Figure 7a shows the contour map derived from the CoMFA PLS model. The contours are mapped on to the most active compound 71. The contour plots may help to identify important regions where any change may affect the activity of molecule. Furthermore, they are helpful in identifying important features contributing to increased activity of molecules. The steric interactions are represented by green, and region where steric bulk decreases the activity is depicted in yellow color. The yellow contour map shown extending from the tri-fluoro carbon to the left lower region indicates that the longer chain substituents toward this spatial distribution decrease the activity. The CoMFA steric contour map was less informative than the electrostatic map as the green contours were missing.

Electrostatic analysis: The electrostatic contour maps obtained from CoMFA analyses were mapped on to the compound 71 and are shown in the Figure 7b. The presence of red contour map at the C13 position of the methyl benzamide scaffold indicates that the substituents having electronegative group attached at this position show improved activity as observed in the compounds 57–72 against the unsubstituted compounds 12–34. The compound 71 is showing highest activity as its electropositive substituent is mapped to the blue contour appeared at this position. Substitution of electronegative group in this position has detrimental effect on activity as observed in compounds 39, 40, and 45.

CoMSIA analysis

Steric contour analysis: Figure 8a shows sterically favored (green) and disfavored (yellow) regions. The yellow contour map shown toward the left lower region of this position indicates that the longer chain substituents toward this spatial distribution decrease the activity. The contour at this position is also observed in CoMFA map [Figure 7]. The yellow contour around the pyrrolidin ring signifies that the steric substituents are disfavored (in case of compound 34). Yellow blocks appeared at fourth and fifth positions of phenyl ring indicate that the bulky groups with longer chain substituents decrease the activity.

Electrostatic contour analysis: The electrostatic contour maps shown by the CoMSIA model [Figure 8b] are more informative than that of CoMFA model [Figure 7b], as they give more detailed picture in the case of substitutions. The blue contour map positioned near CF3 of benzene ring illustrates that electropositive groups are favored at this position.

Figure 8b

CoMSIA STDEV*COEFF contour maps: electrostatic fields; blue contours indicate regions where electropositive groups increase activity, while red contours indicate regions where electronegative groups increase activity. Active compound 71 is displayed in the background for reference

Figure 7b

CoMFA STDEV COEFF contour maps: electrostatic fields; blue contours indicate regions where electropositive groups increase activity, while red contours indicate regions where electronegative groups increase activity. Active compound 71 is displayed in the background for reference

As depicted from the Figure 8b, blue contour map of CoMSIA model shown at the methyl benzamide near the red contour map denotes that electropositive groups are favored at this region. At this position, the longer carbon linker chains show increased activity as noted in compounds 31, 21, and 61. The red contour map observed at methyl group near the blue contour specifies the region for electronegative groups. The electronegative group substituted compounds 58, 59, and 72 are observed to have comparatively higher activity. It was also observed in other compounds having phenyl ring with substituted halogens showing higher activity (56 and 57). The compounds 12, 17, and 15 with electropositive substituted phenyl ring at this position are showing decreased activity.

Hydrogen bond acceptor and donor contour analysis: Figure 8c and Figure 8d depict the hydrogen bond acceptor and donor contour maps of the CoMSIA models. Magenta contours encompass regions where a hydrogen bond acceptor will lead to improved biological activity, while an acceptor located near the red regions will result in impaired biological activity. In donor contour map, 6d, cyan color indicates the regions where hydrogen bond donor acts as favored and orange color refers to the disfavored regions. There are two magenta contour maps in Figure 8c, surrounding the two benzene rings, which supports the requirement of H-bond acceptor. This can be seen in compounds of table 3, 31, and 21. The small red contours in the same figure indicate that any hydrogen bond acceptors are not favored in these areas. This offers an explanation for the worse biological activity of compounds 45, 39, 40, and 15 as their substituents having acceptor group at this position are mapped to this region. The cyan contours in Figure 8d indicate that presence of H-bond donors increases the biological activity. The presence of strong H-bond donor in compounds with first scaffold is responsible for their high activity.

Figure 8c

CoMSIA STDEV*COEFF contour maps: Acceptor fields; the magenta contour for H-bond acceptor group increase activity, red indicates the disfavor region. Active compound 71 is displayed in the background for reference

Figure 8d

CoMSIA STDEV*COEFF contour maps: Donor fields; the cyan contour for H-bond donor favor region, purple indicates the disfavor region. Active compound 71 is displayed in the background for reference

Hydrophobic contour analysis: Yellow and white contours enclose regions favorable for hydrophobic and hydrophilic groups, respectively. The white contours in Figure 8e support the importance of hydrophilic substitutions. This is responsible for the higher biological activities of molecules (64, 68, 69, and 70).This hydrophilic interaction might be very important for binding affinity, since this feature was also observed in CoMFA and CoMSIA steric contour maps. A small yellow contour mapped very near to methyl group of benzene ring illustrates that hydrophobic group at this position is important for increased biological activity.

Figure 8e

CoMSIA STDEV*COEFF contour maps: Hydrophobic fields; the yellow contour for hydrophobic favor region, white indicates the hydrophilic favored region. Active compound 71 is displayed in the background for reference

Conclusions

In this study, we have investigated the CoMFA and CoMSIA models based on a training set of 37 structurally diverse (R)-3-aminopyrrolidine series, followed by validation of the results by an external test set of 13 analogues. Despite the lack of structural information on ccr0 2b, the design of potent inhibitors can be attempted by means of CoMFA and comparative molecular similarity indices, well-established 3D-QSAR techniques. These models demonstrated excellent internal and external predictive ability, which was shown by several strategies including cross-validation, predictive r2, and test set predictions.

Overall, the CoMFA model gave good results. The CoMSIA model was more valuable for the three fields that contributed significantly (hydrophobic, hydrogen bond acceptor, and hydrogen bond donor). The CoMSIA analysis indicated that variations in the activity are dominated by hydrophobic interactions. The excellent correlation with several experimental studies suggests that these 3D-QSAR models are reliable, helping us to understand the binding interaction of these compounds and providing a helpful guideline for further lead optimization. The features derived from the above models bear a close correlation with the structural variations inherent in the training set, so other structurally distinct data results in diverse features causing different conclusions. In summary, our preliminary findings may aid in identifying potent and specific (R)-3-aminopyrrolidine series that may be used as potent antagonists of the CCR2b and offer more significant insights into the overall pharmacology of this system.