Theoretical Characterization of Binding Mode of Organosilicon Inhibitor with p38: Docking, MD Simulation and MM/GBSA Free Energy Approach

Abstract

P38 mitogen activated protein (MAP) kinase is an important anti-inflammatory drug target, which can be activated by responding to various stimuli such as stress and immune response. Based on the conformation of the conserved DFG loop (in or out), binding inhibitors are termed as type-I and II. Type-I inhibitors are ATP competitive, whereas type-II inhibitors bind in DFG-out conformation of allosteric pocket. It remains unclear that how these allosteric inhibitors stabilize the DFG-out conformation and interact. Organosilicon compounds provide unusual opportunity to enhance potency and diversity of drug molecules due to their low toxicity. However, very few examples have been reported to utilize this property. In this regard, we performed docking of an inhibitor (BIRB) and its silicon analog (Si-BIRB) in an allosteric binding pocket of p38. Further, molecular dynamics (MD) simulations were performed to study the dynamic behavior of the simulated complexes. The difference in the biological activity and mechanism of action of the simulated inhibitors could be explained based on the molecular mechanics/generalized Born surface area (MM/GBSA) binding free energy per residue decomposition. MM/GBSA showed that biological activities were related with calculated binding free energy of inhibitors. Analyses of the per-residue decomposed energy indicated that van der Waals and non-polar interactions were predominant in the ligand-protein interactions. Further, crucial residues identified for hydrogen bond, salt bridge and hydrophobic interactions were Tyr35, Lys53, Glu71, Leu74, Leu75, Ile84, Met109, Leu167, Asp168 and Phe169. Our results indicate that stronger hydrophobic interaction of Si-BIRB with the binding site residues could be responsible for its greater binding affinity compared with BIRB.

Introduction

The p38 mitogen-activated protein kinases (p38 MAPK) are serine/threonine specific protein kinases, which participate in a signaling cascade controlling cellular responses.1 They are activated by several stimuli such as cellular stresses or through immune responses. Among the different isoforms (α, β, γ and δ) of p38, α was tested in various animal models of inflammation and has been identified as a crucial antiinflammatory drug target.2-6 It has been observed that most of the kinase inhibitors compete with adenosine triphosphate (ATP), which bind in the active DFG-in conformation, making them type-I inhibitors. Design and development of type-I inhibitors is challenging because of the conserved catalytic domain and same ATP binding site. Therefore, some groups targeted the inactive DFG-out conformations to develop allosteric inhibitors. DFG-out conformation of kinases allows inhibitors to interact with back pocket which is not much conserved among kinases. Inhibitors targeting allosteric binding pocket are termed as type-II inhibitors. Many research groups performed studies on the DFG-out conformation of kinases and developed type-II kinase inhibitors.7-17 Previous reports indicated that allosteric binding pocket of p38 was used to find out potent molecules. This lead to discovery of potent allosteric inhibitor BIRB, currently in phase II clinical trials.6,18

The introduction of bioisostere in medicinal chemistry is a key strategy to improve properties of a molecule and to introduce a new drug-like chemical space into drug discovery process. The incorporation of silicon as a replacement for carbon (C/Si switch) provides an attractive approach for drug design. Sila-substitution (C/Si exchange) of existing drugs enables the search for new drug-like candidates with valuable biological properties. Organosilicon compounds provide opportunities to enhance potency and improve pharmacokinetic and pharmacological properties due to differences in their chemical properties.19,20 Though, some similarities exist between the chemistry of carbon and silicon due to their adjacent position in group IV of the periodic table,21 but the fundamental differences between them can lead to substantial modifications in the physicochemical and biological properties of the silicon-containing analogues.19

Organosilicon agents offer several potential benefits as compared to their relative carbon structures. Higher lipophilicity of organosilicon molecules generally improves cell and tissue permeability and modifies their potency and selectivity. Larger covalent radius of silicon effect the conformation and reactivity of ring structures comprising of a silicon center. Silicon transforms the metabolic process of organosilicon molecules. Also no well-known intrinsic “element-specific” toxicity is related to organosilicon molecules.22

Recently, it has been reported that the C/Si exchange has lead to the identification of potent p38 allosteric inhibitors.23 The present study was performed to understand the binding mechanism of allosteric inhibitors and p38 in DFG-out conformation, by combining molecular docking, dynamics simulation, MM-GBSA free energy calculation and per residue energy decomposition techniques. The differences between the biological activities of BIRB and aryl silane derivative of BIRB (referred as Si-BIRB) were explained on the basis of MM-GBSA binding free energy decomposition. Our results could be exploited for the rational design and development of more potent p38 silicon based isosteric inhibitors.

 

Materials and Methods

Molecular Modeling. Atomic coordinates of BIRB and p38 were retrieved from protein data bank (PDB code: 1KV2, resolution 2.8 Å).24 This structure is a DFG-out conformer and activation loop part is not well defined. Firstly, we modeled the missing loop (Asn115-Leu122) and activation loop (Gly170-Ala184) using MODELLER 9v4 program25 and then refined using loop refinement protocol of MODELLER. Resultant modeled structure was then subjected to biopolymer module of SYBYL 8.1 program.26 Finally this structure was energy minimized using TRIPOS force field for 100 steps with staged minimization protocol in SYBYL. Complexed BIRB was extracted from the relaxed p38 structure and checked for correct atoms and bonds assignment and saved as SYBYL’s mol2 format. Redocking of BIRB into p38 active site was performed using Surflex-Dock in the predefined active site of p38. This step was taken to check performance of docking program. Silicon isostere23 (Si-BIRB) of BIRB was sketched using SYBYL sketching program and minimized. Partial atomic charges were applied using Gasteiger-Hückel method. Figure 1 represents the chemical structures of BIRB and Si-BIRB. The binding pose with the highest score was selected.

Figure 1.Chemical structures of BIRB and Si-BIRB.

MD Simulation. AMBER 12.0 program was used for MD simulations.27 Selected docked conformations of BIRB and Si-BIRB were extracted from p38 active site. These conformations were fed to Antechamber program28 to calculate semi-empirical charges29 using Austin model 1 (AM1) Hamiltonian.30 Force field parameters for ligands were generated using general AMBER force field (GAFF),31 and macromolecules AMBER ff99SB.32 The ligand-protein complexes were placed in a TIP3P33 water box with 10 Å distance along XYZ directions. This system was neutralized by adding appropriate positive counter ion (1 Na+).

Ligand-protein complexes (BIRB-p38; system 1 and Si-BIRB-p38; system 2) were subjected to minimization to remove steric clashes. 2000 steps of total minimization was performed with the initial 500 steps steepest descent followed by remaining steps of conjugate gradient minimization. Time interval was set to 2 femto seconds ( fs). Weak harmonic restraints (2.0 kcal/mol × Å2) were applied to the all atoms of ligand and protein. All the simulations were performed by applying periodic boundary conditions. After minimization, a canonical NVT simulation was performed with the gradual temperature increment from the 0 to 300 K in 50 picoseconds (ps). Langevin thermostat was used to couple temperature with a force constant of 1.0 kcal/mol·Å2 and coupling coefficient of 1.0/ps. Followed by NVT simulation, 50 ps of density equilibration was performed by keeping weak harmonic restraints (2.0 kcal/mol·Å2) on ligand-protein. Constant temperature (300 K) and pressure (1 atm) was maintained during simulation, and pressure relaxation time of 1.0 ps was used. Finally, equilibration was performed at the constant pressure (1 atm) and temperature (300 K) for the 500 ps, by removing all the restraints on ligand-protein. Here, pressure relaxation time was 2.0 ps. After system stabilization, production run of 15 ns was performed. Long range electrostatic interactions were handled using Particle Mesh Ewald (PME) method.34 The SHAKE algorithm was used to constrain all hydrogen atoms.35 Non-bonded cut off was set to 8.0 Å.

MM/GBSA Binding Free Energy Calculations. Kollman et al.36 introduced MM/GBSA method to compute binding free energy of ligand-protein system. Successful applications of MM/PBSA and MM/GBSA methods in binding free energy calculations are documented in literature.37-41 For binding free energy estimation, 500 snapshots were taken from the last 5 ns trajectory. Non-polar solvation free energy is computed using the solvent accessible surface area (SASA) in AMBER 12 using MOLSURF42 module, and polar contribution is calculated using GB model (igb = 2).43

ΔGnonpolar,solv = γSASA + b

Where γ, represents the surface tension and b is a constant. Parameters for γ and b were set to 0.0072 kcal/mol·Å2 and 0, respectively. Probe radius to calculate SASA was set to 1.4 Å. Dielectric constant value for protein ligand complex was set to 1, while external solvent 80.0.

Computationally expensive and time consuming normal mode analysis was performed to estimate solute entropy contribution (−TΔS) upon ligand binding by using NMODE module of AMBER 12. To compute entropy contribution, fifty frames were selected at specific interval from the last 5 ns simulation trajectory. Before normal mode analysis, each of the complex, receptor and ligand systems were minimized within 50,000 steps. The entropy term was decomposed into rotational (ΔSrotat), translational (ΔStransl) and vibrational (ΔSvibrat) contributions.

Free Energy Decomposition. Residues in the 4 Å radius of ligand were selected to compute their contribution towards binding free energy. To compute binding free energy between ligand-residues pairs, fifty frames were selected from the last 5 ns trajectory. Frames were selected at particular interval. MM-GBSA method was used to calculate binding free energy contribution. It allows calculating the free energy of interactions between each ligand and residue pair, as well as contribution of backbone and side chain. MM-GBSA calculated the electrostatic, van der Waals, polar solvation and non-polar solvation energy for each pair of ligand-residue. Sander program was used to calculate the contributions of ΔGele and ΔGvdW, whereas, MM-GBSA decomposition method was used to calculate polar solvation contribution between each ligand and residue pair. SASA program was used to calculate the non-polar contribution of desolvation.

 

Results and Discussion

Molecular Modeling and Docking Analysis. Two missing loops were built using MODELLER program, which are shown in supplementary material S1. This modeled structure was used for the identification of binding modes of BIRB and Si-BIRB using Surflex-Dock module of SYBYL. Co-crystal ligand site was selected as binding site for docking experiment as shown in Figure 2. The binding pose of co-crystal ligand (BIRB) was reproduced using Surflex-Dock to ensure the validity of docking protocol. Figure 3(a) represents the docked mode of BIRB in p38 active site, and Figure 3(b) gives closer look of interactions pattern. Similarly, Figure 4(a) shows binding mode of Si-BIRB into the p38 active site and Figure 4(b) represents closer look of ligand-p38 interactions. It has been found that nitrogen of urea group of both BIRB and Si-BIRB formed hydrogen bonds with the Glu71. However, the tri-methyl-silane of Si-BIRB was docked deep inside a hydrophobic cavity formed by the residues Leu74, Leu75, Met78, Val83, Ile84, Ile166 and Leu167. Para-methoxy phenyl of Si-BIRB was docked in between Arg67, Arg70 and Glu71. Biphenyl moiety of the Si-BIRB and BIRB was accommodated into a pocket formed by the hydrophobic residues Val38, Ala51, Lys53, Ile84, Leu86, Leu104, Val105, Thr106 and Phe169. Ethoxymorpholino group of BIRB and Si-BIRB was docked into the Val38, Leu108, Met109, Gly110 and Asp169. Both BIRB and Si-BIRB interacted through the hydrophobic interactions and partly through hydrophilic interactions. Dynamic interactions between the ligand-receptor were investigated using MD simulations.

Figure 2.Binding pocket residues were shown by the magenta stick, whereas p38 is shown by transparent cyan cartoon.

Figure 3.(a) Stereo-view of superposed re-docked structure of BIRB over co-crystal (BIRB). Co-crystal BIRB is shown by cyan stick, whereas docked BIRB by magenta, p38 by cyan cartoon. (b) Stereo-view of docked mode of BIRB inside p38 active site. BIRB shown by the magenta color and active site residues by the green stick. Hydrogen bonds between Glu71 and urea NH of BIRB were represented by the cyan dashed line. This figure was generated using Pymol program (www.pymol.org).

Figure 4.Stereo-view of docked mode of Si-BIRB inside p38. (a) Docked mode of Si-BIRB superposed over co-crystal (BIRB). Co-crystal BIRB is shown by yellow stick, whereas docked Si-BIRB represented by magenta stick and p38 by brown cartoon. (b) Closer view of docked mode of Si-BIRB inside p38 active site. Si-BIRB was shown by cyan stick and active site residues by green stick. Hydrogen bonds between Glu71 and urea NH of Si-BIRB were represented by yellow dashed line.

Root Mean Square Deviation (RMSD). Docked poses of BIRB and Si-BIRB along with receptor (p38) were used in AMBER for MD simulations. Both the systems were simulated for 15 ns each in an explicitly solvated water box. Protein structure stability was measured as backbone atoms (C, Cα and N) deviation from its starting structure. All atom RMSD was calculated for the BIRB and Si-BIRB throughout the simulation period. Figure 5(a) shows RMSD plot for p38 in two simulated systems, whereas Figure 5(b) represents all atom RMSD for BIRB and Si-BIRB. RMSD plot for p38 in system 1 showed that protein structure took long time to equilibrate (~6000 ps). Continuous rise in RMSD was observed upto equilibration period and later on it formed a plateau for the rest of simulation period. After equilibration average RMSD was 2.65 Å. However, calculated RMSD for p38 in system 2 indicated that protein equilibrated upto ~1500 ps. Later RMSD decreased upto ~4000 ps and elevated again upto ~6000 ps and later formed plateau till the end of the simulation. Average RMSD of 2.81 Å was observed after equilibration period.

Figure 5.(a) RMSDs of the backbone atoms of p38 in system 1 (black) and 2 (red) were plotted as a function of time. (b) RMSDs of BIRB (black) and Si-BIRB (red) as a function of time.

RMSDs for the inhibitors are shown in Figure 5(b). For BIRB, RMSD continuously rose upto ~2600 ps and later fell back throughout the simulation period. However, calculated RMSD for Si-BIRB showed that less deviation were occurred during simulation. The average RMSD for Si-BIRB was 1.45 Å. RMSDs for both the inhibitors were reduced during last 4000 ps. Average temperature and pressure of both the simulated systems were calculated and plotted in supplementary material S2. From the figure, it is evident that the simulations were carried out at constant temperature and pressure.

Potential Energy and Root Mean Square Fluctuation (RMSF). Potential energies of the simulated systems were calculated and shown in Figure 6(a). It shows that the potential energy for system 1 fluctuated from −143227 kcal/ mol to −142110 kcal/mol. However, potential energy for the system 2 fluctuated from −140666 kcal/mol to −139345 kcal/mol. Average potential energies for system 1 and 2 were −142659 kcal/mol and −140078 kcal/mol, respectively. Potential energy plots indicated that systems were stable throughout simulation and formed a plateau.

Figure 6.(a) Potential energy plot of system 1 (black) and 2 (red) as a function of time. (b) RMSF of backbone atoms of p38 during MD simulation. RMSF for p38 in system 1 and 2 are represented by black and red lines.

RMSFs of p38 in two simulated systems were calculated and graphically plotted in Figure 6(b). RMSF plot indicates that very few residues were fluctuated with the intensity > 4 Å. Higher RMSF values were observed for the residues which belong to the loop regions of p38. RMSF profile for p38 simulations suggested similar trends of dynamic feature for both the structures. This observation further suggested similar interaction mechanism for the simulated inhibitors. Active site residues of both the simulated systems showed fluctuations < 3 Å, which indicated that active site regions of p38 were stabilized by ligands.

Binding Mode Analysis After MD Simulation.

Binding Mode of BIRB: To further investigate the effect of MD simulation on binding mode of BIRB inside p38, average structure of the last 3 ns MD simulation was generated for BIRB-p38 that is shown in Figure 7. It can be seen that there is not much deviation after MD simulation, which is evident through the overlapped structures. Closer view of interaction between BIRB and p38 shows that ligand interacted through hydrophilic and hydrophobic interactions. Significant hydrogen bond interactions occurred with important residues inside the binding pocket. It can be seen that the hydrogen bond plays a crucial role in the binding of BIRB to p38. Four hydrogen bonds were observed and the distance profile of hydrogen bonds is given as a function of time in Figure 8(a). It was observed that the hydrogen bond between urea ‘O’ of BIRB and main chain NH of Asp168 was most important and maintained the distance of ~2 Å throughout simulation period. However hydrogen bond between main chain NH of Met109 and morpholino ‘O’ of BIRB identified by docking disappeared after ~1 ns MD simulation, and then reappeared after ~11 ns MD simulation and later maintained the distance of ~2 Å. This observation also indicated the role of hinge region residue to stabilize BIRB inside p38. Glu71 also interacted with the urea NH of BIRB through hydrogen bond interactions. Figure 8(a) shows the distance between OE2 of Glu71 and NH of BIRB urea moiety. It indicates that during MD simulations both the hydrogen bonds were unstable upto ~6 ns, but later on found to be stable throughout simulation with ~2 Å distance upto ~14 ns time. 2D-atomic interactions between BIRB and p38 were generated using Ligplot server and shown in supplementary material S3.

Figure 7.Stereo-view of binding mode of BIRB after MD simulation. (a) Superposed binding mode of BIRB before (white stick) and after MD simulation (cyan stick) inside p38. (b) Closer view of interactions between BIRB and p38 after MD simulation. BIRB is represented by grey stick, active site residues by magenta stick and p38 by yellow cartoon. Hydrogen bonds between BIRB and p38 were depicted by the dashed yellow.

Figure 8.Calculated hydrogen bond distances between (a) BIRB (b) Si-BIRB and p38 during MD simulations. Color index is shown at the right side. BIRB and Si-BIRB are numbered as residue number 349.

Binding Mode of Si-BIRB: Figure 9(a) represents the overlaid structures of complex before and after MD simulation and Figure 9(b) indicates closer interactions after MD simulation. After MD simulation ligand showed similar orientation and was stabilized inside p38 by hydrophilic and hydrophobic contacts. It was observed that OE2 of Glu71 formed two hydrogen bonds with the urea NH. During initial simulation period (~1.5 ns) both hydrogen bonds were fluctuating, but later on found to be stabilized with the ~2 Å distance. This indicated that hydrogen bonds held the ligand tightly inside active site. Similar to BIRB, Si-BIRB formed hydrogen bond with the main chain NH of Asp168 (Fig. 8(b)). This hydrogen bond seemed to be important to hold ligand tightly, because throughout simulation period it showed ~2 Å distance. Similar trend of hydrogen bond as of BIRB was observed between morpholino ring and NH of Met109. During simulation this hydrogen bond was destabilized and later on after ~12 ns it reappeared with the average distance of 2 Å for rest of simulation period. At this region the binding pocket was opened and there was more room for substituent to move. This is the reason, why this hydrogen bond was unstable for most of the simulation period in both the system (Fig. 8). 2D-atomic interactions between Si-BIRB and p38 were calculated using Ligplot server and shown in supplementary material S4. Indeed, hydrogen bonds are well known to break and reform dynamically, and only a little bit of that behavior is observed on this time scale. Morpholine group showed more conformational flexibility than the urea group when not held in place by a hydrogen bond. This can be seen through the fluctuating distance profile between morpholino oxygen and Met109 (Fig. 8). However, urea oxygen is bound firmly to the Asp168 backbone and not allowed it to move during a short 15ns (Fig. 8).

Figure 9.Stereo-view of binding mode of Si-BIRB after MD simulation. (a) Superposed binding mode of Si-BIRB before (sky blue stick) and after MD simulations (yellow stick) inside p38. (b) Closer view of interactions between Si-BIRB and p38 after MD simulation. Si- BIRB is represented by yellow, active site residues by magenta stick and p38 by light brown cartoon. Hydrogen bonds between Si-BIRB and p38 were showed by the dashed yellow.

It was also observed that the salt bridge between Lys53 and Glu71 was maintained throughout simulation period (Fig. 10). In both the simulated system it was observed that average salt bridge distance was 3 Å, while maintaining important hydrogen bonds with the simulated ligands.

Figure 10.Calculated salt bridge distances between Lys53 and Glu71 during MD simulations in (a) system 1 and (b) 2. Color index for the salt bridges were shown at the right side of figure.

MM/GBSA Binding Free Energy Based Interaction Mechanism. To better understand in depth interaction mechanism of simulated inhibitors and p38, we performed MM/ GBSA binding free energy calculations to characterize their binding affinity. It has been previously reported that the MM/GBSA perform well than the MM/PBSA in predicting relative binding free energies.39 Binding free energy approaches like MM/PBSA and MM/GBSA are computationally more efficient36,44,45 than the theoretically rigorous approaches such as thermodynamic integration46 and free energy perturbation.47-49 In this study, the binding free energies were estimated using mmpbsa program in AMBER. Previous report indicates that the binding free energies obtained from MM/PBSA and MM/GBSA methods correlate well with the experimental values, but it cannot reproduce absolute experimental values accurately.24,50 MM/GBSA predicted binding free energies for BIRB and Si-BIRB correlates well with the calculated binding free energies. However, advantage of MM/PBSA and MM/GBSA method is that it can decompose the obtained binding free energies into its individual components such as van der Waal and electrostatic contribution from gas phase and solvent phase. Thereby, it helps us to understand the complex binding process in detail. Table 1 shows the calculated binding free energy for BIRB and Si- BIRB using MM/GBSA method. From Table 1, it can be seen that the unfavorable contribution occurs for electrostatics in solution phase whereas, favorable contribution occurs in gas phase for both the simulated systems. The resulting balance of electrostatic interaction (ΔGele + ΔGele,solv) in both the simulated systems were shifted towards unfavorable side for the binding process. However, favorable contributions occur from the van der Waals interaction and non-polar solvation energy for both the system. The reason to obtain favorable contribution might come from deep burial of hydrophobic group of inhibitors inside hydrophobic crevices of p38.

Table 1.All the units are in kcal/mol. Standard error of mean is shown in parentheses. ΔGvdw is van der Waals contributions from MM. ΔGele is the electrostatic interactions calculated by MM. ΔGnonpol,solv is the non-polar contribution to solvation. ΔGele,solv is the polar contribution to solvation. ΔGGas is total gas phase energy. ΔGSolv is total solvation energy. TΔS is enthalpy change. ΔGpred is calculated binding free energy by MM/GBSA method. ΔGexpt is experimental binding free energy calculated from IC50 values according to ΔG ≈ RT ln IC50

Binding Modes Comparison of BIRB and Si-BIRB. After MD simulation, binding modes of BIRB and Si-BIRB were superposed in p38 and shown in Figure 11. From this figure we could say that both the inhibitors share same space in p38, but para-methoxy-phenyl of Si-BIRB penetrated deep inside the binding pocket than the tolyl of BIRB, respectively. It was observed that the AMBER calculated binding free energy for Si-BIRB (−53.47 kcal/mol) is higher than the BIRB (−49.62 kcal/mol). Our estimated binding free energies correlates with the biological activity values of Si-BIRB (IC50 = 3 nM) and BIRB (IC50 = 9 nM), it also indicates the usage of MM/GBSA to predict relative binding free energies. However, experimental binding free energies for Si-BIRB and BIRB are calculated from IC50 values using following relationship,

Where, R is ideal gas constant, T is temperature in K (298 K is used in this study) and Cenzyme is the concentration of enzyme, which is a very small number after equilibration and can be omitted in most cases.

Figure 11.Overlay of binding modes of BIRB (magenta) and Si-BIRB (yellow) after MD simulation.

Above approximation was used in previous studies to estimate experimental binding free energies from IC50 values.51,52 Experimental binding free energy for Si-BIRB (−11.63 kcal/mol) and BIRB (−10.98 kcal/mol) correlates well with the AMBER calculated binding free energy. From the experimental binding free energy values it is obvious that the energy difference is very less and it might occur from the variable part of Si-BIRB and BIRB. From the decomposed binding free energy values, it was observed that the major differences in activity values were occurred from the residues surrounding 5-tert-butyl-2-(4-methylphenyl)pyrazol-3-yl and 4'-methoxy-4-(trimethylsilyl)-[1,1'-biphyenyl]-2-yl of BIRB and Si-BIRB, respectively. However, energy differences were also observed for the residues surrounding morpholino ring. More van der Waals interactions energy were noted for the 4'-methoxy-4-(trimethylsilyl)-[1,1'-biphyenyl]-2-yl of Si-BIRB than the 5-tert-butyl-2-(4-methylphenyl)-pyrazol-3-yl of BIRB, because former substituent is more hydrophobic than the later, and this could be the reason to have more favorable free energy for Si-BIRB. Moreover, the pockets surrounding this substituent are formed by the hydrophobic residues. It was also observed that the BIRB (−28.83 kcal/mol) complex is entropically more favorable than the Si-BIRB (−26.21 kcal/mol) complex. Major difference in the entropy value was arise from vibrational entropy for BIRB (−3.9 kcal/mol) and Si-BIRB (−1.04 kcal/mol), respectively. Translational entropy was found to be more favorable in both the complexes.

Identification of the Key Residues Based on Free Energy Decomposition. Crucial residues into the binding process of ligand-protein have been identified from the per residue decomposition energy. Figure 12 represents the contribution of residues in terms of van der Waals (ΔGvdw) and non-polar (ΔGnonpol,sol) energy. It has been observed that the interaction spectra of BIRB and Si-BIRB with p38 are quite similar. Overall, favorable van der Waals and non-polar contributions were occurred from the residues Glu71, Leu75, Ile84, Met109, Leu167 and Asp168 for both the complexes. Residues having relatively larger differences in the binding free energy were carefully scrutinized. Residues with the free energy difference of ~0.2 kcal/mol are Tyr35, Ala51, Lys53, Arg67, Arg70 Thr106, Leu108, Met109, Gly110, Ala111 and His148. From the MD simulated mode of BIRB (Fig. 7(b)) and Si-BIRB (Fig. 8(b)), it is evident that the Tyr35 interacts more closely with the para-methoxy-phenyl moiety of Si-BIRB than the BIRB. This is the reason why Tyr35 has more van der Waals and non-polar contribution in Si-BIRB ~(−2.39 kcal/mol) than BIRB ~(−0.08 kcal/mol). It has also observed that Arg67 is the only residue which has unfavorable van der Waals and non-polar contribution in BIRB (~0.24 kcal/mol) system compare to Si-BIRB ~(−1.8 kcal/mol) system. It is evident from the Figure 7(b) that Arg67 is far from the para-methyl-phenyl moiety. Met109 also shows difference in the interaction energy with the BIRB ~(−0.6 kcal/mol) and Si-BIRB ~(−2.4 kcal/mol). It could be seen from the Figures 7(b) and 8(b) that in case of BIRB the side chain of Met109 facing towards morpholino ring, which is away in the Si-BIRB, this could be the reason to occur free energy difference.

Figure 12.Non-polar and van der Waals contributions of the interacting residues with Si-BIRB (light cyan) and BIRB (grey) are represented as histogram.

The sum of electrostatics (ΔGele) and polar solvation (ΔGpol,solv) energy for the residues nearby BIRB and Si-BIRB are plotted in Figure 13. It has been observed that the electrostatic and polar solvation free energies are unfavorable for binding process. Maximum unfavorable polar solvation free energy was obtained for the hydrogen bond forming residues such as Glu71, Met109 and Asp168.

Figure 13.Histogram of electrostatics and polar contributions of the interacting residues with Si-BIRB (light cyan) and BIRB (grey).

Decomposition of Free Energy into Backbone and Side Chain Contribution. More detailed binding process of ligands to p38 was studied by decomposing per residue energy into backbone (Fig. 14) and side chain (Fig. 15) contributions. Estimated van der Waals and electrostatics contribution of backbone atoms (Fig. 14) indicates that all the residues had unfavorable contribution to the total free energy. However, all the studied residues showed favorable backbone van der Waals energy for BIRB and Si-BIRB. In most of the cases residue showed higher van der Waals contribution for Si-BIRB, than the BIRB. Although, residues like Val38, Ala51, Lys53, Arg70, Leu74, Leu104, Val105, Thr106, Leu108 and Met109 had higher energy contribution for BIRB. Highest backbone interaction energy for BIRB and Si-BIRB was summed up from the Glu71 (−0.568 kcal/mol vs −0.86 kcal/mol), Leu167 (−1.15 kcal/mol vs −1.19 kcal/mol) and Asp168 (−1.51 kcal/mol vs −1.57 kcal/mol).

Figure 14.Electrostatics and van der Waals contribution of backbone atoms of p38 with Si-BIRB and BIRB were represented by the green/violet and light blue/red cylinders, respectively.

Figure 15.Binding free energy contribution (van der Waals) of side chains of interacting residues with Si-BIRB and BIRB are represented by the light blue and red cylinders, whereas electrostatics contribution of side chains with Si-BIRB and BIRB are denoted by the green and violet cylinders, respectively.

Side chain contribution of the studied residues towards total free energy is given in Figure 15. From the graph it can be seen that all the residues have positive interaction electrostatic energy, indicating unfavorable contribution. However, favorable van der Waals contribution to free energy was noted for all the residues. Residues such as Tyr35, Val38, Lys53, Glu71, Leu74, Leu75, Ile84, Leu167 and Asp168 have van der Waals contributions more than −1.0 kcal/mol. It has been observed that the side chain of Tyr35 interacts more hydrophobically with para-methoxy-phenyl moiety of Si-BIRB (−1.79 kcal/mol) than the para-methyl-phenyl moiety of BIRB (−0.26 kcal/mol). Lys53 shows more side chain contribution for BIRB (−1. 9 kcal/mol) and Si-BIRB (−2.03 kcal/mol) in terms of the total binding free energy. Glu71 side chain contributed more energy for BIRB (−1.69 kcal/mol) than Si-BIRB (−1.57 kcal/mol). Highest side chain van der Waals contribution was occurred from the Ile84 for Si-BIRB (−2.15 kcal/mol) than BIRB (−1.94 kcal/mol). It is visible from the Figures 7(b) and 9(b) that it interacts with the naphthalene of both the inhibitors. However, it has been observed that the Leu74 and Leu75 has maximum van der Waals interaction and resulted in the free binding energy. It can be seen that these residues interact with the tri-methyl-silyl-phenyl and 3-t-butyl-pyrazole moiety, of Si-BIRB and BIRB, respectively. We can expect more van der Waals interaction of Leu74 and Leu75 with tri-methyl-silyl-phenyl than the 3-t-butyl pyrazole, because of more hydrophobicity of former than the later. This could be the reason for having more binding free energy of Si-BIRB than the BIRB.

Although the activation loop residues (Asp168 and Phe169) also interacted variably to the simulated ligands. It was observed from Figure 15 that the Asp168 has more side chain contribution to the total binding free energy in Si-BIRB than BIRB. MD simulated binding modes (Figs. 7(b) and 9(b)) shows that the side chain of Asp168 interacts closely with the pyrazole and phenyl moiety of BIRB and Si-BIRB, respectively. From this observation we could say that acetic acid moiety of Asp168 interacts with the phenyl of Si-BIRB in CH-π fashion, which is somehow impossible with the pyrazole of BIRB. This is the reason why side chain of Asp168 has more contribution for binding free energy in Si-BIRB (−1.66 kcal/mol) than the BIRB (−1.33 kcal/mol). It has been already reported that this type (CH-π) of weak intermolecular forces are orientation dependent and contribute 1.5-2.5 kcal/mol additively in enthalpy.53 Hinge regions Leu108 has more side chain contribution towards free energy of BIRB (−0.54 kcal/mol) than Si-BIRB (−0.08 kcal/ mol). From Figure 7(b) it can be seen that the side chain of Leu108 facing towards morpholino of BIRB, which is facing away from the morpholino of Si-BIRB (Fig. 9(b)), this might be the reason to have their erratic contribution towards binding free energy. The simulation time in this study is bit short compared to current standard, but we got equilibrated structures with reasonable RMSDs which were used to perform the analysis.

Our group is also involved in the modeling studies of different protein targets using homology modeling, docking and MD simulation techniques.54-62

 

Conclusion

In this study, we performed docking, MD simulation and MM/GBSA binding free energy calculation for the allosteric inhibitors of p38. Docking guided MD simulation inside p38 revealed that the BIRB and Si-BIRB share the same binding site in p38. However, salt bridge contact between Glu71 and Lys53 were found to be important to stabilize the side chains conformations and to interact with the simulated inhibitors. Hydrogen bonds between simulated inhibitors with two amino acid residues (one with Glu71 of C-helix and the other with Asp168 of activation loop) were found to be important to hold ligand inside p38 and to stabilize DFG-out conformation. Hinge region’s hydrogen bond (Met109) might be responsible for the observed morpholino ring conformations in both inhibitors. MM/GBSA estimated binding free energies showed trend with the biological activity (experimental binding free energy). At the same time, estimated binding free energies indicated that the affinity might come from van der Waals and non-polar interactions. That is, higher hydrophobic interactions between inhibitors and active site residues could be responsible for the inhibitory activity of Si-BIRB and BIRB, respectively. Crucial residues identified for hydrogen bond, salt bridge and hydrophobic interactions were Tyr35, Lys53, Glu71, Leu74, Leu75, Ile84, Met109, Leu167, Asp168, and Phe169. Information obtained from the 3D atomic interaction and decomposition energy could be used to design and develops novel allosteric type-II inhibitors targeting DFG-out conformation of p38.