COMPUTER AIDED DESIGN OF SELECTIVE CALCIUM CHANNEL BLOCKERS: USING PHARMACOPHORE – BASED AND DOCKING SIMULATIONS

COMPUTER AIDED DESIGN OF SELECTIVE CALCIUM CHANNEL BLOCKERS: USING PHARMACOPHORE - BASED AND DOCKING SIMULATIONS

Reetu* and Vipin Kumar

Institute of Pharmaceutical Sciences, Kurukshetra University, Kurukshetra-136119, Haryana, India

ABSTRACT:In the present study, 3-D QSAR analysis was performed on the previously synthesized and evaluated derivatives of novel 2-arylthiazolidinones as selective analgesic N-type calcium channel blockers. Calcium Channel blockers is the molecular target responsible for the treatment of neuropathic and inflammatory pain. The 3D-QSAR study based on the principle of the alignment of pharmacophoric features by PHASE module of Schrodinger suite has been carried out on the same set of calcium channel blockers. Statistically significant 3-D QSAR model (R²=0.9288) were generated using 21 molecules in the training set. The predictive ability of model was determined using a randomly chosen test set of 6 molecules which gave predictive correlation coefficients (R²_pred) of 0.946 for 3-D models, indicating good predictive power. PHASE pharmacophore hypothesis AAHR.13 may correspond very closely to the interactions recorded in the active site of the ligand bound complex.  These studies produced models with high correlation coefficient and good predictive abilities. Docking studies were also carried out wherein these analogues were docked into the active sites of COX-2 to analyze the receptor-ligand interactions that confer selectivity for COX-2. Compound 2 have the highest dock score (-7.28).  In the active site, there are some strong hydrogen-bonding interactions observed between residues GLU67, ALA103, ASP96, SER184 and ASP22. Additionally a correlation of the quantitative structure –activity relationship data and the docking results is found to validate each other and suggest the importance of the binding step in overall drug action.

INTRODUCTION:Ion channel blockers are attractive drug target for an expending range of therapeutic indications and voltage dependent calcium channels play important roles in critical biological processes. Based on their pharmacological and electro-physiological properties the channels have been classified into several subtypes:  N, L, P, R and T type calcium channel. These targets for therapeutic intervention play specialized roles in cellular functions ^{1, 2}.

N type channels are located primarily at pre-synaptic nerve terminals and mediate spinal transmission of pain signals from the periphery to the central nervous system by modulating release of nociceptive neuro-transmitters and neuropeptides ³.

Selective calcium channels blockers in particular are now emerging as prospective therapeutics for the treatment of neuropathic and inflammatory pain. Mice lacking N-type channels show suppressed response to painful stimuli that induce inflammation, and show reduced neuropathic pain symptoms. This provides evidence that N-type channels may be essential for development of neuropathic pain, that is pain associated with nerve injury, and control of N-type activity is important in the management of pain ⁴.

The recently approved pain drug Ziconotide is a potent blocker of N-type calcium channels. Pre- clinical and clinical studies of Ziconotide conclude that selective blockade of the N-ype channel is effective in reducing inflammatory and neuropathic pain in humans ⁵.

Cyclooxygenase (COX) is an enzyme (EC 1.14.99.1) that is responsible for formation of important biological mediators called prostanoids, including prostaglandins, prostacyclin and thromboxane. Pharmacological inhibition of COX can provide relief from the symptoms of inflammation and pain. Non-steroidal anti-inflammatory drugs, such as aspirin and ibuprofen, exert their effects through inhibition of COX. The enzyme contains two active sites: a heme with peroxidase activity, responsible for the reduction of PGG₂ to PGH₂, and a cyclooxygenase site, where arachiodonic acid is converted into the hydroperoxy endoperoxide prostaglandin G₂ (PGG₂).

The main aim of this study is to construct a pharmacophore model based on key chemical features of compounds with COX-2 inhibitory activity and then is to find the new lead candidate molecules through databases searching using the pharmacophore model. In our earlier studies, the key chemical groups of thiazolidinones analogues were extensively studied (QSAR) analysis methods and their binding modes were established by molecular docking studies.

Unlike the previous QSAR study, which is based on single scaffold containing molecules, we have generated pharmacophore model that can provide a rational hypothetical picture of the primary chemical features responsible for activity. The best pharmacophore model is expected to provide useful knowledge for developing new potentially active candidates targeting the COX-2.

Molecules used for Pharmacophore modeling and Docking studies:

Computational Details: The chemical structures of the 27 analogues are given in table 1. The biological activities were converted into the corresponding –log (IC₅₀) values. The observed analgesic activities for the training set and test set are given in table 2 and 3, respectively. The various 3D-QSAR studies were carried out using Phase while the docking studies were carried out using Glide respectively.

Pharmacophore Modeling: The common pharmacophore hypothesis (CPH) and alignment based on it was carried out using PHASE, version 3.0, 2008 (Schrödinger, LLC, New York, NY) The structures were imported from the project table in the “Develop Pharmacophore Hypothesis” panel and geometrically refined (cleaned) using Ligprep. Conformations were generated by the MCMM/LMOD method using a maximum of 2,000 steps with a distance-dependent dielectric solvent model and an OPLS-2005 force field. All the conformers were subsequently minimized using truncated Newton conjugate gradient minimization up to 500 iterations.

For each molecule, a set of conformers with a maximum energy difference of 30 kcal/mol relative to the global energy minimum conformer was retained.  Pharmacophore features—hydrogen bond acceptor (A), hydrogen bond donor (D), hydrophobic group (H), aromatic ring (R) were selected for creating sites. Common pharmacophoric features were then identified from a set of variants—a set of feature types that define a possible pharmacophores — using a tree-based partitioning algorithm. These CPHs were examined using a scoring function to yield the best alignment of the active ligands using an overall maximum root mean square deviation (RMSD) value of 1.2 Å with default options for distance tolerance. The quality of alignment was measured by a survival score, defined as:

S = W_siteS_site +W_vecS_vec +W_volS_vol +W_selS_sel +W^m_rew

Where W are weights and S are scores; Ssite represents alignment score, the RMSD in the site point position; S_vec represents vector score, and averages the cosine of the angles formed by corresponding pairs of vector features in aligned structures; S_vol represents volume score based on overlap of van der Waals models of non hydrogen atoms in each pair of structures; and S_sel represents selectivity score, and accounts for what fraction of molecules are likely to match the hypothesis regardless of their activity towards the receptor. W_site, W_vec, W_vol, W_rew have default values of 1.0, while W_sel has a default value of 0.0. In hypothesis generation, default values have been used. W^m_rew represents reward weights defined by m−1 where m is the number of actives that match the hypothesis.

QSAR Model Building: The in vitro biological data of a series of 27 calcium channel blockers was used for this study. 21 molecules were used as training set for developing pharmacophore model and 6 molecules were used as test set. The ligands of training set were used for generating QSAR model. This pharmacophore based model was generated keeping grid spacing 1Å.  

Validation of the Generated Pharmacophore Model: Validation is a crucial aspect of pharmacophore design, particularly when the model is built for the purpose of predicting activities of compounds in external test series (35). A graph was plotted between experimental activity and estimated activity of the test set of the data. The model was validated by determining correlation between experimental activity and estimated activity of the test set of the data.

Structure based Studies: Docking studies were carried out using Glide and affinity for 2-arylthiazolidinones molecules in form of COX-2. The molecular target taken into consideration is the nuclear COX-2 whose crystallographic structure is available on the PDB database as 3HS5. Docking is a method which predicts the preferred orientation of one molecule to a second when bound to each other o form a sable complex. Docking of a ligand is typically achieved by generating a number of orientation ( or poses) of a ligand within the active sites, and scoring of poses to identify one or more that closely approximate the bioactive conformation determined by X-ray crystallography.

The docking analysis was firstly carried out on the cyclooxygenase binding site in the 3HS5 obtained from the RCSB Protein Data Bank, http://www.rcsb.org/pdb, where the residues were bonded more closely to the arachiodonic acid, co-crystallized with COX-2.  The arachiodonic acid was extracted from the complex, and, during the calculation, the active site was determined within a cubic box centered on the co-crystallized ligand (coordinates X, Y and Z).  Within the pre-defined region, the grid of probe atom interaction energies were computed with a default resolution and the ligands were then docked by genetic algorithm and the 50 lowest energy structures were stored for further analysis.

Steps undertaken in Docking Study:

Step I- Data set for Analysis: The analgesic activity of 27 compounds was measuredwith IC₅₀ values of upto 0.2 µM in an in vitro IMR assay. These inhibitors along with conformational energy were shown in Table 2.

Step II- Protein Preparation: Protein preparation was performed by using protein preparation Wizard of Maestro software. The protein structure 3HS5 was taken from PDB. The A chain was treated to add missing hydrogen, assign proper bond orders (i.e., breaking bonds to metal and correct the formal charge 237 on it and neighboring atoms), and to delete water compounds that were more than 5Ǻ from the heterogeneous groups. The H-bonds were optimized using sample orientations. All the polar hydrogens were displayed. Finally, the protein structure was minimized to the default Root Mean Square Deviation (RMSD) value of 0.30.

Step III- Ligand Preparation: Ligand preparation was accomplished on all the 15 compounds using Ligprep module to clean the structure and generate tautomers as described in the PHASE CPH generation studies.

Step IV- Receptor Grid Generation: Receptor was defined and the co-crystallized ligand was differentiated from the active site of receptor A chain. The atoms were scaled by Vander Waals radii of 1.0 Ǻ with the partial atomic charge less than 0.25 defaults. The active site was defined as an enclosing box at the centroid of the workspace ligand as selected in the receptor folder. The ligands similar in size to the workspace ligand were allowed to dock into the active site. No constraints either positional, H-bonding or hydrophobic were defined.

Step V- Ligand Docking: Ligand docking was performed using OPLSAA force field. The receptor grid defined in the receptor grid generation folder was selected for the docking of ligands prepared using Ligprep. Flexible docking was performed using the Extra Precision (XP) feature of Glide module. The Vander Waals radii was scaled using a default scaling factor of 0.80 and default partial cutoff charge of 0.15 to decrease the penalties for close contacts. The core pattern comparison and similarity mode were not used since the aim was to study the binding of ligands to the active site. The constraints to defined ligand – receptor interactions were not set. The structure output format was set to pose viewer file so as to view the output of the resulting docking studies from pose-viewer. The XP-Glide predicted pose of active compounds (4h) in the Scheme are shown in Figure 1.

Step VI- Viewing Docking Results: It was done using pose-viewer. The H-bonds and Vander Waals contacts (good, bad and ugly) to the receptor were visualized using default settings to analyze the binding modes of the ligands to receptor.

TABLE 2: EXPERIMENTAL –LOG (IC₅₀) AND CORRESPONDING MODEL PREDICTED VALUES OF MOLECULES USED IN TRAINING SETS

TABLE 3: EXPERIMENTAL –LOG (IC₅₀) AND CORRESPONDING MODEL PREDICTED VALUES OF MOLECULES USED IN TEST SETS

FIGURE 1: RELATION BETWEEN EXPERIMENTAL AND PREDICTED ACTIVITY VALUES OF TEST SET MOLECULES USING MODEL AAHR-1

FIGURE 2: RELATION BETWEEN EXPERIMENTAL AND PREDICTED ACTIVITY VALUES OF TRAINING T SET MODEL AAHR-1

RESULTS AND DISCUSSION: The model generated by Phase is summarized in Table 4. They were analyzed for their predictive ability for training set as well as test set molecules. Final models were selected primarily based on the values of better cross- validated r², predictive r² and SD values of test set molecules.

A summary of statistical data of the best CPHs, labeled AAHR-1 with their fitness scores is listed in Table 4. Good and consistent external predictivity was observed for AAHR-1. AAHR1 showed a good r² value for the training set (0.615). A good Pearson-R value of 0.6308 was also observed. Actual and predicted values of the training set and test set molecules are given in Table 2 and table 3. The features represented by this model (Figure 3) were two hydrogen bond acceptor (A), and one aromatic ring (R) and one hydrophobic(H).

TABLE 4: PARAMETERS OF HYPOTHESES

 Parameters of 4 Featured Pharmacophore Hypotheses:

Entry    Site1    Site2    Distance

A2        A3        3.616

A2        H6        7.136

A2        R7        5.750

A3        H6        5.094

A3        R7        2.846

H6        R7        2.909

Entry    Site1    Site2    Site3    Angle

A3        A2        H6        42.5

A3        A2        R7        23.8

H6        A2        R7        23.0

A2       A3       H6       108.8

A2       A3       R7       125.3

H6        A3         R7        28.1

A2        H6        A3        28.7

A2        H6        R7        50.7

A3        H6        R7        27.4

A2        R7        A3        30.9

A2       R7       H6       106.3

A3       R7       H6       124.5

FIGURE 3: BEST PHARMACOPHORE MODEL AAHR-1 ALIGNED WITH THE COMPOUND NO.23 OF TEST SET. PHARMACOPHORE FEATURES ARE COLOR CODED: HYDROGEN BOND ACCEPTOR (PINK), AROMATIC RING (ORANGE), HYDROPHOBIC (GREEN)

FIGURE 4: 3-D QSAR MODEL BASED ON COMPOUND NO 23 OF TEST SET ILLUSTRATING HYDROGEN BOND DONOR FEATURE

FIGURE 5: 3-D QSAR MODEL BASED ON COMPOUND NO 27 OF TRAINING SET ILLUSTRATING HYDROPHOBICITY FEATURE

The 3-D QSAR model based on molecule 23 of the test set using hydrogen bond donor feature is shown in figure 4.  In order to gain additional insight about the nature of interaction of these ligands with COX-2 receptor, we visualized and analyzed the 3D-QSAR model of compound no. 23 of test set using hydrogen bond donor and hydrophobic features.

Fig. 4 shows the 3-D QSAR model illustrating the H-Bond donor features. The blue region shows where H bond donor group will increase the COX-2 antagonistic activity where the red region shows where such groups decrease the activity. H-bond donor group at this blue region will increase the activity. The red region at the thiazolidinone ring shows the importance of hydrogen bond acceptor group for activity. The NH₂ and 0H group are hydrogen bond donor and decrease the activity at this red region. All halogen (Cl, Br, and F) are hydrogen bond acceptor in nature. So the compounds 1-4 have halogen at these thiazolidinone rings and have better activity as compare to other compounds. Any bulky group substitution on pyridine ring have lowers the activity. The hydrophobic features shows red region on pyridine moiety. So compounds having bulky group substitution on pyridine ring have low activity.

REFERENCES:

1. Spafford, J.D.; Zamponi, G. W. Curr. Opin. Neurobiol. 2003, 13, 308.
2. Catterall, W. A. Cell Calcium 1998, 24, 307.
3. Cox, B.; Denyer, J. C. Expert Opin. Ther. Patents 1998, 8, 1237; Elmslie, K. S. J. Neurosci. Res. 2004, 75, 733.
4. Saegusa, H.; Kurihara, T.; Zong, S.; Kazuno, A.; Matsuda, Y.; Nonaka, T.; Han,W.; Toriyama, H.; Tanabe,. EMBO J. 2001, 20, 2349.
5. Miljanich, G. P. Curr. Med. Chem. 2004, 11, 3029; Matthews, E.A.; Dickenson, A. H. Pain 2001, 92, 235.
6. Stahl, S. M. J. Clin. Psychiatry 2004, 65, 1033.
7. Hara, H.; Morita, T.; Sukamoto,T.; Cutrer, F. M. CNS Drug Rev. 1995, 1, 204.
8. Pauwels, P. J.; Leysen,J.E.;Janseen, P. A. J. Life Sci. 1991, 48, 1881.
9. Thomas, A. D.; Asokan, C. V. J. Chem. Soc. Perkin trans. 1 2001, 20, 2583.
10. Hariprasad, V., Kulkarni, V.M., J. Comp. Aided Mol. Des., 1996, 10, 284.
11.  Wermuth, C.G.; Ganellin, C.R.; Lindberg, P.; Mitscher, L.A. (1998). "Glossary of terms used in medicinal chemistry". Pure Appl. Chem. (International union of pure and applied chemistry).
12.  www.diabetesselfmanagement.com/article/ppar _agonist/-
13.  Zimmet P Z, Alberti K G M M & Shaw J. Global and social implication of the diabetes epidemic, Nature, 414 (2001) 782-787.
14.  J.W.F. Elte, J.F. Blickle, European Journal of Internal Medicine 18(2007) 18-25.
15. Skyler J S, Diabetes mellitus: Pathogenesis and treatment strategies, J Med Chem, 47 (2004) 4113-4117.
16. Mantri, P.; Witiak, D. Inhibitors of COX and 5-Lipoxygenase. Curr. Med. Chem. 1994, 1, 328-355.
17. Needleman, P.; Isackson, P.C. Selective Inhibitors of COX 2. Sci. Med. 1998, 5, 26-35.
18. SYBYL 6.6 Molecular Modeling Software; Tripos Associates Inc.; 1699 South Hanley Road, St. Louis, MO 63144.
    

    ---