3D QSAR AND DOCKING STUDY OF INDOLE DERIVATIVES AS SELECTIVE COX-2 INHIBITORS

Objective: Non-steroidal anti-inflammatory agents (NSAIDs) continue to be one of the most widely used groups of therapeutic agents. QSAR (quantitative structure-activity relationship) approach is a very useful and widespread technique for drug design. 3D QSAR facilitates evaluation of three-dimensional molecular fields around molecules and generates a relationship of these fields' values with the activity. 
Methods: 3D QSAR study was performed on selected twenty-four compounds from synthesized indole derivatives using the stepwise variable selection k-nearest neighbor (kNN) molecular field analysis approach for indicating the contribution of the steric and electronic field for activity. The docking study was performed to further confirm the binding affinity of synthesized molecules (ligands) to COX-2 enzyme as well as to study binding nature. 
Results: Statistically significant model was generated using VLife Molecular Design Suite 3.5 software with cross-validated correlation coefficient q2 of 0.9461 and high predictive correlation coefficient (Pred_r2) of 0.8782 indicating that the model is robust. The results of docking study suggest that the synthesized compounds have a comparable binding affinity with the COX-2 enzyme. 
Conclusion: The present study may prove to be helpful in the development and optimization of existing indole derivatives as anti-inflammatory agents with selective COX-2 inhibition.


INTRODUCTION
Non-steroidal anti-inflammatory drugs (NSAIDs) inhibit cyclooxygenase (COX), the enzyme responsible for the conversion of arachidonic acid to prostaglandins. COX exists in 2 isoforms. COX-1 is a ubiquitous constitutive isozyme producing prostaglandins responsible for homeostatic functions such as maintenance of gastrointestinal (GI) mucosal integrity. COX-2 is largely a cytokine-induced isozyme producing prostaglandins that mediate pain and inflammation [1]. NSAIDs inhibit both COX-1 and COX-2 to varying degrees. Thus, the therapeutic effects of conventional NSAIDs are derived from inhibition of COX-2, while the adverse effects of these agents, particularly in the upper GI tract, arise from inhibition of COX-1 activity [2][3]. Much recent effort thus has been made to produce selective inhibitors of COX-2 in the belief that these will lack the gastrointestinal damaging effects of traditional NSAIDs [4][5][6]. Diarylheterocycle class of compounds has been investigated extensively as COX-2 inhibitors. Literature survey revealed that indole derivatives, pyrazoline derivatives, and pyrimidine derivatives independently possess good anti-inflammatory, analgesic activity and selective COX-2 inhibitory effects [7][8][9][10][11][12]. Hence we tried to obtain greater selectivity for COX-2 enzymes with the use of the indole nucleus and other structural features of different COX-2 inhibitors in the designed molecules. Thus the concept of chemical hybridization of the indole nucleus with pyrazoline and pyrimidine was attempted.
The computer-aided prediction of biological activity in relation to the chemical structure of a compound is now a commonly used technique in drug discovery. Computational chemistry represents molecular structures as a numerical model and simulates their behavior with the equations of quantum and classical physics. Available programs enable scientists to easily generate and present molecular data including geometries, associated properties (electronic, spectroscopic and bulk). The usual paradigm for displaying and manipulating these data is a table in which compounds are defined by individual rows and molecular properties (or descriptors) are defined by the associated columns [14][15][16][17].
Three-dimensional quantitative structure-activity relationship (3D QSAR) facilitates evaluation of three-dimensional molecular fields around molecules and generates a relationship of these fields' values with the activity. The k-nearest neighbor (kNN) method is one of the methods for generating a relationship between activity and molecular field and provides an interpretation of results thus providing clues for designing new molecules. kNN-MFA (Molecular Field Analysis) requires suitable alignment of a set of molecules. This is followed by generation of a common rectangular grid around the molecules. The steric and electrostatic energies are computed at the lattice points of the grid using methyl probe of charge+1. These interaction energy values at the grid points are considered for relationship generation using kNN method.
The kNN technique is a conceptually simple approach to pattern recognition problems. In this method, an unknown pattern is classified according to the majority of the class memberships of its k nearest neighbors in the training set. The nearness is measured by an appropriate distance metric (e. g. molecular similarity measure calculated using field interactions of molecular structures). The standard kNN method is implemented as follows: (i) calculate distances between an unknown object (u) and all the objects in the training set (ii) select k objects from the training set most similar to object (u), according to the calculated distances (iii) classify object (u) with the group to which a majority of the k objects belong. An optimal k value is selected by the optimization through the classification of a test set of samples or by the leave-one-out crossvalidation. The variables and optimal k values are chosen using different variable selection methods like Simulated Annealing variable selection method.
Drug exerts its biological activity by binding to the pocket of receptor molecule usually protein. In their binding conformations, the molecules exhibit geometric and chemical complementarily, both of which are essential for successful drug activity. The computational process of searching for a ligand that is able to fit both geometrically and energetically into the binding site of a protein is called molecular docking.
Molecular docking helps in studying drug/ligand or receptor/protein interactions by identifying the suitable active sites in the protein, obtaining the best geometry of ligand-receptor complex and calculating the energy of interaction for different ligands to design more effective ligands. The target or receptor is either experimentally known or theoretically generated through knowledge-based protein modeling or homology modeling. The molecular docking tool has been developed to obtain a preferred geometry of interaction of ligand-receptor complexes having minimum interaction energy based on different scoring functions viz. only electrostatics, the sum of steric and electrostatic and dock score. This utility allows one to screen a set of compounds for lead optimization. Genetic algorithm (GA), Piecewise Linear Pairwise Potential (PLP) and Grid algorithms are used to minimize the interaction energy between ligand-receptor.
In the present research work, a data set of twenty-four molecules showing comparable COX 2 inhibitory activity was subjected to 3D QSAR analyses, in search of newer and potent anti-inflammatory agents with selective COX-2 inhibition. Statistically, significant models were generated, and the most robust models for 3D QSAR was obtained using stepwise variable selection kNN-MFA approach using V-Life Molecular Design Suite software version 3.5.

Optimization of molecules' structures
A data set of twenty-four molecules showing comparable COX 2 inhibitory activity measured was chosen for the present 3D QSAR study as shown in table 1. The biological activity was expressed as IC50 values measured on COX-2 enzyme.
Various approaches for the development of COX-2 inhibitors have appeared in the literature over the past several years. The methods employed include 3D QSAR studies such as comparative molecular field analysis (CoMFA), receptor surface analysis (RSA), kNN method.
All models reported are of sufficiently high quality to be used as an estimator of the biological activity of unknown compounds; however, each model has some limitations. We attempted to obtain further insights into the structural requirement of a diarylindole class of COX-2 inhibitors by the application of 3D-QSAR using kNN-MFA method. Molecular alignment was performed which is useful for studying shape variation with respect to the base structure selected for alignment. Further, the aligned molecules were used for 3D QSAR studies. The method used for alignment was Template based. A template structure was defined and used as a basis for alignment of a set of molecules. The geometries of the aligned molecules were stored in the Align Molecules subfolder in the Align folder. All the aligned molecules were opened from Align Molecules subfolder in the Align folder. The 24 aligned molecules are shown in fig. 1.

Calculation of descriptors for 3D QSAR
The 3D QSAR worksheet was launched at the beginning. The subfolder containing all aligned molecules was opened. The data of biological activity (IC50) was inserted in the column. This was followed by the calculation of the molecular local shape field descriptors for finding a relationship with the activity and various parameters for field calculation selected are shown below.

Docking study
Molecular docking techniques consist in finding the low-energy binding modes of a ligand within the active site of a macromolecule and evaluating the binding energy with a score [18]. Explanation of the selectivity of small sets of ligands has been attempted with accurate but time-consuming techniques [19][20][21]. Otherwise, automated docking methods may be used to estimate the COX-2 affinity for large molecular databases. The objective of this work had been to use structural information of the target to further confirm the binding affinity of synthesized molecules to COX-2 enzyme as well as study binding nature.

Optimization of protein
Docking studies were carried out in Vlife molecular docking suite 3.5 by using Biopredicta. All the molecular modeling studies were performed on Pentium Core2Duo workstation using Sybyl. Docking studies were carried out using COX-2 (PROSTAGLANDIN SYNTHASE-2; EC: 1.14.99.1) complexed with a selective inhibitor, SC-558 (1-phenylsulphonamide-3-trifluoromethyl-5bromophenylpyrazole). The 1CX2 isoform of human cyclooxygenase was downloaded from Protein Data Bank website. The tetramer is converted into the monomer. Water molecules, cofactors, and heme were deleted from protein. The SC558 reference molecule was extracted. The bond order, bond angle, peptide bond check performed using the local geometry check option from the analyze menu of the BioPredicta module. Hydrogens were added in molecule and energy minimized using Merck Molecular Force Field (MMFF).

Optimization of ligand
The structures of SC558 and all the compounds used in COX-2 enzyme in vitro assay were energy minimized using MMFF until the root mean square gradient values became smaller than 0.0001 kcal/mol A °. Cavity 2 in 1CX2 (protein and co-crystallized ligand) is the cavity where the co-crystallized ligand is located. This can be ascertained by using the Edit BioPolymer tool. Grid Docking was started. The start and progress of the grid docking process at different grid points in the cavity (docking site is chosen) was observed in the Output Window. On completion of the docking, the final minimum score (dock score interaction/docking energy of receptor-ligand) for the best ligand pose is displayed in the output window together with the location path of the docked receptor-ligand complex. Grid-based docking performed in cavity number 2 with a number of bump 4 and rotation angle 100 °. Docking scores were recorded. After merging and energy minimization the amino acid interactions were analyzed.  a Values are means of two determinations acquired using an ovine COX-1/COX-2 assay kit and the deviation from the mean is<10% of the mean value, b In vitro COX-2 selectivity index (COX-1 IC50/COX-2 IC50).

RESULTS AND DISCUSSION
The results of in vitro COX inhibition assay on the most active twenty-four molecules from all the series shows that these compounds are more selective towards COX-2 than COX-1 as shown in

Interpretation of 3D QSAR model
In the kNN-MFA method, several models were generated for the selected members of the training and test set and the corresponding best model is reported in table 3. In this model, steric and electrostatic fields' contribution was found to be 75: 25 % respectively. kNN results can be viewed as the stereo view of the molecular fields (descriptors) in a rectangular grid around the molecular unit as shown around aligned molecules.

Predictivity of kNN-MFA model
Predictivity of the model was evaluated by predicting the activity of the molecules belonging to the training set (internal predictivity) as well as molecules in the test set (external predictivity) as shown in tables 4 and 5. Graphical representation of the predictive power of the model is shown in fig. 5.    a) Negative range indicates that negative electrostatic potential is favorable for an increase in the activity and hence a more electronegative substituent group is preferred in the region. b) Positive range indicates that positive electrostatic potential is favorable for an increase in the activity and hence a less electronegative substituent group is preferred in that region.
2) Steric field: a) Negative range indicates that negative steric potential is favorable for an increase in the activity and hence a less bulky substituent group is preferred in the region.
b) Positive range indicates that positive steric potential is favorable for an increase in the activity and hence a more bulky substituent group is preferred in that region.
A green sphere (S 832) near the sulfonamide moiety and aryl group at 1 position of indole ring, suggests a more bulky group is required in this region as indicated by more positive values in the paranthesis. Another green sphere (S 720) at the para position of aryl ring at 2-position of indole ring, suggests less bulky groups are favored for activity in this region. Green sphere (S 1207) at the para position on aryl ring at 3-position of indole ring suggests that addition of bulk here may decrease the activity as a value in the parentheses is more towards the negative side.
Only one blue sphere (E 1436) is observed in the model near the para position of aryl ring at position 1 of the indole nucleus. Negative values in the parentheses suggest the more electronegative group is favorable in this position.
Molecules were optimized using MMFF method till gradient convergence 0.05 kcal/mol. For viewing the results of docking, compound IIIB17 has been selected according to the best-scored conformation predicted by the scoring function. The interaction of this ligand has been shown in fig. 7. The reference ligand SC558 docked in cavity 2 of a 1CX2 isoform of the cyclooxygenase enzyme is shown in fig. 8.   The hydrogen bonding distance between NH group of IIIB17 with NH of Lysine 97A was found to be 1.599 A ° (N-----H).
The pyrimidine ring at 3 positions was surrounded by active site amino acid residues His 356A, Gly 350A, His 351A, and Asn 581A. The indole ring found to interact with residues His 356A, Lysine97A.
The benzene ring on the pyrimidine nucleus was found to interact with Asn581A and Gly 192A.
Twenty four molecules from six series viz IA, IB, IIA, IIB, IIIA, and IIIB were studied by docking with an aim to reveal their binding interactions with COX-2 enzyme. Compounds IIIB17, IIIB11 and IIB2 shown good docking score comparable to SC558, confirms the results obtained through COX inhibition assays.

CONCLUSION
The present 3D QSAR study using kNN-MFA method indicated the contribution of steric and electronic fields in the activity of compounds. The model suggests that a more bulky group is required near the sulfonamide moiety/aryl group at 1 position of indole ring, as indicated by positive values in the parenthesis and less bulky groups are favored for activity at the para position of aryl ring at 2position of indole ring while addition of bulk at the para position on aryl ring at 3-position of indole ring may decrease the activity. The model further suggests that the more electronegative group is favorable at the para position of aryl ring at position 1 of the indole nucleus. The results of docking studies indicate that the compound IIIB17 with highest docking score was found to be involved in the hydrogen bonding with a residue Lys 97. The hydrogen bonding distance between NH group of IIIB17 with NH of lysine97A was found to be 1.599 A ° (N-----H). Thus the results obtained can be used for further modification and optimization of the indole derivatives for better anti-inflammatory activity with selective COX-2 inhibition.