J Mol Model (2014) 20:2075 DOI 10.1007/s00894-014-2075-1

ORIGINAL PAPER

Molecular dynamic simulations reveal the mechanism of binding between xanthine inhibitors and DPP-4 Yongliang Gu & Wei Wang & Xiaolei Zhu & Keke Dong

Received: 12 August 2013 / Accepted: 12 November 2013 # Springer-Verlag Berlin Heidelberg 2014

Abstract We apply molecular docking, molecular dynamics (MD) simulation, and binding free energy calculation to investigate and reveal the binding mechanism between five xanthine inhibitors and DPP-4. The electrostatic and van der Waals interactions of the five inhibitors with DPP-4 are analyzed and discussed. The computed binding free energies using MM-PBSA method are in qualitatively agreement with experimental inhibitory potency of five inhibitors. The hydrogen bonds of inhibitors with Ser630 and Asp663 can stabilize the inhibitors in binding sites. The van der Waals interactions, especially the key contacts with His740, Asn710, Trp629, and Tyr666 have larger contributions to the binding free energy and play important roles in distinguishing the variant bioactivity of five inhibitors. Keywords DPP-4 . MM-PBSA . Molecular dynamics simulation . Xanthine inhibitor

Introduction Type 2 diabetes is a chronic disease characterized by high levels of glucose, which results from insulin resistance and impairment of insulin secretion [1]. This metabolic disorder may lead to glucose toxicity caused by their hyperglycemia, which results in eye (or nerve, kidney, cardiovascular) damage [2]. According to the current level of biology and iatrology, people just control glucose to delay or suppress the disease Electronic supplementary material The online version of this article (doi:10.1007/s00894-014-2075-1) contains supplementary material, which is available to authorized users. Y. Gu : W. Wang : X. Zhu (*) : K. Dong State Key Laboratory of Materials-Oriented Chemical Engineering, College of Chemistry and Chemical Engineering, Nanjing University of Technology, Nanjing 210009, China e-mail: [email protected]

emergence, in which there are some treatment methods as follows: increasing insulin supplies with exogenous insulin [3, 4] and insulin secretion with sulfonylureas [5–7], limiting glucose absorption with glucosidase inhibitors [8–10], reducing insulin resistance using glitazones [11–16] and hepatic glucose output with biguanides [17, 18]. The above therapies can control the glucose homeostasis, however, they will bring about undesirable side effects such as hypoglycemia and weight gain [19]. Recently, with the development of genomics and deep understanding of type 2 diabetes, new targets and drugs of the disease have been found. Since 2005, novel treatment strategies of type 2 diabetes focused on hormone glucagonlike peptide 1(GLP-1(7–36) amide) have been proposed [20]. The GLP-1 is a 30-amino acid peptide hormone produced by L-cells of the intestinal mucosa in response to ingestion of food [21, 22]. When foods enter into intestinal canal, GLP-1 is released to loop and glucose gets into blood, in which both of them stimulate insulin [23]. Specifically, GLP-1 stimulates endogenous insulin release [24–26], inhibits glucagon secretion [27], slows gastric emptying [28, 29], and reduces appetite [30], which is favorable to controlling blood glucose. Unfortunately, active GLP-1(7–36) amide will be rapidly truncated to inactive GLP-1(9–36) amide by the cleavage enzyme of dipeptidyl peptidase 4 (DPP-4, E(3.4.14.5) [31]. DPP-4 is 110-kDa a serine protease, and belongs to the proyl oligopeptidase family [32–34]. DPP-4 is a membrane-bound exoprotease with a 22-amino acid hydrophobic transmembrane region and a cytoplasmic tail of only five amino acids at the N-terminus [35]. DPP-4 catalyzes the cleavage of GLP1(7–36) amides from their N-terminus with the sequence HX-Pro-Y or H-XAla-Y (where X, Y = any amino acid; Y ≠ Pro) [36, 37]. Therefore, DPP-4 is rapidly emerging as a new target for the treatment of type 2 diabetes [38]. In recent years, the design and synthesis of potent small molecule inhibitors to target DPP-4 have been the important

2075, Page 2 of 12

topics of many experimental studies [39–51]. The small molecule inhibitors of DPP-4 can be divided into two major classes, that is, non-peptide inhibitors and peptidomimetics inhibitors. The xanthine/phenylethylamine and cyanopyrrolidine are important non-peptide and peptidomimetics inhibitors, respectively. So far, only a few theoretical studies about the detailed binding mechanism of inhibitors with DPP-4 have been reported [50, 51]. Based on the above situation, our research group is going to systematically study the DPP-4 inhibitors mentioned above and compare their bioactivity in terms of molecular dynamics simulations [50–54]. In this work, we want to focus on the binding mechanism of xanthine inhibitors with DPP-4 and reveal the relationship between the structures of xanthine inhibitors and bioactivity, which will be helpful for the design of new potential DPP-4 inhibitors. Later, we will extend this work to the binding mechanism of peptidomimetics inhibitors with DPP-4. Herein, in order to obtain the detailed information about the binding of inhibitors to DPP-4, molecular docking, molecular dynamics (MD) simulations, and binding free energy calculations have been performed to reveal the binding mechanism of the five xanthine inhibitors [43, 49] and DPP-4. Results demonstrate that the van der Waals and net electrostatic energies dominate the binding of inhibitors to DPP-4. It is important to find that variant bioactivity of five xanthine inhibitors can be responsible for the key contacts of His740, Asn710, Trp629, and Tyr666. In current work, we reveal the reason why these five xanthine inhibitors have variant inhibitory potency, which will be helpful for understanding the binding modes of them with DPP-4 and provide some clues to the design of more promising inhibitors.

Computational details Initial structure preparation The initial structure of DPP-4 for docking and MD simulation calculations is obtained from the X-ray crystal structure (PDB ID: 10.2210) of the 1RWQ complex in the RCSB Protein Data Bank (PDB). In addition, five xanthine inhibitors are selected as DPP-4 inhibitors. All the five DPP-4 inhibitors are constructed by 3D graphical software. The inhibitors are optimized at the B3LYP/6-31G(d) level using the Gaussian 09 program [55]. The chemical structures of five inhibitors are shown in Table 1. Molecular docking The molecular dockings are performed by Autodock4.0 program applying the Lamarckian genetic algorithm [56]. The polar hydrogen atoms are added to DPP-4 and Kollman united

J Mol Model (2014) 20:2075

atom charges are assigned before docking [57]. During the docking process, DPP-4 is taken as rigid while positions and torsional bonds of each inhibitor are kept free for flexible docking. Before grid maps are calculated by AutoGrid, the three dimensions of the grid are set to 66 Å×65 Å×52 Å with 0.375 Å spacing value, therefore, there is sufficient space to include most of the protein and the active sites. During the docking process of inhibitors to DPP-4, the empirical free energy function and the Lamarckian genetic algorithm are employed. 200 independent runs are carried out with maximum of energy evaluations to 25,000,000 performed for each ligand and population size to 300. After clustering analysis, the structures with lower mean binding energy and the larger number of conformations are chosen as the preferred docking conformations as shown Fig. SI-1. Molecular dynamics simulation The missing hydrogen atoms of DPP-4 protein are added using the leap module in AMBER 10.0 software package [58]. The restrained electrostatic potential (RESP) method [59] is used to determine its partial atomic charges at the B3LYP/6-31G(d) level. The standard AMBER force field (ff03) [60] and the general AMBER force field (gaff) [61] are used to describe the interactions of protein and ligand, respectively. The counterions (Na+ ions) are added to make the system neutral. Each system is immersed in an octahedron periodic box with TIP3P [62] water model, and the minimum distance between each protein and box walls is set as 10 Å. Before MD simulations, in order to satisfy the convergence criterion, the energy minimization of each system is carried out with steepest descent method for 2000 steps using the sander program in AMBER10.0. After that, a position restrained dynamics simulation is performed on each system. Proteins, ligands, and water molecules are coupled separately to a temperature bath of 300 K with a coupling time of 0.1 ps. Finally, a 20 ns MD simulation is performed on each system. In MD simulations, the SHAKE procedure [63] restrains all covalent bonds containing hydrogen atoms, and particle mesh Ewald (PME) is applied to deal with the long-range electrostatic interaction [64]. Binding free energy calculation The binding free energies (ΔGbind) of the inhibitors with DPP4 are computed using molecular mechanics PoissonBoltzmann/generalized Born solvent accessible surface area MM-PBSA/GBSA [65] procedure in AMBER10. The binding free energy of a protein-ligand complex (ΔGbind) is calculated in terms of the following relations [66, 67]: ΔGbind ¼ ΔGðcomplexÞ−½ΔGðproteinÞ þ ΔGðligandފ ð1Þ

J Mol Model (2014) 20:2075

Page 3 of 12, 2075

Table 1 Structure and in vitro activity of xanthine inhibitors of DPP-4

ΔGbind ¼ ΔE gas þ ΔGsol −TΔS

ð2Þ

ΔE gas ¼ ΔEele þ ΔEvdW

ð3Þ

ΔGsol ¼ ΔGpolar þ ΔGnonpolar

ð4Þ

ΔGnonpolar ¼ γ ðSASAÞ þ β

ð5Þ

where ΔGcomplex, ΔGreceptor, and ΔGligand are the free energies of the complex, protein, and ligand, respectively. ΔEgas is the gas phase interaction energy, ΔGsol is the solvation free energy, and TΔS represents the entropy term. The gas phase interaction energy (ΔEgas) is composed of the electrostatic (ΔEele) and van der Waals (ΔEvdw). The solvation free energy (ΔGsol) includes the polar contribution (ΔGpolar) and the nonpolar contribution (ΔGnonpolar). The polar part (ΔGpolar)

can be calculated by solving the Poisson Boltzmann (PB) equations for MM-PBSA method or generalized Born (GB) model for MM-GBSA. The nonpolar solvation contribution (ΔGnonpolar) is determined using the solvent accessible surface area (SASA). The probe radius of the solvent is set to 1.4 Ǻ. The corresponding solvating parameters γ and β are 0.00542 kcal (mol Å)−1 and 0.92 kcal mol−1, respectively. If the ligands have similar structure, the entropy contributions are not considerably different for ligands binding to a same protein site [68]. Therefore, in current work, the conformational entropy (ΔSbind) is ignored [69].

Results and discussion As mentioned above, five DPP-4 inhibitors are docked to DPP-4 and the most favorable docking conformations are shown in Fig. 1. Clearly, five inhibitors locate in S1 binding pocket and N-terminal lobes of the protease. The detailed binding modes for five DPP-4/inhibitor complexes are

2075, Page 4 of 12

J Mol Model (2014) 20:2075

Fig. 1 The most favorable docking structures. M1∼M5 are represented by yellow, blue, magenta, red, and wheat, respectively

illustrated in Fig. SI-2. The docking details of Ki values, docking energies, binding energies are represented in Table SI-1. As displayed in Fig. 1 and Fig. SI-2, five inhibitors position themselves inside the hydrophobic pocket made of Arg125, Glu205, Glu206, Tyr547, Ser630, Tyr662, Tyr666, and His740, in which there are hydrogen bond interactions between inhibitors and residues Ser630, Glu205 and Glu206.

System stability and residue flexibility during MD simulations In order to examine the dynamics stability of five complexes after 20 ns MD simulations, the root-mean-squared deviations (RMSDs) relative to their starting structures of the backbone Cα atoms for the bound and free DPP-4 are analyzed as shown in Fig. 2. For each complex, the RMSD values rapidly increase during the first 5 ns, and then they keep stable with the RMSD fluctuation values of ∼1.4 Å, implying that the five complexes reach stability. According to the above situation, it is reasonable to do the subsequent analyses. On the other hand, the gyration radii (Rg) of five complexes are computed as shown in Fig. SI-3, where Rg is a parameter of molecular

Fig. 2 The time dependence of RMSDs for the five DPP-4/ inhibitor complexes relative to the initial minimized structures during the MD simulations

size [70]. The average values of five complexes are approximately the same as that (26.7–27.3 Å) of free DPP-4, which reveals that the structures of all complexes are stable [71, 72]. To examine the residue flexibility of the bound DPP-4, the B-factors [73] for each system are calculated and compared with experimental values. Figure SI-4 shows the calculated Bfactors of inhibitor-bound DPP-4 and experimental values obtained from X-ray crystal structure of DPP-4. The residue flexibility of the bound DPP-4 is analyzed based on the root mean square fluctuation (RMSF) of the Cα atoms of each residue as displayed in Fig. SI-5. It is noted that these complexes with the different inhibitors have similar RMSF profiles. However, for DPP-4/M1, the residues near Ile518Pro531, Tyr540-Ala548, Gly587-Val603 and Gln612Gly633 exhibit higher flexibilities.

Binding mechanism of DPP-4/inhibitor complexes and bioactivity of inhibitors In order to understand the different biological activities of five inhibitors and examine the interactions between inhibitors and DPP-4, the binding free energies are calculated for five DPP/ inhibitor complexes using the MM-PBSA method. The results

J Mol Model (2014) 20:2075

Page 5 of 12, 2075

Table 2 Binding free energy computed by MM-PBSA for five complexes (kcal mol−1) Component

DPP-4/M1

DPP-4/M2

DPP-4/M3

DPP-4/M4

DPP-4/M5

ΔEele ΔEvdW ΔEgas ΔGpol ΔGnonpolar ΔGsol ΔGbind IC50 [45, 46]

−21.82 (3.09) −51.55 (2.34) −73.37 (3.63) 56.95 (4.40) −5.71 (0.08) 51.24 (4.37) −22.13 (4.31) 5

−38.13 (4.23) −35.02 (2.50) −73.16 (4.68) 56.35 (4.09) −4.58 (0.08) 51.77 (4.08) −21.39 (3.44) 35

−11.71 (2.53) −44.50 (2.27) −56.21 (3.78) 42.00 (3.37) −5.40 (0.18) 36.60 (3.31) −19.61 (3.59) 284

−34.06 (3.54) −37.03 (2.58) −71.09 (3.51) 56.09 (4.87) −4.76 (0.14) 51.33 (4.85) −19.76 (4.64) 57

−21.01 (3.94) −24.46 (2.13) −45.48 (4.55) 32.12 (4.20) −3.78 (0.18) 28.35 (4.16) −17.13(2.17) >500

Fig. 3 Decomposition of ΔG on a per-residue basis for the protein-inhibitor complexes

2075, Page 6 of 12

J Mol Model (2014) 20:2075

Fig. 4 2D representation of hydrogen bond and hydrophobic interactions by LIGPLOT program analyses. Dashed lines represent hydrogen bonds, and spiked residues account for the hydrophobic interactions of DPP-4 with inhibitors. a DPP-4/M1; b DPP-4/M2; c DPP-4/M3; d DPP-4/M4,e DPP-4/M5

J Mol Model (2014) 20:2075

of computed free energies and energy components of each complex are listed in Table 2. It is not difficult to note from the structures of inhibitors in Table 1 and binding free energies in Table 2 that the bioactivities of five inhibitors are very sensitive to the substituted groups (R1 and R2) or “X” atom. As shown in Table 2, the calculated binding free energies of DPP4/M1 ∼ DPP-4/M5 are −22.13, −21.39, −19.61, −19.76, −17.13 kcal mol−1, respectively, suggesting the biological activity of M1>M2>M3, and M4>M5, which are consistent with the experimental bioactivity data (IC50 values) [43, 49]. On the other hand, the binding free energy decomposition is performed to analyze the contributions of each residue to the binding. The total energies per inhibitor–residue pair of five complexes (DPP-4/M1∼DPP-4/M5 complexes) are shown in

Page 7 of 12, 2075

Fig. 3, which demonstrates that the key residues for binding interactions between DPP-4 and inhibitors are as follows, DPP-4/M1(G205 (Gly205), Y547 (Tyr547), W629 (Trp629), Y666 (Tyr666), N710 (Asn710), V711 (Val711), and H740 (His740)), DPP-4/M2(G205 (Gly205), Y547 (Tyr547), Y662 (Tyr662), Y666 (Tyr666), N710 (Asn710), V711 (Val711) and H740 (His740)), DPP-4/M3(F357 (Phe357), Y547 (Tyr547), S552 (Ser552), Y631 (Tyr631), Y662 (Tyr662) and Y666 (Tyr666)), DPP-4/M4(R125 (Arg125), G205 (Gly205), Y547 (Tyr547), W629 (Trp629) and Y666 (Tyr666)), DPP-4/M5 (Y547 (Tyr547) and W629 (Trp629)). It is worth noting that the residues Ser630, Glu205 and Glu206 seem to have low magnitudes of residue interactions in the plot of the decomposition of ΔG on a per-residue basis

Fig. 5 Van der Waals interaction energy spectra of inhibitor-residue pair in DPP-4/inhibitor complexes

2075, Page 8 of 12

J Mol Model (2014) 20:2075

(Fig. 3) after MD simulations. It is due to the fact that during MD simulations, the initial conformations (the preferred docking conformations) of M1/DPP-4∼M5/DPP-4 will be changed with the time evolution, leading to the changes of binding sites and interactions. To obtain the details about the binding mechanism of DPP/inhibitor complexes, the 2D schemes of protein-inhibitor complexes are generated by LIGPLOT program [74] from MD simulation results, and the hydrogen bonding and hydrophobic interactions are examined as shown in Fig. 4, in which the key residues of each DPP-4/inhibitor complex are similar to those obtained from the binding free energy decomposition as mentioned above (Fig. 3). In order to deeply understand which energy terms are more favorable to the binding of each DPP-4/inhibitor complex, it is important to compare the individual energy components (ΔEele, ΔEvdW, ΔGpol, ΔGnonpolar). Obviously, Table 2 reveals that the van der Waals energy term (ΔEvdW) has larger contribution to the total binding free energy (ΔGbind) and is the main reason for different bioactivity of five inhibitors. However, the polar solvation energy term (ΔGpol) is considerably unfavorable for the binding in five complexes. The gas phase electrostatic energy term (ΔEele) has favorable contribution to binding, but it still does not completely cancel the negative effect caused by the polar solvation free energy (ΔGpol) [75, 76]. The nonpolar solvation free energies (ΔGnonpolar) slightly drive the binding. Interestingly, the correlation coefficient (r2) between the ΔG(ΔEvdW +ΔEele+ ΔGpol) and total predicted free energy (ΔGbind) is up to 0.81, which implies that the van der Waals energy term and electrostatic energy term are dominant for differentiating the binding affinity of five inhibitors. In order to reveal how the van der Waals energy terms (ΔEvdW) impact on the binding of five complexes, the

van der Waals interaction energy spectra of inhibitorresidue pair in DPP-4/inhibitor complexes are analyzed as shown in Fig. 5, suggesting that DPP-4/M1 and DPP4/M2 have the almost same key residues, but the van der Waals interactions (ΔEvdW, −2∼−3 kcal mol−1) between key residues and M1 are larger than those (−1∼−2 kcal mol−1) between key residues and M2, leading to larger bioactivity of M1 than M2. As shown in Fig. 5, the number of key residues for DPP-4/M1 is significantly larger than that of DPP-4/M3, which results in the stronger van der Waals interactions (ΔEvdW) of M1 than those of M3. All DPP-4/M1∼DPP-4/M3 have major residues of Y547, W629, Y662, Y666, V711, and H740 as shown in Fig. 5. The above analyses illustrate that van der Waals energy terms make the biological activity ordering of M1>M2 (or M3). As mentioned above, the electrostatic interactions also affect the binding affinity of five inhibitors. As shown in Fig. 4, there are two stable hydrogen bonds between M1 and its adjacent residues in the DPP/M1 complex. One hydrogen bond is established between the oxygen atom (O1) in the pyrimidine ring of M1 and the oxygen atom of Ser630, and another hydrogen bond forms between the nitrogen atom (N6) of M1 and oxygen atom of Asp663. Hydrogen bonds formed and their occupancies are displayed in Table 3 during the 20 ns simulations. As shown in Table 3, the two hydrogen bonds are stable during the whole MD simulation with average distances of 2.64 and 3.10 Å for O1 (M1)-OG (Ser630) and N6(M4)OD(Asp663), respectively, as represented in Fig. 4. The results of the hydrogen bonds occupancies during the MD simulations (Table 3) indicate that the M1 can form stable hydrogen bonds with Ser630 (66.89 %) and Asp663 (62.41 %), respectively. These stable hydrogen bonds result in stronger electrostatic interactions between M1 and DPP-4.

Table 3 Hydrogen bonds formed between inhibitors and DPP-4 during MD simulations a System

Donor

Acceptor

Occupancy (%)

Distance (Å)

Angle( )

DPP-4/M1

:M1@O1 :663@OD1 :M2@O1 :M2@O1 :663@OD2 :663@OD1 :M2@O1 :710@OD1 :M3@N4 :M4@OD1 :M5@O2 :M5@O2 :M5@O2

:630@OG-:630@HG :M1@N6-:M1@H30 :630@OG-:630@HG :740@NE2-:740@HE2 :M2@N6-:M2@H25 :M2@N6-:M2@H25 :547@OH-:547@HH :M3@N6-:M3@H30 :666@OH-:666@HH :628@H-:628@N :554@NZ-:554@HZ3 :554@NZ-:554@HZ1 :554@NZ-:554@HZ2

66.89 62.41 80.27 45.59 38.73 37.74 18.06 18.80 8.87 16.79 31.71 22.57 21.54

2.784(0.17) 3.085(0.19) 2.717(0.14) 3.048(0.19) 3.012(0.18) 3.021(0.18) 2.996(0.22) 3.089(0.20) 3.248(0.16) 3.090 (0.21) 2.853(0.13) 2.857(0.14) 2.860(0.14)

20.72(12.41) 18.13(10.03) 18.65(9.71) 31.03(14.01) 17.76(9.07) 18.79(9.61) 35.25(13.26) 28.60(13.95) 46.15(10.94) 25.78(12.30) 24.20(12.44) 25.50(12.64) 24.65(12.80)

DPP-4/M2

DPP-4/M3 DPP-4/M4 DPP-4/M5

a

The hydrogen bonds are determined by acceptor/donor atom distance of less than 3.5 Å and acceptor/H–donor angle of larger than 120º

J Mol Model (2014) 20:2075

Page 9 of 12, 2075

M2 also has hydrogen bonding interactions with Ser630 and Asp663 with occupancies of 80.27 % and 38.73 %, respectively. In addition, M2 forms weak hydrogen bonds with residues (H740 and Y547) as shown in Table 3. It leads to ΔEele of M2 is larger than that of M1 as shown in Table 2. As shown in Fig. 4, for DPP/M3 complex, there is one hydrogen bond between the nitrogen atom (N6) of the thiazole ring of M3 and the oxygen atom of Tyr662(M2 (or M3) and M4>M5. Although ΔEgas (ΔEele+ΔEvdW) of M1 is almost same as that of M2 as shown in Table 2, in fact, the larger biological activity M1 than M2 is ascribed to the non-polar interactions between DPP-4/M1 complex and solvent. As shown in Table 2, the ΔGnonpolar of M1 is −5.71 kcal mol−1 while the ΔGnonpolar of M2 is −4.58 kcal mol−1. It is due to the fact that M1 has an extra benzene ring in R1 group, resulting in that M1 has a lower dielectric constant that M2. As shown in Fig. 7, the ΔGnonpolar of the key residues of M1 are significantly favorable for total binding free energy while the ΔGnonpolar of the key residues of M2 are almost unfavorable, which is main reason for larger bioactivity of M1 than M2. To gain some insights into the difference in binding modes of DPP-4/M1∼DPP-4/M3 complexes, the binding modes of inhibitors M1, M2, and M3 in the active site of DPP-4 are displayed in Fig. SI-6. It is not difficult to note from Fig. SI-6(c) that the double bond of R2 in M3 can form van der Waals interaction with the side chain of Tyr666, and there is the π-π interaction of the benzene ring of R1 in M3 with the side chain of Tyr547. Moreover, Tyr631 and Ser552 interact with the

pyrimidine ring of M3 and the benzene ring at R1, respectively. As shown in Fig. SI-6(b), it is can be seen that the key residue His740 is around the double bond of R2 in M2. Since key residue His740 residue occupies corresponding site of Tyr666 of M3, the position change of Tyr666 leads to the orientation change of double bond R2 in M2, strengthening van der Waals interaction with key residue His740. As shown in Fig. 4 and Fig. SI-6(a), there is strong π-π interaction of key residue Trp629 with benzene ring of R1 in M1. Tyr666 has van der Waals interaction with the double bond of R2 in M1. In addition, there are strong van der Waals interactions of key residues Asn710 and His740 with imidazole ring and with pyrimidine ring of M1, respectively. The above analyses demonstrate that total interactions with His740, Asn710, Trp629, and Tyr666 could be the major factor for the variant bioactivity of M1 and M2 (or M3).

Fig. 6 Comparison of interaction energy based on the inhibitor-Ser630 (Tyr662 or Asp663) pair in the five DPP-4/inhibitor complexes

Fig. 8 Comparison of interaction energy based on the major inhibitorresidue pair in the three DPP-4/inhibitor complexes

Fig. 7 Comparison of nonpolar energies based on the major inhibitorresidue pair in the DPP-4/M1 and DPP-4/M2 complexes

2075, Page 10 of 12

To further confirm above inference about biological activity of M1∼M3, the contributions of major residues of three complexes (DPP-4/M1∼DPP-4/M3) to binding are compared in Fig. 8. Clearly, M1 has significantly stronger interactions with W629 and H740 while the interaction energy terms between other major residues (except for Y547) and M1∼ M3 are approximately similar, resulting in biological activity ordering of M1>M2 (or M3). As represented in Fig. 8, except for Y547, inhibitor-residue pair interactions between M2 and other major residues (W629, Y662, Y666, V711, and H740) are significantly stronger than those between M3 and the same residues, making biological activity ordering of M2>M3. As shown in Table 1, from M4 to M5, “X” atom is changed from nitrogen to carbon, and R1 and R2 are fixed. As shown in Fig. 5, the van der Waals energy terms of the residue Y547 are similar for DPP-4/M4 and DPP-4/M5. However, the major residue number of M4 is significantly larger than that of M5, leading to the stronger interactions between M4 and major residues and stronger biological activity. Figure SI-7 represents the decomposition of ΔGgas+sol for the key residues of the five complexes. Although the ΔGgas+sol values of residue Y547 are similar for DPP-4/M4 and DPP-4/M5, ΔGgas+sol values of other key residues of M4 are remarkably larger than those of M5, which reveals the reason for ten-fold biological activity difference between M4 and M5. On the other hand, as illustrated in Table 2, ΔG(ΔEvdW +ΔEele+ ΔGpol) values of DPP-4/M4 and DPP-4/M5 are −15.00 and −13.36 kcal mol−1, respectively, and ΔGnonpolar values are −4.76 and −3.78 kcal mol−1, resulting in larger bioactivity of M4 than M5. It should be noted that for M3 and M4, both R1 and R2 are different. Their IC50 values are drastically different while the binding free energies of M3 and M4 are very close, that is, the current MD simulations cannot distinguish their bioactivities, which may be ascribed to the immature potential functions used in this work.

Conclusions Summarily, molecular docking, molecular dynamics simulation, and binding free energy analysis are performed on five DPP-4/inhibitor complexes to reveal the binding mechanism of five xanthine inhibitors and DPP-4 and distinguish the bioactivity from high active inhibitors to low ones. The computed binding free energies of DPP-4/inhibitor complexes are qualitatively consistent with experimental biological data(IC50) of five inhibitors. Results display that the van der Waals and net electrostatic energies are dominant for the binding of inhibitors to DPP-4. The hydrogen bonds with Ser630 and Asp663 slightly derive binding of inhibitor/DPP-4. Interestingly and importantly, the bioactivity variation of five inhibitors can be ascribed to the key contacts of His740, Asn710, Trp629, and

J Mol Model (2014) 20:2075

Tyr666. The current work will be significant for the design of novel DPP-4 inhibitors in future. Acknowledgments This work is supported by grants from the National Science Foundation of China (Nos. 21276122, 21136001, and 20876073) and State Key Laboratory of Materials-Oriented Chemical Engineering, College of Chemistry and Chemical Engineering, Nanjing University of Technology of China (No. ZK201212).

References 1. Feng J, Zhang Z, Michael B (2007) Discovery of alogliptin: a potent, selective, bioavailable, and efficacious inhibitor of dipeptidyl peptidase IV. J Med Chem 50:2297–2300 2. Wright S, Mark J (2006) cis-2,5-Dicyanopyrrolidine inhibitors of dipeptidyl peptidase IV: synthesis and in vitro, in vivo, and x-ray crystallographic characterization. J Med Chem 49:3068–3076 3. Rosenstock J, Schwartz SL, Clark CM Jr, Park GD, Donley DW, Edwards MB (2001) Basal insulin therapy in type 2 diabetes: 28week comparison of insulin glargine (HOE 901) and NPH insulin. Diabetes Care 24:631–636 4. Barnett AH, Owens DR (1997) Insulin analogues. Lancet 349:47–51 5. King P, Peacock I, Donnelly R (1999) The UK Prospective Diabetes Study (UKPDS): clinical and therapeutic implications for type 2 diabetes. Br J Clin Pharmacol 48:643–648 6. UK Prospective Diabetes Study Group (1990) UK prospective diabetes study 7: response of fasting plasma glucose to diet therapy in newly presenting type II diabetic patients. Metabolism 39:905–912 7. Gerich JE (1989) Oral hypoglycemic agents. N Engl J Med 321: 1231–1245 8. Lillelund H, Jensen H, Liang X, Bols M (2002) Recent developments of transition-state analogue glycosidase inhibitors of non-natural product origin. Chem Rev 102:515–553 9. Berecibar A, Grandjean C, Siriwardena A (1999) Synthesis and biological activity of natural aminocyclopentitol glycosidase inhibitors: mannostatins, trehazolin, allosamidins, and their analogs. Chem Rev 99:779–844 10. Jacob S (1995) Glycosylation inhibitors in biology and medicine. Curr Opin Struct Biol 5:605–611 11. Saltiel R, Olefsky M (1996) Thiazolidinediones in the treatment of insulin resistance and type II diabetes. Diabetes 45:1661–1669 12. Schoonjans K, Martin G, Staels B, Auwerx J (1997) Peroxisome proliferator-activated receptors, or-phans with ligands and functions. Curr Opin Lipidol 8:159–166 13. Pittas G, Greenberg S (2002) Thiazolidinediones in the treatment of type 2 diabetes. Expert Opin Pharmacother 3:529–540 14. Rami K, Smith A (2000) Synthetic ligands for PPARç s review of patent literature 1994–1999. Expert Opin Ther Pat 10:623–634 15. Van Gaal L, Scheen J (2002) Are all glitazones the same diabetes. Metab Res Rev 18:S1–S4 16. Hauner H (2002) The mode of action of thiazolidinediones. Diabetes Metab Res Rev 18:S10–S15 17. Witters L (2001) The blooming of the French lilac. J Clin Invest 108: 1105–1107 18. Bailey J, Path C, Turner C (1996) Drug therapy-metformin N. Engl J Med 334:574–579 19. Wagman AS, Nuss JM (2001) Current therapies and emerging targets for the treatment of diabetes. Curr Pharm Des 7:417–450 20. Deacon CF (2011) Dipeptidyl peptidase-4 inhibitors in the treatment of type 2 diabetes:a comparative review. Diabetes Obes Metab 13:7–18

J Mol Model (2014) 20:2075 21. Drucker DJ (2001) Minireview: the glucagon-like peptides. Endocrinology 142:521–527 22. Ross SA, Gulve EA, Wang M (2004) Chemistry and biochemistry of type 2 diabetes. Chem Rev 104:1255–1282 23. Bo A (2009) Islet G protein-coupled receptors as potential targets for treatment of type 2 diabetes. Nat Rev Drug Discov 8:369–385 24. Mojsov S, Weir GC, Habener JF (1987) Insulinotropin: Glucagonlike peptide I (7–37) co-encoded in the glucagon gene is a potentstimulator of insulin release in the perfused rat pancreas. J Clin Invest 79:616–619 25. Kreymann B, Ghatei MA, Williams G, Bloom SR (1987) Glucagon-like peptide-1 7–36: a physiological incretin in man. Lancet 2:1300–1304 26. Rskov C, Holst JJ, Nielsen OV (1988) Effect of truncated glucagonlike peptide-1 [proglucagon-(78–107) amide] on endocrine secretionfrom pig pancreas, antrum, and nonantral stomach. Endocrinology 123:2009–2013 27. Nauck MA, Heimesaat MM, Behle K (2002) Effects of glucagonlike peptide 1 on counterregulatory hormone responses, cognitive functions, and insulin secretion during hyperinsulinemic, stepped hypoglycemic clamp experiments in healthy volunteers. J Clin Endocrinol Metab 87:1239–1246 28. Wettergren A, Schjoldager B, Mortensen PE (1993) Truncated GLP1 (proglucagon 78–107-amide) inhibits gastric and pancreatic functions in man. Dig Dis Sci 38:665–673 29. Nauck MA, Niedereichholz U, Ettler R, Holst JJ, Orskov C, Ritzel R, Schmiegel WH (1997) Glucagon-like peptide 1 inhibition of gastric emptying outweighs its insulino-tropic effects in healthy humans. Am J Physiol 273:E981–E988 30. Flint A, Raben A, Ersboll AK, Holst JJ, Astrup A (2001) The effect of physiological levels of glucagons-like peptide-1 on appetite, gastric emptying, energy and substrate metabolism in obesity. Int J Obes 25:781–792 31. Magnin DR, Robl JA (2004) Synthesis of novel potent dipeptidyl peptidase IV inhibitors with enhanced chemical stability: interplay between the N-terminal amino acid alkyl side chain and the cyclopropyl group of α-aminoacyl-L-cis-4,5-methanoprolinenitrilebased inhibitors. J Med Chem 47:2587–2598 32. Hiramatsu H, Yamamoto A, Kyono K (2004) The crystal structure of human dipeptidyl peptidase IV(DPPIV) complex with diprotin. A Biol Chem 385(6):561–564 33. Lambei AM, Seharpe, DeMeester (2008) DPP-4 inhibitors for diabetes- what next. Biochem Pharmacol 76(12):1637–1643 34. Mentlein R, Gallwitzand B, Schmidt WE (1993) Dipeptidylpeptidlase IV hydrolyses gastric inhibitory polypeptide, glueagonlike PePtide- l (7–36) amide peptide histidine methionine and 15 responsible for their degradation in human serum. EurJBioehem 214(3):829–835 35. Chien C-H, Tsai C-H, Lin C-H (2006) Identification of hydrophobic residues critical for DPP-IV dimerization. Biochemistry 45(23): 7006–7012 36. Hansen L, Deacon CF, Orskov C, Holst JJ (1999) Glucagon- like peptide-1 (7–36) amide is transformed to glucagon-like peptide-1 (9– 36) amide by dipeptidyl peptidase IV in the capillaries supplying the L-cells of the porcine intestine. Endocrinology 140:5356–5363 37. Kieffer TJ, McIntosh CHS, Pederson RA (1995) Degradation of glucose-dependent insulinotropic polypeptide and truncated glucagon-like Pep- tide-1 in vitro and in vivo by dipeptidyl peptidase IV. Endocrinology 136:3585–3597 38. Gwaltney SL, Stafford JA (2005) Inhibitors of dipeptidyl peptidase. Annu Rep Med Chem 40:149–165 39. Hiramatsu H, Yamamoto A, Kyono K (2004) The crystal structure of human dipeptidyl peptidase IV(DPPIV) complex with diprotin A. Biol Chem 385:561–564 40. Chien C-H, Tsai C-H, Lin C-H (2006) Identification of hydrophobic residues critical for DPP-IV dimerization. Biochemistry 45:7006–7012

Page 11 of 12, 2075 41. Engel M, Hoffmann T, Manhart S (2006) Rigidity and flexibility of dipeptidyl peptidase IV: crystal structures of and docking experiments with DPIV. J Mol Biol 355:768–783 42. Barnett A (2006) DPP-4 inhibitors and their potential role in the management of type 2 diabetes. Int J Clin Pract 60:1454– 1470 43. Eckhardt M, Langkopf E, Mark M (2007) 8-(3-(R)-aminopiperidin1-yl)-7-but-2-ynyl-3-methyl-1-(4-methyl-quinazolin-2-ylmethyl)-3, 7-dihydropurine-2,6-dione (BI 1356), a highly potent, selective, long-acting, and orally bioavailable DPP-4 inhibitor for the treatment of type 2 diabetes. J Med Chem 50:6450–6453 44. Li C, Weiqiang L, Chunhua L (2012) Identification of diverse dipeptidyl peptidase IV inhibitors via structure-based virtual screening. J Mol Model 18:4033–4042 45. Li C, Shen J, Li W (2011) Possible ligand release pathway of dipeptidyl peptidase IV investigated by molecular dynamics simulations. Wiley-Liss Inc 79:1800–1809 46. Pradip Jadav, Rajesh Bahekar, Shailesh R. Shah (2012) Long-acting peptidomimetics based DPP-IV inhibitors. Bioorg Med Chem Lett 47. Janey J, Hsiao Y, Joseph Armstrong III (2006) Proline-catalyzed, asymmetric mannich reactions in the synthesis of a DPP-IV inhibitor. J Org Chem 71:390–392 48. Green BD, Flatt PR, Bailey CJ (2006) Dipeptidyl peptidase IV (DPP IV) inhibitors: a newly emerging drug class for the treatment of type 2 diabetes. Vasc Dis Res 3:159–165 49. Gue Hao, Yanzheng Cai, Weicheng Zhou (2009) Progress in the research of dipeptidyl peptidase-4 inhibitors. World Clinical Drugs 30 NO.8 50. Liang D, Gao J, Cheng Y, Cui W, Zhang H, Ji M (2011) Molecular dynamics simulations and MM/GBSA methods to investigate binding mechanisms of amino methyl pyrimidine inhibitors with DPP-IV. Bioorg Med Chem Lett 21:6630–6635 51. Wei C, Desheng L, Jian G, Fang L, Lingling G, Mingjuan J (2013) Molecular dynamics and free energy studies of chirality specificity effects on aminobenzo[a]quinolizine inhibitors binding to DPP-IV. J Mol Model 19:1167–1177 52. Otyepka M, Skopalík J (2007) What common structural features and variations of mammalian P450s are known to date? Biochim Biophys Acta 1770:376–389 53. Vincent BC, Tan BZ, Lim KM, Tay TE (2010) Explaining the inhibition of cyclin-dependent kinase 5 by peptides derived from p25 with molecular dynamics simulations and MM-PBSA. J Mol Model 16:1–8 54. Chen Q, Cui W, Ji M (2009) Studies of chirality effect of 4(phenylamino)-pyrrolo[2,1-f] [1,2,4]triazine on p38a by molecular dynamics simulations and free energy calculations. J Comput Aided Mol Des 23:737–745 55. Frisch MJ, Trucks GW, Schlegel HB, Scuseria GE, Robb MA, Cheeseman JR, Scalmani G, Barone V, Mennucci B, Petersson GA, Nakatsuji H, Caricato M, Li X, Hratchian HP, Izmaylov AF, Bloino J, Zheng G, Sonnenberg JL, Hada M, Ehara M, Toyota K, Fukuda R, Hasegawa J, Ishida M, Nakajima T, Honda Y, Kitao O, Nakai H, Vreven T, Montgomery JA Jr, Peralta JE, Ogliaro F, Bearpark M, Heyd JJ, Brothers E, Kudin KN, Staroverov VN, Kobayashi R, Normand J, Raghavachari K, Rendell A, Burant JC, Iyengar SS, Tomasi J, Cossi M, Rega N, Millam JM, Klene M, Knox JE, Cross JB, Bakken V, Adamo C, Jaramillo J, Gomperts R, Stratmann RE, Yazyev O, Austin AJ, Cammi R, Pomelli C, Ochterski JW, Martin RL, Morokuma K, Zakrzewski VG, Voth GA, Salvador P, Dannenberg JJ, Dapprich S, Daniels AD, Farkas O, Foresman JB Ortiz JV, Cioslowski J, Fox DJ (2009) Gaussian 09 Revision.A.02. Gaussian Inc, Wallingford, CT 56. Morris GM, Goodsell DS, Halliday RS, Huey R, Hart WE, Belew RK, Olson AJ (1998) Automated docking using a Lamarckian genetic algorithnm and an empirical binding free energy function. J Comput Chem 19:1639–1662

2075, Page 12 of 12 57. Morris GM, Goodsell DS, Huey R, Olson AJ (1996) Distributed automated docking of flexible ligands to proteins parallel applications of AutoDock 2.4. J Comput Aided Mol Des 10:293–304 58. Case DA, Cheatham TA, Simmerling CL, Wang J, Duke RE, Luo R, Crowley M, Walker RC, Zhang W, Merz KM, Wang B, Hayik S, Roitberg A, Seabra G, Kolossvary I, Wong KF, Paesani F, Vanicek J, Wu X, Bronzell SR, Steinbrecher T, Gohlke H, Yang L, Tan C, Mongan J, Hornak V, Cui G, Mathews DH, Seetin MG, Sagui C, Babin V, Kollman PA (2008) AMBER 10. University of California, San Francisco, CA 59. Bayly CI, Cieplak P, Cornell WD, Kollman PA (1993) A wellbehaved electrostatic potential based method using charge restraints for determining atom-centered charges: the RESP model. J Phys Chem 97:10269–11028 60. Cornell WD, Cieplak P, Bayly CI, Gould IR, Merz KM Jr, Ferguson DM, Spellmeyer DC, Fox T, Caldwell JW, Kollman PA (1995) A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. J Am Chem Soc 117:5179–5197 61. Wang J, Wolf RM, Caldwell JW, Kollman PA, Case DA (2004) Development and testing of a general Amber force field. J Comput Chem 25:1157–1174 62. Jorgensen WL, Chandrasekhar J, Madura JD, Impey RW, Klein ML (1983) Comparison of simple potential functions for simulating liquid water. J Chem Phys 79:926–935 63. Ryckaert JP, Ciccotti G, Berendsen HJC (1977) Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J Comput Phys 23:327–341 64. Darden T, York D, Pedersen L (1998) Particle mesh Ewald: an Nlog(N) method for Ewald sums in large systems. J Chem Phys 98: 10089–10092 65. Srinivasan J, Cheatham TE, Cieplak P, Kollman PA, Case DA (1998) Continuum solvent studies of the stability of DNA, RNA, and Phosphoramidate-DNA Helices. J Am Chem Soc 120:9401–9409 66. Saíz-Urra L, Cabrera MA, Froeyen M (2011) Exploring the conformational changes of the ATP binding site of gyrase B from Escherichia colicomplexed with different established inhibitors by using molecular dynamics simulation: Protein – ligand interactions in

J Mol Model (2014) 20:2075

67.

68.

69.

70.

71.

72.

73.

74.

75.

76.

the light of the alanine scanning and free energy decomposition methods. J Mol Graph Model 29:726–739 El-Barghouthi MI, Jaime C, Al-Sakhen NA, Issa AA, Abdoh AA, Omari MM, Badwan AA, Zughul MB (2008) Molecular dynamics simulations and MM—PBSA calculations of the cyclodextrin inclusion complexes with 1-alkanols, para-substituted phenols and substituted imidazoles. J Mol Struct THEOCHEM 853:45–52 Cong XJ, Tan JJ, Liu M, Chen W, Wang CX (2010) Computational study of binding mode for N-substituted pyrrole derivatives to HIV-1 gp41. Prog Biochem Biophys 37(8):904–91566 Andricioaei I, Karplus M (2001) On the calculation of entropy from covariance matrices of the atomic fluctuations. J Chem Phys 115: 6289–6292 Lobanov MY, Bogatyreva NS, Galzitskaya OV (2008) Radius of gyration as an indicator of protein structure compactness. Mol Biol 42:623–628 Vivier I, Marguet D, Naquet P, Bonicel J, Black D, Li CXY, Bernard AM, Gorvel JP, Pierres M (1991) Evidence that thymocyte-activating molecule is mouse cd26 (dipeptidyl peptidase-iv). J Immunol 147: 447–454 Tanaka S, Murakami T, Horikawa H, Sugiura M, Kawashima K, Sugita T (1997) Suppression of arthritis by the inhibitors of dipeptidyl peptidase IV Int. J Immunopharmacol 19:15–24 Bristol LA, Sakaguchi K, Appella E, Doyle D, Taka’cs L (1992) Thymocyte costimulating antigen is CD26 (Dipeptidyl peptidaseIV)–costimulation of granulocyte, macrophage, and T-lineageceliproliferation via CD26. J Immunol 149:367–372 Wallace AC, Laskowski RA, Thornton JM (1995) LIGPOLT: a program to generate schematic diagrams of protein-ligand interactions. Protein Eng 8:127–134 Hu GD, Zhu T, Zhang SL, Wang D, Zhang QG (2010) Some insights into mechanism for binding and drug resistance of wild type and I50V V82A and I84V mutations in HIV-1 protease with GRL-98065 inhibitor from molecular dynamic simulations. Eur J Med Chem 45:227–235 Wu EL, Han K, Zhang JZ (2008) Selectivity of neutral/weakly basic P1 group inhibitors of thrombin and trypsin by a molecular dynamics study. Chem Eur J 14:8704–8714