Comparative amino acid decomposition analysis of potent type I p38α inhibitors

Background and purpose of the study p38α is a member of mitogen-activated protein kinases (MAPK) considered as a prominent target in development of anti-inflammatory agents. Any abnormality in the phosphorylation process leads to the different human diseases such as cancer, diabetes and inflammatory diseases. Several small molecule p38α inhibitors have been developed up to now. In this regard, structural elucidation of p38 inhibitors needs to be done enabling us in rational lead development strategies. Methods Various interactions of three potent inhibitors with p38α active site have been evaluated in terms of binding energies and bond lengths via density function theory and MD simulations. Results Our comparative study showed that both ab initio and MD simulation led to the relatively similar results in pharmacophore discrimination of p38α inhibitors. Conclusion The results of the present study may find their usefulness in pharmacophore based modification of p38α inhibitors.


Introduction
In the process of drug discovery, lead generation is a key bottleneck. The costly experimental testing of so many compounds leads to a real challenge in high throughput screening (HTS) method and makes it critical to perform virtual screening techniques to reduce the size of chemical collection richen in active compounds [1]. Computer based prescreening of chemical databases has found its key role in lead identification and is known as in silico drug design. Generally speaking, in silico drug design falls into four categories that are related to each other depending on the structural information on targets and their ligand. These methods are structure-based design, ligand-based design, combinatorial chemistry-based design and de novo design [2].
De novo design techniques are used in the case of known receptor structure and unknown ligand structure. One of the most efficient and rational methods to afford this challenge is fragment based drug design [3,4]. In fragment based drug design, binding of small molecule fragments (MW < 300, No. H-donors < 3, No. H-acceptors < 3, ClogP < 3) [5] to specific domain of active site is evaluated. Based on the binding energies, best fragments are selected and bridged together with appropriate linker(s) to generate new scaffolds. The reverse process, i.e. fragmentation of ligands to constructing fragments, can be used for modification of known ligands. By fragmentation, the chemical diversity of fragment database decreases and the chance of achievement to new lead compound increases. In this process, assessment of interaction between fragments and receptor is the rate limiting step. Estimating the contribution of individual amino acid-ligand interaction energies in total binding energy, i.e. Amino acid Decomposition Analysis (ADA) [6][7][8], would be a very useful trend in fragment development. ADA is based on receptor structure and could be applied to different types of scaffolds. The power of ADA in predicting the effect of individual residues on ligand-receptor interactions can be used as supporting information in drug design. In this regard, estimation of the optimum binding geometry could assist in choosing the best fragment(s) leading to the improved ligand potency profiles.
The phosphorylation of proteins by protein kinases, the largest family of signaling proteins, regulates cell life. More than 500 protein kinases are encoded by the human genome and it is no surprise that any abnormality in the phosphorylation process would lead to the different human diseases such as cancer, diabetes and inflammatory diseases. Different types of these regulating enzymes are introduced as therapeutic target. The active site conservation between protein kinases makes it a real challenge to design selective agents [9][10][11][12]. Therefore evaluation of structural features of these protein kinases and the role of fragments to achieve selectivity may be regarded as an important topic.
Considering the important role in inflammatory pathways, p38 can be regarded as an attractive target to design and develop anti-inflammatory agents. Indeed, p38α is a distinguished target in development of anti-inflammatory agents. Different classes of p38α inhibitors have been developed up to now and their pharmacophore were evaluated in detail [22][23][24][25][26][27]. In the present contribution, we used MD simulations and ab initio method to evaluate pharmacophore model of three potent type Ι p38α inhibitors comprehensively. The results of both MD and ab initio methods were reported and compared with each other. Three different inhibitors, diarylimidazole [28], dihydroquinazolinone [29] and 2-arylpyridazin-3-one [30] scaffolds were selected for our study (Figure 1). These inhibitors are direct ATP-binding site inhibitors with sub-micromolar to nanomolar activity. SB203580 inhibits p38α and β with almost similar potency. This compound is~10 times selective towards p38α/β compared to p38γ/δ [31].
In the case of SB203580, crystallographic studies (PDB code 1A9u, resolution: 2.50 Ǻ) demonstrated that pyridyl nitrogen formed a hydrogen bond with Met109 ( Figure 2). Moreover; 4-fluorophenyl ring occupied the hydrophobic pocket adjacent to the Met109. These two types of interactions have been observed in most of the ATP-binding inhibitors [32]. Nitrogen atom of imidazole ring interacts with Lys53 via hydrogen bond and electrostatic forces. Electrostatic forces are long range interactions between ligand and receptor and have determinant effect on ultimate ligand-receptor complex stability.
In this project we used ab initio method and MD simulations to pharmacophore modeling of these three potent inhibitors via ADA strategy. B3LYP functional [33] and BP86 functional [34] together supplemented with triple-ζ basis set (TZV) [35] were used with the aim of: -Evaluation of ligand-p38α stability via molecular dynamic simulations -ADA through MD simulations and ab initio method to detect hot spots in p38α active site -Comparative evaluation of MD and ab initio results

Materials and methods
Preparation of data set X-ray crystallographic structures of p38α with its cognate ligands were obtained from Brook Haven Protein databank (www.rcsb.org). The Swiss PDB Viewer program (http://www.expasy.org/spdbv) [36] was used to rebuilt and add missed atoms of Lys15 and Arg173 in 1M7Q and 2I0H PDB structures. The crystallographic holo structures were used as a starting point for MD simulations. The force field parameters of ligands were obtained using PRODRG, an automated topology generation tool web server [37]. The electrostatic potential (ESP) derived charges of ligands were all recalculated according to Breneman and Wiberg [38] by using B3LYP/TZV(2d,2p) method and were adjusted in topology file.
The evaluated amino acid residues were selected on the basis of information from schematic 2D representations of ligand-receptor interactions produced by LIGPLOT program [39].
In ab initio studies, all amino acids were considered in their real electrostatic state. Each residue under study was truncated at the C-terminal and N-terminal. N-terminal was acetylated and C-terminal was methyl amidated to mimic the original electron density profile. All conformational and configurational features were the same as the X-ray structure. Molecular images were produced using VMD program [40].

Molecular dynamics simulations
Molecular dynamics simulations were carried out using GROMACS 4 package with the standard GROMOS96 force field [41]. In each case, the p38-ligand complex was solvated in a cubic box with the dimension of 95 × 95 × 95 (all in Ǻ). Explicit simple point charge model (SPC216) was used to represent water molecules [42]. Na + atoms were added to neutralize the total charge of the systems. Short-range interactions were evaluated using a twin-rang cutoff with van der Waals and electrostatic interactions truncated at 12 and 10 Ǻ, respectively. Particle mesh Ewald (PME) method was used to evaluate long-range electrostatic interactions. The protein-ligand complex and waters along with ions were coupled in a temperature bath at 300K, separately (τ T = 0.1 ps). Berendsen barostat (τ p = 1 ps) was used to maintain pressure in 1.0 atm [43,44]. Linear constraint solver (LINCS) was applied to constrain all bonds [45].
In the first step, energy minimization was performed using steepest descent integrator realized in GROMACS package. After energy minimization, 100 ps NVT and NPT ensembles were used to equilibrate system. During NVT and NPT ensembles a harmonic position restrain (1000 kJ mol -1 nm -1 ) was applied to all the heavy atoms of the p38α and ligand. Following equilibration step, production of MD simulation was conducted for 20 ns without any constrains. Finally, analysis programs included in GROMACS package were used to evaluate trajectories.

Ab initio method
All calculations were performed on a structure that obtained by averaging over last 10 ns of MD simulations. Heavy atom-hydrogen bonds were optimized by using heavy atom fixing approximation (constrained optimizations). Constrained optimization was used to prevent irrational movement of the side chains. Otherwise the position of ligand and relevant residue could be changed drastically. Neese and Bykov evaluated optimization errors of this type [46]. According to their evaluation, obtained results are reliable. BP86 functional together in association with triple-ζ basis set (TZV) was used in optimization process. Resolution-of-identity (RI) approximation [47] together with fitting auxiliary basis set TZV/J [48] was applied for all atoms.
Energy calculations were done using B3LYP (Becke-Lee-Yang-Parr) functional in association with triple-ζ basis set (TZV) on optimized structures. For these calculations, chain-of-spheres (RIJCOSX) approximation [49] was invoked. Two set of first polarization functions were applied on hydrogen and non-hydrogen atoms [TZV(2d,2p)]. To consider long-range dielectric effect of protein in our calculation, COSMO model [50] with a dielectric constant of 4.8 was applied [51]. All calculations were done using the ORCA quantum chemistry package [52].
Ligand-residue binding energies (ΔE b ) were calculated using the previously introduced equation [53]. Counterpoise correction [54] was used to take into account basis set superposition error (BSSE). In the case of SB203580, Potential Energy Surface (PES) scan were performed in the direction of hydrogen bond with Met109 in forty steps considering 0.05Å step sizes. The PES calculations were performed by the same method and basis set as mentioned above.

MD simulations
Crystallographic structure of p38α with its cognate ligands enabled us to perform MD simulations and evaluate the role of individual amino acids in total binding energy. This structure was used as starting conformation for our simulations. In the first step, we performed a 20 ns MD simulation to reach a stable trajectory. The stabilities of trajectories were confirmed by assessment of total energy, temperature and RMSD (Figures 3 and 4).
The Average temperature during 20 ns MD simulation at 300.0K was found to be 300K (±1.32, 1.34 and 1.33 respectively) for these systems (Figure 3). These results represented that the obtained equilibriums and energy conservations for the studied systems were desirable.
The RMSDs of ligands in the active site of p38α with respect to their initial structures were used to evaluate the stabilities of these three complexes.
In the case of SB203580 ( Figure 4A), RMSD rose up to 1.13 Ǻ at the beginning of MD simulations and fluctuated around this value (± 0.15) for the rest of simulation. This distribution pattern demonstrated us that ligand achieved to the equilibrium state just after 1 ns distinguished by the RMSD profile. For dihydroquinazolinone scaffold (1M7Q), RMSD increased to 0.71 Ǻ and leveled off to nearly 3 ns ( Figure 4B). At this point of simulation, RMSD rose up again to 1 Ǻ and leveled off (± 0.13). In the case of 2-arylpyridazin-3-one scaffold, RMSD increased to 0.76 Ǻ and fluctuated around (± 0.16) over the course of MD simulation time ( Figure 4C). These data showed that evaluated ligands reached to an equilibrium state after preliminary fluctuations.
During MD simulations the average of 1.2 H-bond(s) could be detected between SB203580 and p38α active site residues ( Figure 5A). The variation of donor-acceptor  distance during MD simulation could be used to evaluate the forming and breaking of H-bonds. Over the whole process, the donor-acceptor distances less than 3.5 Ǻ demonstrated hydrogen bond formation. As can be seen in Figure 5A, pyridine nitrogen-Met109 NH distance remained less than 3 Ǻ for the whole simulation time (green). This fact could convince us in considering a permanent H-bond between these two moieties. But hydrogen bond between imidazole nitrogen and quaternary amine hydrogen of Lys53 was significantly less detectable. This hydrogen bond was formed and broken frequently during MD simulations (Figure 5A blue).
In the case of dihydroquinazolinone scaffold, an average of 1.5 H-bond(s) with p38α active site residues ( Figure 5B) could be detected. Results showed that Hbond between Met109 NH and ligand O18 atoms (green) existed during whole MD simulation period. The distance between His107 O and ligand HN13 atoms remained less than 3.5 Ǻ during MD simulation (blue). According to obtained results, these two hydrogen bonds are permanent during 20 ns MD simulations. The Gly110 NH-ligand O18 distance (red) fluctuated between 5 ns and 20 ns. However the averaged distance remained higher than 0.35 Ǻ for 98.8% of the simulation period (Temporary H-bond, Table 1).
2-arylpyridazin-3-one scaffold formed an average of 1.2 H-bond(s) with Met109 and Gly110 during MD simulation ( Figure 5C). In this case, the distance between Met109 NH and ligand O18 atoms was almost under 0.35 Ǻ in the whole period (green). But the distance between Gly110 NH and ligand O18 atoms was higher than 0.35 Ǻ in 49.5% of simulation period (blue). These outcomes showed that Met109-ligand and Gly110ligand H-bonds were of permanent and temporary types, respectively.
On the basis of results (Table 1), it might be concluded that hydrogen bond between ligand and Met109 is the key structural point in binding to the receptor. This interaction is the common structural feature of all type Ι p38α inhibitors. More detailed analysis of H-bonds between p38α active site residues and evaluated ligands is summarized in Table 1.
After obtaining an equilibrium system, ADA was carried out as follow: participation of each amino acid in total binding energy was computed by evaluation of Lennard-Jones (LJ) and coulombic interaction energies between each amino acid and ligands via performing an additional 1 ns MD simulation in each case. The results of ADA are shown in Figure 6 and Table 2.
Negative energy shifts showed that the residue made favorable contribution to ligand-receptor interactions.
LIGPLOT program [39] was used to detect residues that interact with ligand in each case. Based on the obtained data, same binding pattern to p38α active site could be detected in all the scaffolds. Interaction energies with hinge region residues (Met109, Gly110, and Ala111) are significant and in each case at least, there is one interaction with these amino acids. Residues constructing hydrophobic pocket in the proximity of Met109 were almost involved in interactions with ligand.
In SB203580, Lys53 was found to be the most significant residue in ligand-receptor interactions (ΔE: -5.77 kcal/mol). Nitrogen atom of an imidazole ring participated in H-bond with quaternary amine hydrogen of Lys53. In fact electrostatic forces between these groups made it a favorable interaction. Lys53 had maximum coulombic and LJ interaction energies in these series (−1.64 and −4.13 kcal/mol, respectively). Electrostatic interactions are important forces in primary approach of ligand and receptor to each other. These types of interactions are of long-range type  (Table 2). Due to the reported pharmacophore models of diverse classes of p38 MAPK [22], interactions with Met109 and this hydrophobic pocket are the chemical features designated for type Ι p38α inhibitors. Tyr35 participated in π-π stacking interaction with para-methylsulfinyl phenyl ring of SB203580 (ΔE: -3.59 kcal/mol).
By using ADA we could model the binding mode of three different p38α inhibitors. Obtained results were in good agreement with evaluated pharmacophore models in the literature [22].

Ab initio evaluation
In this part we used ab initio method to evaluate contribution of individual amino acid-ligand interaction energies in total binding energy and compare obtaining results with MD simulations.
The constrained optimization process was done using BP86/TZV method on a structure which was obtained by averaging over last 10 ns MD simulations. All ab initio studies were done on achieved optimized structure.
Various interaction energies between studied p38α inhibitors and selected residues in the active site were obtained independently. The relevant data are shown in Figure 7.
Interactions between imidazole nitrogen and quaternary amine of Lys53 in SB203580 had the most considerable interaction energy (−4.66 kcal/mol) ( Figure 7A). This strong interaction occurred as a result of electrostatic forces between positive nitrogen (cationic Lys53, atomic partial charge: +0.51) and partially negative imidazole N1 atom (atomic partial charge: -0.62). This ionic dipole interaction had determinant participation in total ligandreceptor binding energy.
Another important interaction could be recorded between Met109 and pyridine nitrogen (−2.30 kcal/mol). Interestingly, residues participated in hydrophobic interactions exhibited repulsive interaction with evaluated inhibitor. In the case of Tyr35, the repulsive interaction may be interpreted on the basis of inappropriate orientation of ligand para-methylsulfinyl phenyl ring versus Tyr35 phenyl ring. It should be noted, p38α inhibitors lacking this moiety might not have any significant effect on ligand potency [28].
Asp168 carboxylic moiety interacts via electrostatic forces with quaternary amine in dihydroquinazolinone ligand. This major interaction had prominent binding energy in this series of residues. Hydrogen bond between Met109 backbone NH and ligand O18 atom had binding energy equal to −8.78 kcal/mol. Negative binding energies could be detected between His107 backbone NH and HN18 (−6.34 kcal/mol), Gly110 backbone NH and ligand O18 atom (−4.16 kcal/mol) but like the other ones all hydrophobic interactions had positive contribution in binding energy. Proximity of Lys53 and ligand quaternary amines made this interaction inefficient.
In the case of 2-arylpyridazin-3-one scaffold, Cation-π interaction could be detected between Lys53 and 4flouro-2-methilphenyl moiety. This interaction had maximum binding energy. Hydrogen binding could be Owing to the important role of hydrogen bond with Met109 in type Ι inhibitors, we decided to optimize the geometric position of the involving functional group in the SB203580 ligand. For this purpose, hydrogen bond distance between Met109 backbone hydrogen and pyridyl nitrogen in SB203580 was scanned in the original direction. The maximum interaction energy was found in the 2.25 Å bond length (data was not shown). The optimum distance between pyridine nitrogen and Met109 backbone hydrogen after 20 ns MD simulation was calculated to be 2.14 Å. These findings interestingly showed that ab initio method and MD simulations converged to the same results. Moreover, it was demonstrated that crystallographic structures may not be appropriate starting points for ab initio calculations in all cases.

Comparison of the two methods
MD simulations and ab initio methods were used to calculate the involvement of each amino acid in total binding energy. The results of applied methods were compared to reveal the accuracy and efficiency levels. Our calculations revealed that MD simulations and ab initio based studies led to the similar trends in estimation of amino acid-ligand binding energies. In both methods residues responsible for major interactions in the p38α active site could be recognized with adaptable level of reproducibility (Figures 6 and 7).
For p38α active site, ab initio method resulted in more repulsive hydrophobic and more attractive electrostatic interactions when compared to MD simulations ( Figure 8). This effect seemed to be probably related to the solvent effect and also interactions among adjacent residues. Moreover B3LYP method tended to produce more polarized wave function in electrostatic interactions [55] leading to false positive values.
For instance in p38α, Lys53 interacted with Asp168 and this electrostatic interaction decreased the attractive interaction between Lys53 and SB203580 in MD simulations. But in ab initio study, just the interaction between Lys53 (truncated at N-terminal and C-terminal) and ligand was considered. Similar binding patterns for nearly all residues could be detected while in the case of charge-assisted interactions (Lys53, Asp168), significant deviations were seen ( Figure 8). However, relatively similar binding energies were estimated for Lys53 in SB203580. Two rationales might be envisaged for this different trend: -Solvent effects were attended in MD simulations while ab initio studies were performed in vacuum -In SB203580, Lys53 interacted with partially negative imidazole sp 2 nitrogen atom -Probable repulsive forces between dihydroquinazolinone piperazinium and Lys53 quaternary amine -2-arylpyridazin-3-one scaffold contributed in cation-π interaction with Lys53

Conclusion
We applied totally 60 ns MD simulations and ab initio method to evaluate and compare the accuracy of these methods in predicting pharmacophore models of three different p38α MAPK inhibitors. Both methodologies were able to unravel key interactions with different p38α inhibitors. One advantageous feature of DFTbased calculations is their relatively adaptable outputs regarding significantly shorter processing times due to the incorporated approximations. Results indicated that LJ interactions contributed significantly to binding of SB203580, dihydroquinazolinone and 2-arylpyridazin-3one scaffolds despite the important role of electrostatic interactions in initial approach of ligands to the receptor. We used enzyme structure that was obtained by averaging over last 10 ns of MD simulations for our ab initio studies. This technique led to the results indicating that crystallographic primary structures may be used with care in such calculations. Further investigations would be required to rule out widespread structure activity relationships of p38α inhibitory activity. Finally the results of the present study may find their usefulness in pharmacophore based modification of p38α inhibitors.