Docking and 3D QSAR Studies on p38 α MAP Kinase Inhibitors

: The p38 signaling cascade has emerged as an attractive target for the design of novel chemotherapeutic agents for the treatment of inflammatory diseases. Three dimensional quantitative structure-activity relationship (3D-QSAR) studies were performed on a series of 25, 2-aminothiazole analogs as inhibitors of p38 α mitogen activated protein (MAP) kinase. The docking results provided a reliable conformational alignment scheme for the 3D-QSAR model. The 3D-QSAR model showed very good statistical results namely q 2 , r 2 and r 2pred values for both comparative molecular field analysis (CoMFA) and comparative molecular similarity indices analysis (CoMSIA). The CoMFA and CoMSIA models & docking results provided the most significant correlation of steric, electrostatic, hydrophobic, H -bond donor, H -bond acceptor fields with biological activities and the provided values were in good agreement with the experimental results. The information rendered from molecular modeling studies gave valuable clues to optimize the lead and design new potential inhibitors.


Introduction
p38 MAP kinase (MAPK) is an intracellular, soluble proline targeted serine-threonine kinase which belongs to a large family of proteins that include extra cellular regulated kinases (ERKs) and c-Jun N-terminal kinases (JNKs) [688537], [688545], [688587].Till date only four p38 isoforms have been identified, namely p38α 1 , p38β, p38γ and p38δ.The ubiquitous p38α has been the most extensively studied and is believed to be the most physiologically relevant in the regulation of the inflammatory response.p38 MAP kinase is a member of Ser/Thr kinase family of the mitogen activated protein (MAP) kinase super family.MAP kinase plays a very important role in diseases, such as asthma, osteoarthritis, rheumatoid arthritis (RA) and chronic inflammatory autoimmune diseases.Therefore, the inhibition of MAP kinase would potentially prevent the underlying pathophysiology in the inflammatory diseases 2,3 , which makes it an attractive target for drug discovery.
Rheumatoid arthritis (RA) leads to destruction of cartilage and deformation of bones 4,5 .The injury of inflammation causes release of inflammatory mediators such as Tumor necrosis factor-α and Interlucin-1β. 6p38α regulates the biosynthesis of these two proinflammatory cytokins and these two cytokins 7 are associated with the progression of rheumatoid arthritis.The inhibitors of p38α cease their progression through the action of competitive inhibition.

Experimental
A series of 25, p38α inhibitors, published by John Hynes et al 8 were considered to derive 3D-QSAR models 9 .The reported IC 50 values were converted into their corresponding pIC 50 values (Table 1).All the molecular modeling calculations were performed on a Linux operating system.Three dimensional structure building and all modeling studies were performed using the SYBYL 6.9 molecular modeling program package 10 .Gasteiger-Hückel 11,12 charges were assigned and then energy minimization of each molecule was performed using the conjugate gradient method and Tripos FF standard force field with a distance-dependent dielectric function.The minimization was terminated when the energy gradient convergence criterion of 0.001 kcal mol −1 •Å −1 was reached 13 .The crystal structure of MAP kinase in complex with BMS-640994 inhibitor (pdb id: 3bx5) was downloaded from the protein data bank.The protein was minimized using the conjugate gradient method by applying gasteiger-hückel charges and tripos FF standard force field with a distance-dependent dielectric function.The minimization was terminated when the energy gradient convergence criterion of 0.05 kcal mol −1 •Å −1 was reached.FlexX in SYBYL 6.9 was used for molecular docking of the 25 molecules into the protein active site containing Lys53, Glu71, Thr106, Met109 and Asp168.The crystal structure ligand was also docked and its RMSD was calculated to validate the docking process.The analysis of dock poses of all the molecules showed similar hydrogen bond interaction with the active site residues.

Alignment
In the 3D-QSAR studies, alignment is the most sensitive and vital part since the interaction energies depend upon the positioning of molecules in the 3D fixed lattice.Two alignment methods were employed, in the first method the most active compound 2g was selected as template.Molecules were aligned by the ALIGN DATABASE command available in SYBYL 6.9 using common substructure (Figure 2).This adjusts the geometry of the molecules such that it's steric and electrostatic fields match the fields of the template molecule.The resulting alignment is shown in Figure 3(a).In the second method, the active conformation obtained from docking process of most active molecule in the series compound 2g was used as template for alignment in ALIGN DATABASE taking the substructure that is common to all (Figure 2).The resulting alignment is shown in Figure 3(b).This alignment process is a standard rigid RMSD overlay of selected common structural motive, where the docked poses of the most active molecule was used as a template on to which the dataset was aligned.Standard tripos force field was employed for the CoMFA 14,15 and CoMSIA 16 analysis.A 3D cubic lattice overlapping all entered molecules and extended by at least 4 Å in each direction with each lattice intersection of a regularly spaced grid of 2.0 Å was created.The steric and electrostatic parameters were calculated in case of the CoMFA fields, while hydrophobic, H-bond acceptor and H-bond donor parameters in addition to steric and electrostatic were calculated in case of the CoMSIA fields at each lattice.A sp 3 hybridized carbon atom was used as a probe atom to generate steric (Lennard-Jones potential) field energies and a charge of +1 to generate electrostatic (Coulombic potential) field energies.A distance dependent dielectric constant of 1.00 was used.The steric and electrostatic fields were truncated at +30.00 kcal mol −1 .The similarity indices descriptors were calculated using the same lattice box employed for CoMFA calculations, using sp 3 carbon as a probe atom with a +1 charge, +1 hydrophobicity and +1 H-bond donor and +1 H-bond acceptor properties.
A partial least squares 17 regression was used to generate a linear relationship that correlates changes in the computed fields with changes in the corresponding experimental values of biological activity (pIC 50 ) for the data set of ligands.Biological activity values of ligands were used as dependent variables in a PLS statistical analysis.The column filtering value (s) was set to 2.0 kcal mol −1 to improve the signal-to-noise ratio by omitting those lattice points whose energy variations were below this threshold.Cross-validations were performed by the leave-one-out (LOO) procedure to determine the optimum number of components (ONC) and the coefficient q 2 .The optimum number of components obtained is then used to derive the final QSAR model using all of the training set compounds with noncross validation and to obtain the conventional correlation coefficient (r 2 ).To validate the CoMFA and CoMSIA derived models, the predictive ability for the test set of compounds (expressed as r 2 pred ) was determined by using the following equation: r 2 pred = (SD -PRESS)/SD Where, SD is the sum of the squared deviations between the biological activities of the test set and mean activities of training set molecules and PRESS is the sum of squared deviation between predicted and actual activity for every molecule in test set

Results and Discussion
The accuracy of a docking procedure was evaluated by determining how closely the lowest energy pose (binding conformation) predicted by the object scoring function (dock score), resembles the experimental binding mode as determined by x-ray crystallography for the cocrystallized ligand.The co-crystallized ligand was redocked into the active site to find out the difference between the localization of the inhibitor upon docking from the crystal structure and that was in very good agreement.The root mean square deviations between the predicted conformation and the observed x-ray crystallographic conformation of cocrystallized ligand 0.375 Å (Figure 4a), a value that suggests the reliability of FlexX docking in reproducing the experimentally observed binding mode for MAP Kinase inhibitor and the parameter set for the Glide docking is reasonable to reproduce the x-ray structure.The 3D QSAR -CoMFA and CoMSIA analysis were carried out using 2-aminothiazols series of MAP kinase inhibitors reported by Hynes, Jr. et.al. 8A total of 25 molecules were used for derivation of model, these were divided into a training set of 20 molecules and test set of five.The CoMFA and CoMSIA statistical analysis is summarized in Table 2. Statistical data shows q 2 loo 0.72 and 0.76 for ligand and receptor based CoMFA, 0.76 and 0.69 for ligand and receptor based CoMSIA models respectively, r 2 ncv of 0.97 and 0.99 for ligand and receptor based CoMFA, 0.92 and 0.93 for ligand and receptor based CoMSIA models respectively, which indicates a good internal predictive ability of the models.To test the predictive ability of the models, a test set of five molecules excluded from the model derivation was used.The predictive correlation coefficient r 2 pred of 0.95 and 0.88 for ligand and receptor based CoMFA, 0.81 and 0.69 for ligand and receptor based CoMSIA models respectively indicating good external predictive ability of the models.The graph for the actual and predicted pIC50 values for training set and test of CoMFA and CoMSIA studies are shown in Figure 5.The predicted activity and glide scores for the molecules are provided in Table 3.
To visualize the information content of the derived 3D-QSAR model, CoMFA and CoMSIA contour maps were generated.The contour plots are the representation of the lattice points and the difference in the molecular field values at lattice point is strongly connected with difference in the receptor binding affinity.Molecular fields define the favorable or unfavorable interaction energies of aligned molecules with a probe atom traversing across the lattice grid points surrounding the molecules.The 3D colored plots suggest the modification required to design new molecules.
The contour maps of CoMFA denote the region in the space where the aligned molecules would favorably or unfavorably interact with the receptor, while the CoMSIA contour maps denote those areas within the specified region where the presence of a group with a particular physicochemical activity binds to the receptor.The CoMFA/ CoMSIA results were graphically interpreted by field contribution maps using the 'STDEV * COEFF' field type.Receptor Based CoMSIA   Steric and electrostatic contour map 12,13 for both ligand and receptor based CoMFA illustrates that the n-propyl group of thiazole ring is located in between two big green regions which is favored region for sterically bulkier groups, suggesting an increase in bulkiness at this position improves the activity.A polyhedron of yellow contour is seen around thiazole ring suggesting other substitutions on this ring will reduce the potency, similar information is obtained from receptor based CoMFA.Red contour for the electrostatic field is observed at the carbonyl group of the amide linkage showing the evidence that the more electronegative substituents would be favorable at this region for improving the inhibitory activity.This prediction is consistent with the observation that the hydrogen bond can be built between the carbonyl group and the NH group of Asp168 The blue contour around the two amino groups suggest that the electron donating group at this position is beneficiary this can be explained by the hydrogen bond interactions of these amino groups with Glu71.The steric and electrostatic fields of CoMSIA, as shown in Figure 7, were generally in accordance with the field distribution of the CoMFA maps (Figure 6).Besides the structural features already mentioned in the CoMFA steric field analysis there was a yellow contour seen near the propyl group from where the propyl group is moving away from it.A green contour is seen near the ethyl substitution at the phenyl ring suggesting the increase in bulkiness will improve potency.Purple and cyan contours represent favorable and unfavorable hydrogen bond acceptor groups, orange and blue contours represent favored and disfavored hydrogen bond donor groups respectively.Purple contour on the carbonyl group on the terminal amide present on the phenyl ring represents a hydrogen bond acceptor favored region where it also forms a hydrogen bonding with amino group of Asp168 suggesting a hydrogen bond acceptor group is more favored at this region.The NH group of the same amide is overlapped with a cyan contour in the receptor based CoMSIA representing an acceptor disfavored region.An orange region near the NH of the same amide in ligand based contour shows a hydrogen bond donor favored region and it is confirmed by the hydrogen bond interaction with carbonyl of Glu71 of the active site amino acid (Figure 4b).

Ligand Based CoMFA
An orange contour is seen on the other side of the molecule near the NH of the amide group adjacent to thiazole ring in both ligand and receptor based contours depicting a hydrogen bond donor favored group will increase the inhibitory potency and this is evident from docking studies where a hydrogen bond interaction is seen with the carbonyl group of Met 108 active site amino acid.As depicted in the Figure 9a for ligand based and 9b for receptor based CoMSIA using compound 2g as reference molecule, where magenta and gray indicate the hydrophobic and hydrophilic group favored regions.A magenta area is placed near the thiazole and the phenyl ring indicating that hydrophobic substituents may be well tolerated in that region.Minor grey hydrophobic contour near the amide on the phenyl ring indicates that groups with hydrophilic characters are preferred in this position.From QSAR and docking studies it is evident that the amide group on the phenyl ring and the NH group adjacent to the thiazole ring have to be incorporated while designing new molecules.

Conclusion
3D QSAR is widely employed to develop new molecules that have an improved biological property.CoMFA and CoMSIA methodologies were used to build models for P38α MAP Kinase inhibitory activity of the 2-aminothiazole analogs.Based on the detailed contour map analysis, improvement in P38α MAP kinase binding affinity can be achieved through conformationally restricted substitution at n-propyl group of thiazole ring, maintaining the hydrophobic character with less steric hindrance at these regions.The CoMFA and CoMSIA contour maps along with the docking results offered enough information to understand the structure-activity relationship and identified structural features influencing the inhibitory activity.The correlation of the results obtained from 3D-QSAR and docking studies can serve as a useful guideline for further modification of 2-aminothiazole analogs as inhibitors of P38α mitogen activated protein (MAP) kinase.

Figure 1 .
Figure 1.Core molecules of MAP kinase inhibitorsThe crystal structure of MAP kinase in complex with BMS-640994 inhibitor (pdb id: 3bx5) was downloaded from the protein data bank.The protein was minimized using the conjugate gradient method by applying gasteiger-hückel charges and tripos FF standard force field with a distance-dependent dielectric function.The minimization was terminated when the energy gradient convergence criterion of 0.05 kcal mol −1 •Å −1 was reached.FlexX in SYBYL 6.9 was used for molecular docking of the 25 molecules into the protein active site containing Lys53, Glu71, Thr106, Met109 and Asp168.The crystal structure ligand was also docked and its RMSD was calculated to validate the docking process.The analysis of dock poses of all the molecules showed similar hydrogen bond interaction with the active site residues.

Figure 2 .Figure 3 .
Figure 2. Common substructure used for the alignment

Figure 4
a). Superimposition of crystal structure pose (orange) on docked pose (blue) of cocrystallized ligand.The RMS deviation is 0.375Å.b) The hydrogen bond interaction of aminoacids with best active ligand

Figure 6 (
Figure 6 (a,) shows the contour maps derived from the CoMFA PLS model generated from ligand based alignment and 6(b) from Receptor based alignment.The most potent analogue, compound 2g, was embedded in the map to demonstrate its affinity for the steric and electrostatic regions of inhibitors.The areas of yellow indicate regions of steric hindrance to activity, while green areas indicate a steric contribution to potency.The blue regions indicate positive electrostatic charge potential associated with increased activity, while red region show negative charge potential.All of the contours represented the default 80 and 20% level contributions for favored and disfavored regions, respectively.Steric and electrostatic contour map12,13 for both ligand and receptor based CoMFA illustrates that the n-propyl group of thiazole ring is located in between two big green regions which is favored region for sterically bulkier groups, suggesting an increase in bulkiness at this position improves the activity.A polyhedron of yellow contour is seen around thiazole ring suggesting other substitutions on this ring will reduce the potency, similar information is obtained from receptor based CoMFA.Red contour for the electrostatic field is observed at the carbonyl group of the amide linkage showing the evidence that the more electronegative substituents would be favorable at this region for improving the inhibitory activity.This prediction is consistent with the observation that the

Figure 6 .
CoMFA steric standard deviation (S.D. * coefficient) contour maps illustrating steric and electrostatic features in combination with compound 2g (a) ligand based (b) Receptor based.Green contours show favorable bulky group substitution at that point while yellow regions show disfavorable bulky group for activity.Red contours indicate negative charge favoring activity, whereas blue contours indicate positive charge favoring activity.

Figure 7 .
CoMSIA S.D. * coefficient contour maps illustrating steric, electrostatic features in combination with compound 2g.a) Ligand base and (b) Receptor base.The green contour indicates a sterically favored region; yellow maps calls for a reduction of this potential to improve activity.Blue indicates a positive charge preferred region to improve activity; Red indicates a negative charge preferred region to improve activity.The hydrogen bond acceptor and donor field contour map of CoMSIA showing Figure 8a for ligand based and Figure 8b for receptor based using compound 2g as a reference molecule.

Figure 8
CoMSIA S.D. * coefficient contour maps illustrating acceptor and donor features in combination with compound 2g.a) Ligand based and (b) Receptor based.The purple contour for H-bond acceptor group increases activity, cyan indicates the disfavored region; the orange and blue contours represent favored and disfavored hydrogen bond donor groups respectively

Figure 9 .
CoMSIA S.D. * coefficient contour maps illustrating hydrophobic features in combination with compound 2g.a) Ligand based and (b) Receptor based.The magenta contour for hydrophobic favored region, gray indicates the hydrophilic favored region

Table 1 .
Dataset used for 3D-QSAR analysis with corresponding IC50 and pIC50 values

Table 2 .
Statistical summary of PLS analysis

Table 3 .
Data set used for 3D-QSAR analysis with corresponding actual and predicted