3D-QSAR and Molecular Docking Studies on Fused Pyrazoles as p38α Mitogen-Activated Protein Kinase Inhibitors

The p38α mitogen-activated protein kinase (MAPK) has become an attractive target for the treatment of many diseases such as rheumatoid arthritis, inflammatory bowel disease and Crohn’s disease. In this paper, 3D-QSAR and molecular docking studies were performed on 59 p38α MAPK inhibitors. Comparative molecular field analysis (CoMFA) and comparative molecular similarity indices analysis (CoMSIA) were applied to determine the structural requirements for potency in inhibiting p38α MAPK. The resulting model of CoMFA and CoMSIA exhibited good r2cv values of 0.725 and 0.609, and r2 values of 0.961 and 0.905, respectively. Molecular docking was used to explore the binding mode between the inhibitors and p38α MAPK. We have accordingly designed a series of novel p38α MAPK inhibitors by utilizing the structure-activity relationship (SAR) results revealed in the present study, which were predicted with excellent potencies in the developed models. The results provided a useful guide to design new compounds for p38α MAPK inhibitors.


Introduction
The p38 mitogen-activated protein kinase (MAPK), a serine/threonine kinase, plays a crucial role in biosynthesis of pro-inflammatory cytokines including tumor necrosis factor α (TNFα) and

OPEN ACCESS
interleukin-1 β (IL-1β) [1]. The p38 MAPK exists in four isoforms (α, β, γ and δ). Expression of these isoforms varies across cell types of the immune system, and it was revealed by previous studies that the predominant isoform involved in inflammation is p38α [2,3]. The p38α MAPK is activated by a range of environmental stimuli such as TNFα, IL-1β and stress. Activation of p38α leads to the up-regulation of both TNFα and IL-1β, resulting in many chronic inflammatory diseases, such as rheumatoid arthritis (RA), inflammatory bowel disease, Crohn's disease, and chronic obstructive pulmonary disease (COPD), as well as the other inflammatory disorders [4][5][6][7][8].
The proven ability of p38 MAPK to efficiently regulate both the release and the activity of these pro-inflammatory cytokines has prompted many pharmaceutical groups to pursue inhibitors of p38 MAPK for the potential treatment of various inflammatory diseases [9][10][11]. A series of fused pyrazole derivatives with potent, selective and orally available inhibition activities towards p38 MAPK were reported [1][2][3]. In this paper, 3D-QSAR and molecular docking studies were performed on these inhibitors. Along with the molecular docking, 3D-QSAR approaches including comparative molecular field analysis (CoMFA) and comparative molecular similarity analysis (CoMSIA) could offer more insight into understanding the structure-activity relationship of these inhibitors, and thus could more effectively direct the design of new potential inhibitors. Based on the structure-activity relationship revealed in the present study, a series of novel p38α MAPK inhibitors were accordingly designed.

Dataset for Analysis
All the compounds and associated biological activities were selected from literature [1][2][3] reported by the same research group. These fused pyrazole derivatives were divided into a training set of 46 compounds and a test set of 13 compounds. The inhibitory data was reported as IC 50 towards p38α MAPK, and the IC 50 values were converted into corresponding pIC 50 by taking Log (1/IC 50 ). The pIC 50 values were used as the dependent variables in all the models subsequently developed. Structures and associated inhibitory activities are shown in Tables 1 and 2.

Molecular Modeling and Alignment
All computational studies were performed using SYBYL 8.1 molecular modeling software from Tripos, Inc., US [12]. 3D structures of all compounds were constructed using the Sketch Molecule module. Structural energy minimization was performed using the standard Tripos molecular mechanics force field and Gasteiger-Hückel charge [13]. The fragment that was used as the common structure is shown in Figure 1. The compound 9 was selected as a template molecule, which was one of the most active derivatives in the dataset. The aligned molecules of the training set are shown in Figure 2.

CoMFA and CoMSIA Setup
To calculate the CoMFA and CoMSIA fields, a 3D cubic lattice with grid spacing of 1 Å and extending to 4 Å units beyond the aligned molecules in all directions was created automatically by SYBYL [12]. In the CoMFA method, a sp 3 hybridized carbon atom with a charge of 1e served as the probe atom to calculate steric and electrostatic fields, in which their energy values were truncated at 30 kcal/mol [14]. The CoMSIA method, incorporating steric, electrostatic, hydrophobic, hydrogen bond donor and acceptor fields, was carried out using a probe atom with radius 1.0 Å, +1.0 charge, and hydrophobic and hydrogen bond properties of +1. The attenuation factor was set to the default value of 0.3 [15].

Regression Analysis and Models Validation
Partial Least Squares (PLS) was used to linearly correlate the CoMFA and CoMSIA fields to the pIC 50 values. The performance of both the CoMFA and CoMSIA models was evaluated using the leave-one-out (LOO) method [16]. PLS was conjunct with the cross-validation option to determine the optimum number of components (ONC), which were then used in deriving the final CoMFA and CoMSIA model without cross-validation. The ONC was the number of components resulting in the highest cross-validated correlation coefficient (r 2 cv ). After obtaining the optimum number of components, a PLS analysis was performed with no validation and column filtering 2.0 to generate the highest correlation coefficient (r 2 ) [17][18][19].

Predictive Correlation Co-Efficient (r 2 pred )
The predictive abilities were determined from a test set of 13 compounds that were not included in the training set. These molecules were aligned to the template and their pIC 50 values were predicted. The predictive correlation coefficient (r 2 pred ), based on the molecules of test set, was calculated by the equation shown below.
SD is the sum of the squared deviations between the inhibitory activities of the test set and mean activities of the training molecules and PRESS is the sum of squared deviations between predicted and actual activity values for each molecule in the test set [20][21][22].

Molecular Docking
To investigate the protein-ligand interactions, compound 9 was docked into the ATP-binding site of p38α MAPK. The Surflex-Dock, using an empirical scoring function and a patented search engine to dock ligands into a protein's binding site, was applied to study molecular docking [12]. The crystal structure of p38α MAPK was retrieved from the RCSB Protein Data Bank (PDB entry code: 3LHJ). The p38α MAPK structure was utilized in subsequent docking experiments without energy minimization. The ligands were docked into the corresponding protein's binding site by an empirical scoring function and a patented search engine in Surflex-Dock [12]. All ligands and water molecules were removed and the polar hydrogen atoms were added. Protomol, an idealized representation of a ligand that makes every potential interaction with the binding site, was used to guide molecular docking. The establishment of protomol demonstrates three behaviors: (a) Automatic: Surflex-Dock finds the largest cavity in the receptor protein; (b) Ligand: a ligand in the same coordinate space as the receptor; (c) Residues: specified residues in the receptor [10]. In this paper, the automatic docking was applied. To visualize the binding mode between the protein and ligand, the MOLCAD (Molecular Computer Aided Design) program was employed. MOLCAD calculates and displays the surfaces of channels and cavities, as well as the separating surface between protein subunits [12]. The MOLCAD program provides a variety of properties and maps onto protein surfaces. In this paper, the cavity depth, electrostatic and lipophilic potential surfaces were established by MOLCAD. Other parameters were established by default in the software.

CoMFA and CoMSIA Analysis
The statistical parameters obtained from the CoMFA model are listed in Table 3 Table 2. The relationship between experimental and predicted pIC 50 values of the training set and test set compounds in CoMFA and CoMSIA are illustrated in Figure 3(a),(b).

Graphical Interpretation of CoMFA and CoMSIA
One of the attractive features of the CoMFA and CoMSIA models is the visualization of the results as 3D coefficient contour maps. To visualize the information content of the derived 3D-QSAR model, CoMFA and CoMSIA contour maps were generated to rationalize the regions in 3D space around the molecules where changes in each field were predicted to increase or decrease the activity. The CoMFA steric and electrostatic contour maps that are shown in Figure 4 use compound 9 as a reference structure. In Figure 4(a), the green contours represent regions of high steric tolerance (80% contribution) while the yellow contours represent regions of low steric bulk tolerance (20% contribution). In Figure 4(b), the electrostatic field is indicated by blue (80% contribution) and red (20% contribution) contours, which reveal the regions where electron-donating group and electron-withdrawing group would be favorable, respectively.
As shown in Figure 4(a), the yellow contour near R 1 position indicates that bulky groups would decrease the potency. Comparing compound 27 with 24, their activity discrepancies can be explained by this yellow contour. A huge green contour around the R 3 and R 4 position suggested that bulkier groups would be favored. Most of the derivatives possessed a relatively bulkier methyl substituent at the R 3 position, compounds 17-19, which had minor groups (e.g., H, F), showed significantly decreased activities. In Figure 4(b), the blue contour near the R 3 position indicates that electron-donating groups may increase the activity. This may explain why compounds 17-19, without an electron-donating substituent at this position, were the most inactive derivatives. A red contour near the R 4 position demonstrated that electron-withdrawing groups would benefit the activity, compound 9 with an electron-withdrawing substituent (-F) at R 4 showed significantly increased activity. Two blue contours around the R 5 position strongly revealed that an electron-donating group would be favorable. Most of the compounds in the database possessed an electron-donating cyclopropyl substituent at R 5 , compounds 17-22 and 26-without an electron-donating group at this position-showed significantly decreased activities.
The CoMSIA steric, electrostatic, hydrophobic, hydrogen bond donor and acceptor field contour maps are shown in Figure 5 using compound 9 as a reference structure. The CoMSIA steric and electrostatic field contour maps were almost similar to the corresponding CoMFA contour maps.
The hydrophobic field contour map is shown in Figure 5(c), white (20% contribution) and yellow (80% contribution) contours highlight areas where hydrophilic and hydrophobic properties were favored. The yellow contour around the C-6 position revealed the extreme importance of the phenyl group at this position; replacing it with a hydrophilic group may result in decreased activity. Two white contours near R 4 and N-15 positions indicated that hydrophilic substituents may be favored.
The hydrogen bond donor field contour map is illustrated in Figure 5(d), the cyan (80% contribution) and purple (20% contribution) contours indicate favorable and unfavorable hydrogen bond donor groups. The purple contour near the C′-2 position of the phenyl group in R 2 indicates that a hydrogen bond donor would be unfavorable. In general, compounds with a hydrogen bond acceptor substituent (e.g., -F) at the C′-2 position showed better activities than those without a hydrogen bond acceptor group. The carbonyl group in the C-14 position was surrounded by a huge purple contour, suggesting it may serve as a hydrogen bond acceptor and that removing it may result in decreased potency.
The hydrogen bond acceptor field contour map is depicted in Figure 5(e), the purple (80% contribution) and red (20% contribution) contours favorable and unfavorable positions for hydrogen bond acceptors. The purple contour near the carbonyl group in the C-14 position revealed a hydrogen bond acceptor substituent at this position would benefit the activity. A red contour near N-15 position revealed the importance of the hydrogen bond donor -NH group.

Docking Analysis
To explore the interaction mechanism between these inhibitors and the receptor, compound 9 was selected for more detailed analysis. The MOLCAD Robbin surfaces structure of the ATP pocket of the p38α MAP within the compound 9 is shown in Figure 6. The -F at C′-2 position acted as a hydrogen bond acceptor and formed a H-bond with the -NH of the Met 109 residue, the carbonyl group at C-14 position also acted as a hydrogen bond acceptor by forming a H-bond with the -NH of the ASP168. The observations obtained from Figure 6 are in agreement with that of CoMSIA hydrogen bond donor and acceptor contour maps. The MOLCAD Robbin and Multi-Channel surfaces structure of the ATP-binding site of the p38α MAP were also developed and displayed with cavity depth, electrostatic potential as well as lipophilic potential to explore the ligand-receptor interactions, furthermore, to examine the 3D contour maps obtained by CoMFA and CoMSIA. Figure 7(a) depicts the MOLCAD cavity depth potential surface of the ATP pocket within compound 9, the cavity depth color ramp ranges from blue (low depth values = outside of the pocket) to light red (high depth values = cavities deep inside the pocket). As shown in Figure 7(a), the phenyl group of R 2 position was in the relative lower depth, while the other parts of the compound 9 were anchored deep inside the ATP pocket. Figure 7(b) demonstrates the MOLCAD electrostatic potential surface of the ATP-binding region, the color ramp for EP ranges from red (most positive) to purple (most negative). In Figure 7(b), the R 3 position was found in a yellow area, which indicated that electron-withdrawing properties were essential for potency; the R 5 position was in a blue area, which suggested that electron-donating properties were crucial for activity. The observations obtained from this electrostatic potential surface were satisfied according to the corresponding CoMFA and CoMSIA electrostatic contour maps. Figure 7(c) shows the MOLCAD lipophilic potential surface of the ATP pocket, the color ramp for LP ranges from brown (highest lipophilic area of the surface) to blue (highest hydrophilic area). The phenyl group at C-5 position was in a brown area, which revealed the importance of hydrophobic properties for this position; the N-15 position was in a white area, which indicated the significance of hydrophilic properties for this position. These observations are in agreement with that obtained from the CoMSIA hydrophobic contour map.

Summary of the Structure-Activity Relationship
The structure-activity relationship revealed by 3D-QSAR and molecular docking studies is illustrated in Figure 8. In detail, the minor substituent in the R 1 position would be favored; the bulky, electron-donating groups in R 3 position would increase the activity; the bulky, electron-withdrawing, hydrophilic groups in R 4 position would benefit the potency; the minor, electron-donating substituent in R 5 position would be favorable; the -F at the C′-2 position and the carbonyl group at C-14 position are crucial for binding to the ATP pocket.

Design for New Inhibitors
Based on the structure-activity relationship revealed by this study, we have designed a series of novel inhibitors, these molecules were aligned to the database and their activities were predicted by the CoMFA and CoMSIA models previously established. The chemical structures and predicted pIC 50 values of these compounds are shown in Table 4. The predicted pIC 50 values of these compounds versus that of the most active compound 9 are illustrated in Figure 9. Most of the designed molecules showed better pIC 50 values than compound 9, which validated the structure-activity relationship obtained in this study.

Conclusion
In conclusion, the 3D-QSAR and docking models established in our present study are quite reliable to efficiently guide further modification of the fused pyrazole derivatives for obtaining better inhibitors. Both the CoMFA and CoMSIA models provided the significant correlations of biological activities with steric, electrostatic, hydrophobic, hydrogen bond donor and acceptor fields. In comparison to CoMSIA, the CoMFA method was found to afford a slightly better predictability. The structure-activity relationships revealed by 3D-QSAR and docking were validated by newly designed derivatives. The results served as a useful guideline for developing novel p38α MAPK inhibitors.