key: cord-0862089-s0cjlflq authors: Chen, Haifeng; Li, Qiang; Yao, Xiaojun; Fan, BoTao; Yuan, Shengang; Panaye, A.; Doucet, J. P. title: CoMFA/CoMSIA/HQSAR and Docking Study of the Binding Mode of Selective Cyclooxygenase (COX‐2) Inhibitors date: 2004-02-18 journal: QSAR Comb Sci DOI: 10.1002/qsar.200330844 sha: eb35a443c49335c3ee31963c7d79e0676bc01c4f doc_id: 862089 cord_uid: s0cjlflq The intermolecular interaction between four types of anti‐inflammatory inhibitors (oxazoles, pyrazoles, pyrroles and imidazoles) and COX‐2 receptor was studied. The results of docking suggest that they have similar interaction mechanism. The most active compounds of these four types of inhibitors could both form several hydrogen bonds with residues His90, Arg513, Leu352 and Arg120, and develop hydrophobic interaction with residues Phe518, Leu352 and Leu359. This is consistent with the investigation reported by R. G. Kurumbail et al. (Nature. 1996, 384, 644‐648). A common 3D‐QSAR model could be constructed with these four categories of COX‐2 inhibitors using the method of docking‐ guided conformer selection. The cross‐validated q(2) values are found as 0.741 and 0.632 for CoMFA and CoMSIA respectively. And the non‐cross‐validated r(2) values are 0.887 and 0.885. 54 inhibitors constitute the test set used to validate the model. The results show that this model possesses good predictive ability for diverse COX‐2 inhibitors. Furthermore, a HQSAR model was used to evaluate the influence of substituents on anti‐inflammatory activity. Compared with the results of previous works, our model possesses significantly better prediction ability. It could help us to well understand the interaction mechanism between inhibitors and COX‐2 receptor, and to make quantitative prediction of their inhibitory activities. In human body inflammation appears when immune system is attacked by virus (for example coronavirus, influenza ...). It is necessary to investigate the anti-inflammatory mechanism for the goal of designing new drugs. COX-1 and COX-2 are two cyclooxygenase (COX) isoforms [1] . Their principal pharmacological effect is inhibiting prostaglandin synthesis. This discovery led to the hypothesis that side effects such as ulcers and renal failure associated with the clinically used nonsteroidal anti-inflammatory drugs (NSAIDs) are caused by the inhibition of COX-1. Whereas COX-2 is an inducible enzyme, it is mainly produced during inflammation processes [2] . Then the study of selective inhibition of COX-2 led to a new class of anti-inflammatory, analgesic and antipyretic drugs with significantly reduced side effects. Recent works suggest that inhibiting COX-2 could also be an important anti-cancer strategy [3] and could be used to delay or slow down the clinical expression of Alzheimer×s disease [4] . But extended use of NSAIDs, such as ibuprofen, may increase the risk of developing Alzheimer disease up to 80% [5] . Therefore we need to develop more specific and efficient COX-2 inhibitors with improved safety profile [6] . A large number of research studies aimed at finding selective COX-2 inhibitors have been reported [7 ± 10] . Many studies have been carried out using computer simulations to develop protocols and methods for designing new COX-2 inhibitors such as oxazoles, pyrazoles, pyrroles and imidazoles [11 ± 15] . But these investigations did not clearly explain how substituents in these inhibitors influence the anti-inflammatory activity. In this paper, we report the binding mode between four types of selective inhibitors mentioned above and the COX-2 receptor using automated molecular docking methodology. The study will show that a unique docking model can be built up on a population of 227 diarylheterocyclic derivatives, suggesting a possible common inhibitory mechanism for the four types of above mentioned compounds. This allows us for attempting to construct 3D-QSAR models by using the conformers obtained from molecular docking investigation. These 3D-QSAR models can help us to better understand the interaction mechanism between inhibitors and COX-2 receptor, and make quantitative prediction of their inhibitory activities. The structures and activities IC 50 (IC 50 is the concentration in mM for 50% inhibition of COX-2 enzyme) for four types of COX-2 inhibitors (oxazoles, pyrazoles, pyrroles and imidazoles), extracted from literature [16 ± 20] , are gathered in tables 1 to 7. The structures marked with ™*∫ symbol constitutes the test set, and the others the training set. The complex of cyclooxygenase (COX-2) receptor and SC558 was extracted from Brookhaven Protein Databank (PDB code: 1CX2). The structure of ligand SC558 is shown in figure 1 . The most active structures for each type (A9, B2, C13 and D59) were optimized using Gaussian98 (at HF/6-31G level) [21] . Other structures were built based on these geometries, and then optimized using Tripos force field [22] with SYBYL 6.9 [23] . All calculations were performed on a SGI Octane2 workstation. Because the four types of COX-2 inhibitors investigated have a skeleton similar to that of SC558 (Figure 1 ), we only need to align their five-membered heterocycle with that of SC558. The active site was determined based on the cocrystallized SC558-protein complex. Firstly, the ligand SC558 was extracted from the cyclooxygenase complex, and then the remaining COX-2 receptor structure was completed by adding the missing polar hydrogen atoms and residues. This structure was optimized with constrained CoMFA/CoMSIA/HQSAR and Docking Study of the Binding Mode of Selective Cyclooxygenase (COX-2) ... molecular dynamics. Partial atom charges of COX-2 receptor were calculated with Kollman-all-atom approximation [24] . Gasteiger-H¸ckel charges were calculated for all inhibitor molecules, as recommended in the AutoDock 3.0 package which was used for performing automated docking of inhibitors to COX-2 receptor. The development and the principle of this software have been described elsewhere [25 ± 26] . During the docking process, a series of the docking parameters were set on. The number of generations, energy evaluations, and docking runs were set to 370 000, 1 500 000, and 50, respectively. CoMFA results may be extremely sensitive to a number of factors, such as conformers, alignment rules, overall orientation of aligned compounds, lattice shifting, step size and the probe atom type [27] . Where S atom is the total number of atoms over which the RMSD is measured and d i the distance between the coordinates of atom i in the two conformers which are overlaid. The conformation with minimum RMSD value is selected in CoMFA analysis. After consistently aligning the molecules within a lattice of 2 ä mesh which extend 4 ä units beyond the aligned molecules in all directions; a sp 3 carbon atom with 1 net charge was employed as probe. The steric and electrostatic interactions between the probe and the atoms of the molecule were calculated. Electrostatic interactions are modeled using a Coulomb potential and Van der Waals interactions by Lennard-Jones potential. The regression analysis was carried out using the partial least-squares (PLS) [28] method. The final model was developed with the optimum number of components (that yielding the highest q cv 2 ). The total set of inhibitors was divided into two groups in an approximate ratio 3 : 1 (for example, 173 in the training set and 54 in the test set). The selection of test set and training set compounds was done manually such that low, moderate and high activity compounds are present in roughly equal proportions in both sets. CoMSIA similarity indices [29 ± 30] were derived within the same lattice box of the CoMFA calculations. Five physicochemical properties related to steric, electrostatic, hydrophobic, hydrogen bond donor, and hydrogen bond acceptor fields were evaluated using the probe atom. A Gaussiantype distance dependence was employed to describe the attenuation of the fields with distance, using the default value (0.3) as attenuation factor. This leads to a much smoother sampling of the fields around the molecules compared with CoMFA. Cross-validated q 2 usually serves as a quantitative measurement of the predictive power of CoMFA and CoMSIA. However, Cho et al. reported that q 2 value was sensitive to the orientation of aligned molecules on the computer terminal and might vary with the orientation by as much as 0.5 q 2 units [27] . So all molecules under investigation were geometrically aligned, then the molecular aggregate was rotated systemati- cally every 308 for X, Y, Z axis and translated every 0.2 ä. At every place, the cross-validated q 2 of PLS was evaluated, until the maximum q 2 value was found [31] . Molecular fingerprinting (hologram) is a method that represents a compound as a unique string of numbers or ™bins∫ [32] . The bins represent all of the unique fragments included within a particular molecule and are assigned by a cyclic redundancy check (CRC) algorithm [33] . A number of holograms with differing bin lengths are created for each molecule. The hologram length that leads to the best PLS analysis is used to construct HQSAR model. A set of factors such as fragment size, number of fragments, atom types, bond types, atom hybridization, hydrogen bond donor and acceptor, were modified to search the best cross-validation r 2 of the model. For each molecule, 50 conformers were selected to dock with COX-2 receptor. The most active molecules for four types of inhibitors (oxazoles, pyrazoles, pyrroles and imidazoles) were selected to study the binding mode with COX-2 receptor. The complexes formed by A9 of oxazoles, B2 of pyrazoles, C13 of pyrroles and D59 of imidazoles with COX-2 receptor were then investigated, as shown in Figure 2 . For all complexes, two oxygen atoms of sulphonamide are linked by hydrogen bond to residues His90 and Arg513. This is consistent with the investigation reported by R. G. Kurumbail et al. [6] . Beside these hydrogen bond interactions, there are also hydrophobic effects around sulphonamide group bonded to phenyl with the residues Phe518, Leu352, Leu359, Trp387 and Met522. This hydrophobic interaction seems to be favourable to the activity. The other residues such as Tyr355, Tyr385 and Ser530 are polar. The electrostatic interaction between these residues and phenyl ring could stabilise the formed complex. Except the complex of C13, there is one hydrogen bond between the nitrogen atom of sulphonamide and carbonyl oxygen of Leu352. This may explain that the sulphonamide containing compounds are more active than those containing sulphomethyl. In the complexes of A9 and B2, there is also electrostatic interaction between halogen atom and residue Arg120. It is perhaps why the anti-inflammatory activity of A9 and B2 is higher than that of D59. For B2-COX complex, there is a hydrogen bond between nitrogen atom of pyrazole ring and À OH of residue Tyr355 with 127.88 for N ¥¥¥ H À O angle and 2.85 ä for N ¥¥¥ H distance. These values (angle < 1408 and distance > 2.0 ä) suggest that this hydrogen bond should be very weak and does not significantly contribute to stabilise the complex. By contrast, for A9 complex, electrostatic interaction between the chlorine atom located in phenyl para position and the residue Trp387 could increase the interaction between ligand and receptor. This leads to the higher activity for A9 compared with B2. The docking study is in agreement with experiments. It seems that the antiinflammatory activity increases with the number of hydrogen bonds for these four compounds, as shown in the following sequences: Number of H-bonds: A9 (5) B2 (6) D59 (4) C13 (3) Activity: 9.00 8.77 8.52 7.70 The docking results suggest that the oxazoles, pyrazoles, pyrroles and imidazoles derivatives here investigated, might have similar interaction mechanism with COX-2. There are at least three hydrogen bonds, favourable hydrophobic interactions between sulphone bonded to phenyl and residues Phe518, Leu352, Leu359, Trp387 and Met522, and electrostatic interactions with residues Tyr355, Tyr385 and Ser530. Figure 3 illustrates the superposition results of A9, B2, C13, and D59 with ligand SC558. Their RMSD values with SC558 are 1.38, 1.32, 1.36 and 1.39 ä, respectively. Their benzenesulphonyl and five-membered heterocycle could be aligned together. We could then attempt to construct a common 3D-QSAR model for these four types of COX-2 inhibitors. At the same time, we investigated these compounds using the method of hierarchical clustering. 34 descriptors were calculated with default methods implemented in SYBYL 6.9, including 27 Charged Partial Surface Area (CPSA) descriptors [34] , 4 polarity, volume, surface and hydrophobicity descriptors. The aim of this analysis was to examine whether the 227 compounds can be gathered into a unique model or must be split into different sub populations according to their structural features. As shown in Figure 4 , cluster analysis indicates that the 227 compounds were classified into two categories. As an example, the black dots illustrate one branch of the cluster, including compounds B1, B32, B35, D18, D61, D97, A9, A10 and A11. This mixture of different types of compounds (A, B and D) suggests that these four categories of molecules cannot be really separated according to their structure features. Therefore all these compounds could be treated together in a common 3D-QSAR model. According to the minimum RMSD of the training set, the most probable conformer of every compound was selected. The alignment diagram of the 173 compounds of the training set within the active site of COX-2 is shown in Figure 5 . For each type of COX-2 inhibitors, the structures were divided into training and test sets. Then, the structures of these training sets were put together for building a new model. Based on this model, two other models were constructed by the method of all-orientation and all-placement searches and random orientation. For all-orientation and all-placement searches, training set was rotated systematically every 308 for X, Y, Z axis and translated every 0.2 ä, then cross-validated q 2 of PLS was evaluated. Because the number of training set of pyrrole is not enough to build a 3D-QSAR model, only six models were constructed. Model 1 is for oxazole, model 2 for pyrazole, model 3 for imidazole, model 4 corresponds to the mixture of all training sets described above. Model 5 results from all-placement searching based on model 4. CoMFA/CoMSIA/HQSAR and Docking Study of the Binding Mode of Selective Cyclooxygenase (COX-2) ... We examine now the three correlations between experimental and predicted anti-inflammatory activities for the three classes of inhibitors (oxazoles, pyrazoles and imidazoles). correlations between EA and PA are shown in Figure 9 . This model could predict the pIC 50 values accurately for most compounds of test set, including D96 and D102, except compounds A33 and B33, which have relatively large error values (1.13 and 1.02 log units, respectively). Model 1 can not predict the pIC50 value accurately for compound A33 (error 1.02 log units), and model 2 is also incapable of accurately predicting the activity for molecule B33 (error 0.79). Analyzing the conformers used in the models, it is found that the phenyl ring of these two compounds is rotated about 808 with respect to the five-membered ring, whereas for other compounds, the phenyl group and five-membered ring are located on a same plane. This may be the major cause of the failure of the models. Comparing these four models, the sequences of the corresponding correlation coefficients r 2 for test sets are 0. Model 6 is the result of random orientation of the 173 compounds of the training set. Random orientation of training set was chosen by random method. The crossvalidated q 2 values of CoMFA and CoMSIA models of training set are 0.693 and 0.651, respectively. Model 5 is obtained by the method of all-orientation and all-placement searching. The cross-validated q 2 values of CoMFA and CoMSIA models of training set for this model are 0.748 and 0.638, respectively. The parameter of CoMFA for model 5 is much higher than that of model 6. It suggests that CoMFA model is sensitive to the orientation of investigated training set, while the orientation of aligned compounds has little influence to CoMSIA model [30] . The parameters of model 5 for CoMFA and CoMSIA approaches are similar to those of model 4. This shows that the method of docking-guided conformer selection could eliminate the influence of orientation for CoMFA model. This docking-guided conformer selection seems to be a valuable method to search for the active conformers and build 3D-QSAR model. Hydrogen bond interaction between NH 2 and residue Leu352 could stabilize the complex of ligand and COX-2. In general, the compounds with a group sulphonylamide at position 4 of the phenyl ring have higher activity than that of molecules bearing a methylsulphonyl at that position. The favorable negative charge (red contour) is near para substituent of benzene in position 5; green contour is near meta substituent in the same position. Because the para substituent of D2, D1 is F and Cl, respectively, their activity is larger than that of D3-D9 with para electropositive substituents ( À H, À Me, À OMe, À NHMe, À NMe 2 , À SMe and À SO 2 Me). The PLS analysis parameters and predicted log(1/IC 50 ) are given in tables 8, 9 and 10. The steric, electrostatic, hydrophobic, hydrogen bond donor and hydrogen bond acceptor field contributions are 0.062, 0.304, 0.227, 0.259 and 0.149, respectively. We remark that the electrostatic, the hydrophobic and the hydrogen donor contributions are about 80%. This shows that these three fields are the main factors affecting anti-inflammatory activity in CoMSIA models. Figure 11 (A and B) are the contour plots of CoMSIA model. The meanings of the different contours for steric and electrostatic fields are the same as in CoMFA model. In figure 11 -B, yellow contours indicate the regions where hydrophobic groups enhance the activity, and the white contours show the regions where hydrophilic groups increase the activity. In cyan regions hydrogen bond donor groups are favourable for activity, whereas purple contours delineate regions where they are unfavourable. Magenta contours indicate the zones where hydrogen bond acceptor group could increase the activity, and in red coloured regions, hydrogen bond acceptor groups would decrease the activity. The substituents at position 2 of the five-membered heterocycle near yellow-coloured, red-coloured and redcoloured (hydrogen bond acceptor field) regions, indicate that electronegative and small bulk groups are favourable to activity, while hydrogen bond acceptor groups will decrease it. This could explain that A7, A15, A21, A37, A39, A40, A44, A46, A51, B16, B19, C11, C17, C20, D75, D82 and D102 with À CH 2 OH, À CH 2 COOH, À SH, À OH, À NH 2 , À CONH 2 , À CHO substituents at position 2 of the fivemembered cycle are less active than molecules bearing À CF 3 groups. In our docking investigation, we find that these former substituents with polar hydrogen at position 2 could form hydrogen bond with the residue Arg120 of COX-2, while the sum number of hydrogen bond between ligand and receptor is less than the compound with substituents À CF 3 . Therefore their binding free energies are lower than the latter one. Previous literature did not report why these compounds have low activity [11 ± 12] . Electrostatic (redcolored) and hydrogen bond acceptor fields near the substituent of five-membered heterocycle at position 3 show that electronegative groups could increase the activity and hydrogen bond acceptor groups decrease it. Compounds of B1 ( À Cl), B2 ( À F) and B17 ( À Cl) have higher activity than B6 ( À OH) and B14 ( À NH 2 ). This agrees with the indication of CoMSIA model. There are red-colored, cyan-colored and magenta-colored regions near the para substituent of the six-membered aromatic ring in position 4. It shows that polar, hydrogen bond donor and acceptor groups are favorable to the activity. This suggests that these molecules with substituents À SO 2 NH 2 on position 4 are more active than those with À SO 2 CH 3 . The result is in agreement with our docking investigation discussed above. Five-membered heterocycle is buried in a white-coloured contour, showing that hydrophobic interaction between residues and this ring is unfavourable to activity. The para substituent of phenyl at position 5 is near magenta-coloured, green-coloured and yellow-coloured (hydrophobic field) contours, suggesting that hydrogen bond acceptor, hydrophobic and bulky groups could increase anti-inflammatory activity. This could explain that compounds A2(4-Cl), A3(4-Br), A9(4-Cl), A10(3,4-Cl 2 ), A29(4-F), B40(4-SMe), B44(4-NMe 2 ), D10(4-Cl) and D11(4-F) are more active than A19, A23, A26, A28, B5-B15, molecules lacking substituent at the same position. There are blue-coloured, green-coloured and yellow-coloured (hydrophobic field), red-coloured (hydrogen bond acceptor field) areas near the meta substituent of phenyl at position 5. This shows that electropositive, bulky and hydrophobic groups are favourable to activity while groups with hydrogen bond acceptor are unfavourable. The CoMSIA model indicates that bulky substituents at ortho position of 5-phenyl are unfavourable to the anti-inflammatory effect. The activity of compounds have the order 21(2-F)