The Molecular Modeling of Novel Inhibitors of Protein Tyrosine Phosphatase 1B Based on Catechol by MD and MM-GB (PB)/SA Calculations

Abstract

Binding modes of a series of catechol derivatives such as protein tyrosine phosphatase 1B (PTP1B) inhibitors were identified by molecular modeling techniques. Docking, molecular dynamics simulations and free energy calculations were employed to determine the modes of these new inhibitors. Binding free energies were calculated by involving different energy components using the Molecular Mechanics-Poisson-Boltzmann Surface Area and Generalized Born Surface Area methods. Relatively larger binding energies were obtained for the catechol derivatives compared to one of the PTP1B inhibitors already in use. The Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) free energy decomposition analysis indicated that the hydroxyl functional groups and biphenyl ring system had favorable interactions with Met258, Tyr46, Gln262 and Phe182 residues of PTP1B. The results of hydrogen bound analysis indicated that catechol derivatives, in addition to hydrogen bonding interactions, Val49, Ile219, Gln266, Asp181 and amino acid residues of PTP1B are responsible for governing the inhibitor potency of the compounds. The information generated from the present study should be useful for the design of more potent PTP1B inhibitors as anti-diabetic agents.

Introduction

Protein tyrosine phosphatase is an exclusive phosphatase that appears to be associated with the regulation of several growth factor signalling pathways. Biochemical and genetic studies have established that PTP1B plays a vital role in regulating body weight, glucose homeostasis, and energy expenditure by acting as a key negative regulator of insulin receptors and leptin receptors mediated signalling pathways. Several experiments have also demonstrated that the knock-out of PTP1B in mice can result in insulin hypersensitivity, even in a high-fat diet.1-6 Therefore; small molecular inhibitors of PTP1B are promising drugs for curing type II diabetes. The full length of PTP1B consists of 435 amino acids, con-stituting the major cellular form; however, only a shorter length of 298 or 321 residues is typically considered in biochemical studies. The model used in this study is the one coded 1wax.pdb in the protein Data Bank which contains 298 amino acids. The active site of PTP1B consists of residues His214-Arg221 and loops WPD (Thr177-Pro185), R (Val113-Ser118). The other important parts of PTP1B are S loops (Ser201-Gly209), R3-helix (Glu186-Glu200), R6-helix (Ala264-Ile281), and R7-helix, (Val287-Ser295) which take part in the catalysis substrate binding.7-13 The second binding active site close to the conserved primary active site (also referred to as site B), which was determined by Zhang et al., is not so conserved and therefore it is believed that it should be exploited in the design of PTP1B enzyme inhibitors with good selectivity.14

Experimental studies have shown that catechol demon-strates an inhibition activity against the PTP1B enzyme. However, the binding mode and inhibition effect of this type of inhibitor against the enzyme has not been well-document-ed. The aim of this paper is to design some new catechol derivatives and to determine their actions as inhibitors against the PTP1B enzyme and also to understand the main driving forces behind the interaction between the inhibitors and the enzyme by computational modeling.

 

Methodology

All molecular dynamic calculations were performed using the Assisted Model Building with Energy Refinement (AMBER) suite of programs (version 11).15-17 3D structures were displayed using Chimera (UCSF),18 VMD19 and DSV.20 RMSD graphics were shown by the XMGRACE package program. Xleap as implemented in AMBER was employed to prepare parameter/topology and coordinate the files and solvate and also to neutralize the system for the MD simula-tions. Structures were solvated with a TIP3P21 water model by creating an isometric water box where the distance of the box is 10 Å from the periphery of the protein-ligand com-plex. An ff99SB22 force field was used for the protein. Atomic partial charges were determined by the antechamber module of the AMBER package using AM1-Bcc (Austian model with Bond and charge correction)23 for the ligands and the General AMBER Force Field (GAFF)24 was adopted in simulation for the ligands because it handles small organic molecules. Hydrogen bond analyses between the protein and the ligands and RMSD changes with time during MD simulations were calculated by ptraj module as implemented in the AMBER programs package.

The crystal structure of the protein was obtained from the Protein Data Bank (1wax.pdb)25 and it was used as the protein model. Crystallographic water molecules were remov-ed from 1wax prior to use. All the ligands were designed by DSV, followed by optimization. Minimizations were perform-ed using the AMBER v11 program package with the parm99SB parameter set for protein atoms and the GAFF force field for ligands. AM1-Bcc charges (Austin model 1-bond charge correction) were computed for the ligands using the ANTECHAMBER module within the AMBER package. 1 ns MD simulation was performed for each ligand in a vacuum, at 300 K.

Initial complexes were located in a cubic box of explicit TIP3P water molecules with a maximum distance between the protein and the edge of the box of 10 Å. Finally, the necessary counterions were added in order to neutralize the systems. Periodic boundary conditions were used in all simulations with the Particle Mesh Ewald method26 to com-pute long electrostatic interactions. A cut off distance of 10 Å was chosen to compute van der Waals (VdW) non-bonded interactions.

Molecular Docking. Dock 6.0 module27 allows the performing of all stages of a docking process, generation of ligand conformations, ligand docking, and scoring of the binding modes. As in this case, where a rigid receptor ap-proximation was used, it is expected that the different receptors considered will lead to different ligand-binding modes depending on the initial size of the PTP1B-binding cavity. Thus, the four new PTP1B inhibitors were docked on the available receptor following a multi step procedure. In order to describe receptor-binding properties, a grid of potential energy was calculated for atoms taking part in the binding pocket. These atoms were obtained from the analysis of each protein–ligand complex. In this step, default para-meters were used. Then, the ligand was docked using the calculated grid to place it into the cavity and score the pro-posed binding modes.

Molecular Dynamics Simulations. Coordinates of the protein mentioned above were taken from the Protein Data Bank (1wax PDB). Crystallographic water molecules were removed from all the structures. In 1wax, the missing coordinates of the atoms were modeled using XLEAP and an ff99SB force field. Atoms on PTP1B were assigned the PARM99 charges, and all ionizable residues were set at their default protonation states at neutral pH. All structures were further processed by the XLEAP module of AMBER. Struc-tures were solvated with a TIP3P water model by creating an isometric water box, where the distance of the box is 10 Å from the periphery of protein. The molecular systems were neutralized by the addition of counterions. The systems and then the energy were minimized in two steps; in the first step, the protein and ligand were kept fixed, only the water molecules were allowed to move, and in the second step, all atoms were allowed to move. For the first step, the energy minimization was performed in 500 and 2500 steps with the steepest descent and conjugate gradient methods, respective-ly. For the second step, the energy minimizations were performed in 500 and 2500 steps using the steepest descent and conjugant gradient methods, respectively. Heating was performed with an NVT ensemble for 200 ps where the protein-ligand complex was restrained with a force constant of 10 kcal/mol/Å. Equilibration was performed for 200 ps on an NPT ensemble restraining the protein-ligand complex by 1 kcal/mol/Å2. Final simulations, the production phase, were performed for 10 ns on an NPT ensemble at a 300 K temperature and 1 atm pressure. Step size was 2 fs for the entire simulation. A Langevin thermostat and barostat were used for coupling the temperature and pressure. A SHAKE algorithm was applied to constrain all bonds containing hydrogen atoms. The nonbonded cut off was kept at 10 Å, and long range electrostatic interactions were treated by the particle mesh Ewald (PME) method with a fast Fourier transform grid with approximately 0.1 nm space. Trajectory snapshots, which were finally used for analysis, were taken at each 1 ps.

MM-PB/SA. Free energy calculations have been widely explored by Molecular Mechanics/Poisson Boltzmann Surface Area (MM/PBSA) and Molecular Mechanics/Generalized Born Surface Area (MM/GBSA), which combine molecular mechanics energy and implicit solvation models.28-30 They are more computationally efficient compared with the rigor-ous methods such as free energy perturbation (FEP) and thermodynamic integration (TI) methods.31 Another related method is the linear interaction energy (LIE), which aver-ages interaction energy from the MD simulations to estimate the absolute binding free energy. This approach restricts the simulations to the two end points of ligand binding, which is similar to MM/PBSA and MM/GBSA. MM/PBSA and MM/GBSA do not use a large training set to fit different parameter for each energy term32-37 which are different from most empirical scoring functions employed in molecular docking and both also allow for rigorous free energy de-composition into contributions originating from different groups of atoms or types of interaction.38-40

The MM-PBSA methodology is held to be one of the more computationally tractable means of obtaining reason-able estimates of the free energy of a complex system. In essence, it is quite straightforward: one performs a conven-tional molecular dynamics simulation of the complex in a periodic water box with counterions and the resulting tra-jectory is then post processed by removing the solvent and the periodicity and calculating the average of free energy over a series of static frames or “snapshots” according to the formula below.

G = EMM + Gpolar + Gnonpolar − TSMM

Here, EMM is the average sum of molecular mechanical energy terms, Ebond + Eangle + Etorsion + EVdW + Eelectrostatic. Gpolar and Gnonpolar describe the free energy of the solvent continuum. The Gpolar term may be obtained either through the solution of the Poisson–Boltzmann equation, or using an equivalent Generalized Born approximation. The Gnonpolar part is usually obtained by scaling the solvent accessible surface area by an appropriate surface tension. TS is the entropy of the solute either from quasi-harmonic analysis of the trajectory or from normal mode calculations on a (limited) number of snapshots. Then binding free energies are obtained from

ΔGbinding = Gcomplex − Greceptor − Gligand

All energy components were calculated using 100 snapshots from 5 ns to 10 ns. The snapshots from 5 to 10 ns molecular dynamics (MD) simulation trajectories were taken from the calculation of free energy. For PBSA and GBSA calcula-tions, dielectric constants solvent was taken as 80.0, respectively.41,42

Free Energy Decomposition Analysis. Free energy was decomposed to estimate the contribution of each residue in the binding process and was performed using MM/PBSA. Free energy was calculated by the MM/GBSA method. The energy of each residue-inhibitor interaction is given by the following equation:

ΔGGB is the free energy due to the solvation process of polar contribution calculated using the generalized Born model. MD simulation trajectories in the 10 ns were taken, and all energies were calculated for each frame. Total energy, the average energy of backbone and side chains for each residue was calculated. In order to get the contribution of each residue to the total binding energy, free energy decomposition is calculated.

ΔGGB is the free energy due to the solvation process of polar contribution calculated using the generalized Born model, ΔGSA is free energy due to the solvation process of nonpolar contribution and was calculated from SASA. MD simulation trajectories at the range of 5-10 ns were taken, and all the energies were calculated for each frame. The average energy of backbone and side chains for each residue was separately calculated, and the total energy was calculat-ed as well.

 

Results and Discussion

This work aimed to investigate the binding affinity of some catechol compounds against the PTP1B second active site. In the literature, different catechol derivatives showed a significant inhibitory effect on protein. Targeting this site might be an alternative challenge in the development of drugs in the fight against diabetes and obesity.

The distribution coefficient of a drug largely influences how easily it can reach its projected target in the body, how strong a consequence it will have once it arrives at its target, and how long it will remain in the body in active form. The distribution coefficient of a drug is measured by LogP, which is one criterion used in medicinal chemistry to evaluate the drug likeness of a given molecule, and to estimate lipophilic efficiency, a function of potency. LogP also assesses the quality of research compounds. For a given compound, lipo-philic efficiency is defined as the pIC50 (or pEC50) of interest minus the LogP of the compound.43,44 LogP and lipophilic efficiency calculated for the compounds are given in Table 1.

All the ligands were docked into the active site of the enzyme as displayed in Figure 2. They produced docking scores ranging from −28 to −39 kcal mol−1 as shown in Table 1.

MD calculations were performed to explore the dynamic stability of the complexes of PTP1B with the ligands. The calculations showed that all the ligands strongly maintained their position over a period of 10 ns. The potential energies as a function of time during the simulation for each complex are shown in Figure 3.

Table 1.*kcal mol−1

Figure 1.The 3-D structures Catechol derivative inhibitors.

Figure 2.Ribbon diagrams of the complex structures PTP1B (1wax). The inhibitor molecules with the best docked energy is shown as a ball and stick model.

To ensure the rationality of the sampling strategy, the root-mean- square-displacement (RMSD) values of the protein backbone atoms (N, CA, C) were calculated based on the starting snapshot as plotted in Figure 4. The RMSD plot in Figure 4 indicates that the conformations of PTP1B/L1 complex achieve equilibrium around 500 ps and fluctuate around 1.5 Å while for the PTP1B/L2 complex, the equilibrium time is around 400-500 ps and the conformations fluctuate around 1.5-2.5 Å. Both trajectories are stable after 600 ps, so it is reasonable to perform the binding free energy calcula-tion and free energy decomposition based on the snapshots extracted from 4 to 10 ns. As to the complex of PTP1B/L3, the equilibrium time is around 400 ps and the conformations fluctuate around 1.5-2 Å while for the complex of PTP1B/L4, the equilibrium time is around 300-400 ps and the conformations fluctuate around 1.5-2 Å.

A more detailed analysis of the root-mean-square fluctua-tion (RMSF) versus the protein residue number for all the complexes is illustrated in Figure 5. In this figure, it can be seen that four inhibitor/protein complexes had similar RMSF distributions, indicating that all inhibitors could have similar interaction modes as PTP1B on the whole.

Figure 3.Black line represents kinetic energy and red line represents potential energy for the 10 ns MD simulation of PTP1B-inhibitor complexes.

In the MD simulations, we found that all bound complexes undergo various degrees of conformation and orientation adjustments. Compared to the energies obtained by the docking procedure, the binding energies ΔGMM-PBSA calculated by the MM-PB/SA method with trajectories from MD simula-tions are significantly compatible (Table 2). This means that the docking tends to overestimate binding affinities, like many other docking programs. The energy contributions to the binding energy are calculated as ΔEvdw = −10.87 and ΔEelec = −17.48 kcal mol−1 for the complex of 1, ΔEvdw = −15.14 and ΔEelec = −27.44 kcal mol−1 for the complex of 2, ΔEvdw = −23.44 and ΔEelec = −21.83 kcal mol−1 for the complex of 3 and ΔEvdw = −17.67 and ΔEelec = −27.06 kcal mol−1 for the complex of 4.

Figure 4.Root-mean square deviations (RMSD) of the backbone atoms (CA, N, C) of the complexes with respect to the first snapshot as a function of time

Figure 5.Root-mean square fluctuations (RMSF) of the backbone atoms (CA, N, C) versus residue number for (1-4 inhibitors) the PTP1B complex.

Table 2.ΔH, Enthalpy of binding calculated by MM/PBSA method in kcal mol−1; ΔEvdw, Van der Waals contribution; ΔEele elektrostatic contribution; ΔGPB, the polar contribution of desolvation; at 300 K, Kia calculated inhibition constant (ΔG = kTlnki), k = 1.987 kcalmol−1, T = 300 K, *% inhibition at 20 M or IC50 (M)

Hydrogen bond analyses between the enzyme and the ligands are listed in Table 3. They indicate that L1 has favorable hydrogen bonds with the residues of 48, Asp181 and Gln262, L2 forms hydrogen bonds with Tyr46 and Asp181, L3 forms hydrogen bonds with Tyr46 and Asp181, while L4 has hydrogen bonds with Asp181 and Gln262. It is obvious that all the ligands form hydrogen bonds through phenolic functions and it is evident that this molecular segment plays a crucial role in the binding and recognition.

For the purpose of obtaining a detailed presentation of the catechol derivatives and PTP1B interactions, free energy decomposition analysis was employed to decompose the total binding free energies into inhibitor-residue pairs. The quantitative information of each residue’s contribution is extremely useful to interpret the binding modes of compounds. Calculated decomposition energy for the complexes of catechol derivatives by the MM/GBSA method demon-strates that L1 in the active site of the protein interacts with Gln262, Val49 and electrostatic interactions with Asp181. L2 interacts with Tyr46, Val49 and Asp181 parts of the protein. L3 is involved in this type of interaction with Gln262, Ile218, Asp181, Arg257 and Met258, while for L4 a similar mode of interaction is observed between Val49 and Asp181 units of residues in the enzyme (Figure 6 and Table 4).

Figure 6.The orientation of the residues shown as sticks around inhibitor compounds in the final snophot of MD simulation.

Table 3.Visible percentage of hydrogen bonds during MD simulations between Human-PTP1B and catechol derivatives

Table 4.Calculated decomposition energy for complexation of catechol derivatives by MM/GBSA method

The 10 ns molecular dynamics (MD) simulations on complexed PTP1B were fulfilled with the aim of revealing the possible mechanisms of ligand recognition and inhibi-tion. Based on our dynamics simulation and decomposition analysis, many useful results were obtained. The most notable fact is the movement of the second active site, and the types of interactions associated with it. Second, according to our calculated Van der Waals and electrostatic energies between PTP1B and the inhibitor catechol derivatives, it is evident that the more contributions to the Van der Waals to PTP1B-inhibitor complexes conformation, binding site, the more electrostatic interaction results, which is consistent with other available experimental studies. Third, by analyzing the interactions between PTP1B and the high-affinity inhibitor compound 2, it was found that the residues adjacent to the second active site, including Asp181, Gln262, Ile219, Val49 and Phe182, might be partially responsible for the high inhibitory activity and selectivity of PTP1B. According to the MD simulations, MM/PBSA free energy calculations, and MM/GBSA free energy decomposition analyses, we can draw the following conclusions: Our simulation results sug-gest that the Sanggenon derivatives PTP1B inhibition, the effects of which are not mentioned in the literature, might be functioning as potent and selective PTP1B inhibitors. In the design of more effective medicine the development of such molecules and targeting the surface residues, for example, the region containing Met258, Arg254, Gln102 and Asp29 of the second phosphate binding site might be advantageous.

Figure 7.The structures were superimposed on the PTP1B for the comparison of the binding mode. Ligand 1 is colored in pink, Ligand 2 is in Blue, Ligand 3 is in yellow and Ligand 4 is in red, respectively.

In conclusion, ligands L2, L3 and L4, being potent inhibitors of PTP1B enzyme, were selectively docked into the active site of the protein. Docking and binding energy calculations indicate that all the ligands have favorable binding with the enzyme. In order to compare the location of individual ligands and to see their binding modes, they were superimposed as shown in Figure 7. This theoretical out-come may provide basic information, which might be further used in the development of anti-diabetic drugs treating large numbers of the world population.

 

Conclusion

The 10 ns molecular dynamics (MD) simulations on com-plexed PTP1B were carried out with the aim of revealing the possible mechanisms of ligand recognition and inhibition. The binding affinity was investigated using MD simulations, MM-PB/SA free energy calculations and MM-GB/SA free energy decomposition analysis. Based on our dynamics simulation and conformation analysis, many useful results were obtained. As a sum of this study, the available calcu-lated evidence shows that the molecules derived from the catechol molecule, which is known as an effective molecule, can hang up the related enzyme more easily. It is especially remarkable that the compound 4 is tied with the lowest energy to the active area of protein and its stability in that area is also remarkable. For all the trained models, the connection areas of protein are correlative. In particular, the activity of amino acids, which are: Val49, Ile218, Gln262, Tyr46, Asp181, Met258, Arg257 and Phe 182, was deter-mined with the calculations of decomposition. Our simula-tion results suggest that potent and selective PTP1B inhibitors may be designed by targeting the surface residues.