Insights into the Interaction Mechanisms of the Proviral Integration Site of Moloney Murine Leukemia Virus (Pim) Kinases with Pan-Pim Inhibitors PIM447 and AZD1208: A Molecular Dynamics Simulation and MM/GBSA Calculation Study

Based on the up-regulation of the proviral integration site of the Moloney murine leukemia virus (Pim) kinase family (Pim1, 2, and 3) observed in several types of leukemias and lymphomas, the development of pan-Pim inhibitors is an attractive therapeutic strategy. While only PIM447 and AZD1208 have entered the clinical stages. To elucidate the interaction mechanisms of three Pim kinases with PIM447 and AZD1208, six Pim/ligand systems were studied by homology modeling, molecular docking, molecular dynamics (MD) simulation and molecular mechanics/generalized Born surface area (MM/GBSA) binding free energy calculation. The residues of the top group (Leu44, Val52, Ala65, Lys67, and Leu120 in Pim1) dominated the pan-Pim inhibitors binding to Pim kinases. The residues of the bottom group (Gln127, Asp128, and Leu174 in Pim1) were crucial for Pims/PIM447 systems, while the contributions of these residues were decreased sharply for Pims/AZD1208 systems. It is likely that the more potent pan-Pim inhibitors should be bound strongly to the top and bottom groups. The residues of the left, right and loop groups were located in the loop regions of the binding pocket, however, the flexibility of these regions triggered the protein interacting with diverse pan-Pim inhibitors efficiently. We hope this work can provide valuable information for the design of novel pan-Pim inhibitors in the future.


Introduction
The proviral integration site of the Moloney murine leukemia virus (Pim) oncogene family plays an essential role in tumor growth and development which consists of three constitutively active serine/threonine kinases (Pim1, 2, and 3) [1]. The PIM1 gene has been found to function as a proto-oncogene. When it is overexpressed, it induces lymphomas in transgenic mice, though at a low frequency. However, when Pim1 is co-expressed in mice with the oncogene c-Myc, 100% of the transgenic mice die of lymphomas in utero [2]. Interestingly, the knockout of Pim1 in mice is not lethal nor does its absence induce any immediately obvious phenotype. The absence of Pim1 may be compensated for by Pim2, although not in all cases. In one case where the Pim2 kinase did compensate, it appeared to contribute to cell survival, indicating that in some instances it functions a molpdf: MODELLER minimizes the objective function with respect to Cartesian coordinates of~10,000 atoms (3D points) that form a system. b DOPE score: assesses the quality of the selected atoms in the model using the discrete optimized protein energy (DOPE) method and is a statistical potential optimized for model assessment. c GA341 score: assesses the quality of the model using the GA341 method.
From the protein structure analysis (ProSA) plot ( Figure 1D), the Z-score of Pim3.B99990003 model was −7.18, indicating a high quality model was obtained. A Ramachandran plot ( Figure 1E) validated the Pim3.B99990003 model with 80.7% of the total residues in the most favored region, 17.2% in the additional allowed region and 1.5% in the generously allowed region (Leu21, Arg169, Asp297, and Val300). Only two residues (Ser35 and Thr310) were located in the disallowed region, which constituted 0.7% of the total protein, which represented a good quality of the predicted model. Moreover, the Pim3.B99990003 model showed 92.9825% of quality factor using error atom type (ERRAT) ( Figure 1F). Thus, the Pim3.B99990003 model was used for the following molecular docking and MD simulations. optimized for model assessment. c GA341 score: assesses the quality of the model using the GA341 method.

Overall Structure and Dynamics
To explore the dynamic stability of these six systems, the root mean square deviation (RMSD) values of the whole protein backbone atoms were calculated (Figure 2A). The RMSD plot indicated that all six systems achieved equilibrium in a short time. The RMSD values of all systems except for the Pim3/PIM447 system fluctuated around 3 Å. The lower RMSD values implied that the binding patterns of Pim/ligand complexes derived from the molecular docking studies changed little during 20 ns MD simulations. To further assess the stability of six systems after 10 ns, the same RMSD calculations were conducted on the whole protein backbone atoms with respect to the structure at t = 10 ns ( Figure 2B). It was clear that all six systems already equilibrated after 10 ns and fluctuated around 1.5 Å, which implied that 20 ns MD simulations were enough to relax the Pim/ligand complexes. It is noteworthy that the RMSD value of Pim2/AZD1208 was more or less larger than other systems, which may be attributed to the bigger fluctuations of residues Phe49 and Val52 obtained by a root mean square fluctuation (RMSF) calculation ( Figure S1). Therefore, it was reasonable to perform the following binding free energy calculation and free energy decomposition analysis based on the last 10 ns trajectories.

Overall Structure and Dynamics
To explore the dynamic stability of these six systems, the root mean square deviation (RMSD) values of the whole protein backbone atoms were calculated ( Figure 2A). The RMSD plot indicated that all six systems achieved equilibrium in a short time. The RMSD values of all systems except for the Pim3/PIM447 system fluctuated around 3 Å. The lower RMSD values implied that the binding patterns of Pim/ligand complexes derived from the molecular docking studies changed little during 20 ns MD simulations. To further assess the stability of six systems after 10 ns, the same RMSD calculations were conducted on the whole protein backbone atoms with respect to the structure at t = 10 ns ( Figure 2B). It was clear that all six systems already equilibrated after 10 ns and fluctuated around 1.5 Å, which implied that 20 ns MD simulations were enough to relax the Pim/ligand complexes. It is noteworthy that the RMSD value of Pim2/AZD1208 was more or less larger than other systems, which may be attributed to the bigger fluctuations of residues Phe49 and Val52 obtained by a root mean square fluctuation (RMSF) calculation ( Figure S1). Therefore, it was reasonable to perform the following binding free energy calculation and free energy decomposition analysis based on the last 10 ns trajectories.

Binding Free Energies Estimated by MM/GBSA
The calculated MM/GBSA binding free energies (ΔGbind) of Pim1/AZD1208, Pim1/PIM447, Pim2/AZD1208, Pim2/PIM447, Pim3/AZD1208, and Pim3/PIM447) were −31. 17 Table 2). The correlation between the calculated binding free energies and the experimental pKi (i.e,. -logKi ) was high (R 2 = 0.795). Meanwhile, the MM/GBSA calculation also performed well in investigating the binding affinities of three Pim kinases to GDC-0039 and hispidulin (Table S1). Table 2 also showed that the van der Waals interactions played a dominant role not only in the total binding free energies of the six systems, but in ranking the origin of difference in the binding affinities of the systems. The net electrostatic free energy (ΔEele + ΔGGB) of each system had a negative contribution on the binding of PIM447 (or AZD1208) to Pim kinases. It may be attributed to the fact electrostatic interaction is fully counteracted by the solvation effect. It is noteworthy that the differences in entropy contribution (TΔS) of the six systems were not obvious, which implied the normal-mode analysis might be not

Binding Free Energies Estimated by MM/GBSA
The calculated MM/GBSA binding free energies (∆G bind ) of Pim1/AZD1208, Pim1/PIM447, Pim2/AZD1208, Pim2/PIM447, Pim3/AZD1208, and Pim3/PIM447) were −31. 17 Table 2). The correlation between the calculated binding free energies and the experimental pKi (i.e., −logKi) was high (R 2 = 0.795). Meanwhile, the MM/GBSA calculation also performed well in investigating the binding affinities of three Pim kinases to GDC-0039 and hispidulin (Table S1). Table 2 also showed that the van der Waals interactions played a dominant role not only in the total binding free energies of the six systems, but in ranking the origin of difference in the binding affinities of the systems. The net electrostatic free energy (∆E ele + ∆G GB ) of each system had a negative contribution on the binding of PIM447 (or AZD1208) to Pim kinases. It may be attributed to the fact electrostatic interaction is fully counteracted by the solvation effect. It is noteworthy that the differences in entropy contribution (T∆S) of the six systems were not obvious, which implied the normal-mode analysis might be not accurate enough for Pim kinases. The molecular mechanisms of the Pim kinases lower inhibitory activity of AZD1208 compared to PIM447, and the better binding affinities of AZD1208 to Pim1 than to the other two kinases (Pim2 and 3), were illuminated in following the free energy decomposition and binding mode analysis.

Decomposition Analysis of the Binding Free Energies
MM/GBSA free energy decomposition analysis was employed to decompose the total binding free energies into ligand-residue pairs, which would provide more detailed information about the contribution of each residue of Pim kinases on the interacting with inhibitors. To investigate the direct interactions in the Pim/ligand complexes, 14 amino acid residues which contributed significantly (<−1 kcal/mol) to the binding of inhibitors to Pim kinases were selected ( Figure 3). Considering that the 14 amino acid residues could form the binding pocket of the inhibitor, these residues were naturally divided into five groups (    Alignment of amino acid sequences of the Pim kinases (Pim 1, 2, and 3) showed that among the 14 residues only Val126 (numbered in Pim1) in the left group was different ( Figure 5). The residue valine of Pim1 was replaced by alanine in Pim2 and 3, which resulted in the decreased contribution of the residue (Figure 3). However, the total contributions of the left group on six Pim/ligand systems were not significantly different ( Figure 6, Table 3). In other words, the weaker contribution of the residue alanine in Pim2 and 3 was nearly complemented by the other two residues (Arg and Pro) in Alignment of amino acid sequences of the Pim kinases (Pim 1, 2, and 3) showed that among the 14 residues only Val126 (numbered in Pim1) in the left group was different ( Figure 5). The residue valine of Pim1 was replaced by alanine in Pim2 and 3, which resulted in the decreased contribution of the residue (Figure 3). However, the total contributions of the left group on six Pim/ligand systems were not significantly different ( Figure 6, Table 3). In other words, the weaker contribution of the residue alanine in Pim2 and 3 was nearly complemented by the other two residues (Arg and Pro) in the left group (Figures 3 and 7). The difference of residue at the position 126 (numbered in Pim1) had little effect on the binding interactions between Pim kinases and pan-Pim inhibitors. the left group (Figures 3 and 7). The difference of residue at the position 126 (numbered in Pim1) had little effect on the binding interactions between Pim kinases and pan-Pim inhibitors.   First-order sequence alignment of three Pim kinases. Fourteen key residues calculated by MM/GBSA free energy decomposition were divided into five groups: top (green), left (blue), bottom (red), right (yellow), and loop group (pink). The blue arrows indicate beta strands and orange spirals indicate alpha helices. The black arrow indicate that the residues at this site differ among the three Pim kinases.     In order to further analyze the reasons for the different activity of Pim inhibitors, the differences in contribution values of the five groups were studied ( Figure 6, Table 3). As shown in Figures 6 and7, the top group had the largest contribution in all of the six systems, and the top, right, and left groups had a greater contribution to the sum binding free energies both of Pims/PIM447 and Pims/AZD1208. In order to further analyze the reasons for the different activity of Pim inhibitors, the differences in contribution values of the five groups were studied ( Figure 6, Table 3). As shown in Figures 6 and 7, the top group had the largest contribution in all of the six systems, and the top, right, and left groups had a greater contribution to the sum binding free energies both of Pims/PIM447 and Pims/AZD1208. It was obviously that the contribution of these three groups in the six Pims/ligands systems were similar, while the bottom and loop groups contributed differently to the sum binding free energies of Pims/PIM447 and Pims/AZD1208. The contribution value of each group in the three Pims/PIM447 was in the same order. The order from high to low was top, bottom, right, left, and loop group ( Figure 6, Table 3). However, the contribution of AZD1208 in Pim1 and Pim3 were ranked as top > right > Left > loop > bottom group. The order of contribution value of five groups in Pim2/AZD1208 system was top > right > left > bottom > loop group. It was notable to conclude that the difference between the contribution values of the bottom and loop groups was the main factor affecting the activity ( Table 3). The total contribution of the bottom group in each Pim/PIM447 system was about three times larger than that of the corresponding Pim/AZD1208 system (Table 3), which indicated that the bottom group was crucial for improving the activity of pan inhibitors. Although the loop group had only one amino acid residue in each system, the contribution of it varied greatly in the six systems. The maximal contribution value of the loop group was −2.47 kcal/mol in the Pim3/AZD1208 system and the minimal contribution value was −0.05 kcal/mol in the Pim2/AZD1208 system, which meant that the residue in the loop group had a significant influence on the activity of inhibitors.

Homology Model of Pim3
Considering that the crystal structure of Pim3 was not available, a homology model of Pim3 was carried out firstly. The amino acid sequence of th ePim3 protein was retrieved from UniProt (http://www.uniprot.org) with the UniProt entry of Q86V86. For template selection, the build_profile script was performed by searching against the pdball database using Modeller 9.20 [14]. Results from the template search showed Pim1 was the best template for construction of the 3D structure of Pim3 due to its high amino acid sequence identity (78%) to Pim3. The 3D structure of the template protein Pim1 can be obtained from the Protein Data Bank (PDB) (www.rcsb.org/pdb) with the PDB ID of 3UIX.
Aligning of amino acid sequences and construction of homology models were conducted by align2d script [15] and model-single script with default values in Modeller 9.20, respectively. The best model can be selected by picking the model with the highest value of the GA341 score (the highest score was 1.00000) and the lowest values of the molpdf and DOPE scores [16][17][18].

Optimization of Homology Model of Pim3
The newly built homology models often produce unfavorable atomic distances, bond angles, van der Waals radius overlap, and undesirable torsion angles. Therefore, it was essential to minimize the energy to regularize local bond and angle geometry as well as to relax close contacts in the geometric chain. The FG-MD (https://zhanglab.ccmb.med.umich.edu/FG-MD/) method [13] was employed to optimize the generated structure. Structural quality of the optimized protein was analyzed by ProSA web. Various other parameters including buried protein, quality factor, three-dimensional score, and non-bonded interaction were analyzed by SAVES (https://servicesn.mbi.ucla.edu/SAVES/).

Construction of Six Pim/Ligand Systems
The Surflex-Dock module of Sybyl X2.1 software (SYBYL_X2.1 is available from Tripos Associates Inc., S. H. R.; St. Louis, MO 631444, USA.) was used for molecular docking studies to predict the binding modes of AZD1208 and PIM447 to Pim kinases. The parameters of docking were set as default. Herein, six systems, i.e., Pim1/AZD1208, Pim1/PIM447, Pim2/AZD1208, Pim2/PIM447, Pim3/AZD1208, and Pim3/PIM447 were constructed. The crystal structure of Pim1/AZD1208 could not be obtained directly, while the ligand (7li) in the crystal structure of 4DTK was similar to AZD1208. Based on the binding pose of ligand 7li, the Pim1/AZD1208 system was constructed via the sketch module in Sybyl X2.1. The crystal structure of Pim1/PIM447 could be obtained directly (PDB ID: 5DWR). The crystal structure of Pim2 which was used for molecular docking was 4X7Q. The 3D structure of Pim3 was generated from the above homology modeling study. Then, AZD1208 and PIM447 were docked to Pim2 and 3, respectively.

Molecular Dynamics Simulations
MD simulations were performed on the above six systems using the Amber12 software package [19]. The tleap module was used to added missing hydrogen atoms of the six complexes. The proteins and the ligands were respectively parameterized using the standard amber force field (ff03) and general amber force field (gaff). Each ligand was minimized using the HF/6-31 * optimization in the Gaussian 09 program [20] and the atomic charges were subsequently fitted by using the restrained electrostatic potential (RESP). Then, each system was immersed into a cubic TIP3P water box extended 12 Å from any solute atom, and appropriate amounts of Na + were added to neutralize the system. Prior to MD simulations, the systems were optimized by the conjugate gradient and steepest descent methods as described in previous papers [21,22]. With the integration step time was given as 2 fs and the constant temperature was run at 310 K, 20 ns MD simulations were performed. The entire coordinate file was saved at each 1 ps for the following binding free energy calculation and free energy decomposition analysis.

MM/GBSA Calculation and MM/GBSA Free Energy Decomposition Analysis
For the purpose of elucidating the binding affinities of the six Pim/ligand systems, the MM/GBSA binding free energy calculation was conducted, and the method was calculated according to the following equation [23]: where ∆E MM is the gas-phase interaction energy that is composed of van der Waals (∆E vdw ) and electrostatic (∆E ele ) energies; ∆G GB is polar desolvation free energy that is calculated by the generalized Born (GB) model [24], and ∆G SA is the non-polar desolvation free energy which is estimated from the accessible surface area (SASA) model by the LCPO method: ∆G SA = 0.0072 × ∆SASA [25]. −T∆S is the conformational entropy contribution at temperature T. In the GB model, the solvent and the solute dielectric constants were set to 80 and 4, respectively. The binding free energy of each system was calculated based on 1000 snapshots evenly extracted from 10 to 20 ns MD trajectories. In view of the high computational demand of normal-mode analysis, only 100 snapshots extracted from the last 10 ns were used to estimate the entropy contribution. To validate the feasibility of the MM/GBSA methods for Pim/ligand systems, the MD simulation and MM/GBSA calculation studies were also applied to the other two Pim1 inhibitors GDC-0039 [8] and hispidulin [26], the potencies of which are lower than the compounds AZD1208 and PIM447. The detailed results about MD simulations and free energy calculations are described in the Supplementary Data. To uncover the key residues in the binding interactions between Pim kinases and inhibitors, MM/GBSA free energy decomposition analysis was carried out using the mm_pbsa module in Amber12. After the decomposition process, the total binding free energy contribution could be allocated to each residue from the association of the receptor and the ligand. All 400 snapshots generated for the binding free energy calculation were also used for the free energy decomposition analysis.

Conclusions
We studied the molecular recognitions of two pan-Pim inhibitors (PIM447 and AZD1208) with three Pim kinases via homology modeling, molecular docking, MD simulation and MM/GBSA binding free energy calculation. Fourteen amino acids which bound strongly to the inhibitor binding pocket were divided into five groups (top, left, bottom, right, and loop). The top group had the greatest contribution on the binding of pan-Pim inhibitors to Pim kinases. The bottom group contributed greater to Pims/PIM447 than to Pims/AZD1208, which gave rise to the higher activity of PIM447. The contribution of the loop group varied greatly in the six systems, which suggested it had an important impact on the activity of different inhibitors. The amino acids of the left group and the right group were located in the loop regions of the proteins, however, the flexibility of these regions gave the protein an advantage to accommodate diverse pan-Pim inhibitors. Not surprisingly, the contributions of the left and right group were similar both in Pims/PIM447 and Pims/AZD1208. It is likely that the more potent pan-Pim inhibitors should be bound strongly to the top and bottom groups. Our work could be helpful in clarifying the distinguishing roles of the 14 amino acids in improving the activity of inhibitors and designing novel inhibitors, which will be valuable for Pim-driven human cancers.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
Pim proviral integration site of the Moloney murine leukemia virus MD molecular dynamics MM/GBSA molecular mechanics generalized Born/surface area