Structural Basis of Fullerene Derivatives as Novel Potent Inhibitors of

Sep 20, 2016 - The corresponding structures are represented as cartoon (light pink), and the WPD loop is highlighted in dark blue. The transition of t...
0 downloads 0 Views 6MB Size
Subscriber access provided by CORNELL UNIVERSITY LIBRARY

Article

Structural Basis of Fullerene Derivatives as Novel Potent Inhibitors of Protein Tyrosine Phosphatase 1B: Insight into the Inhibitory Mechanism through Molecular Modeling Studies Mengdan Qian, Yaming Shan, Shanshan Guan, Hao Zhang, Song Wang, and WeiWei Han J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.6b00482 • Publication Date (Web): 20 Sep 2016 Downloaded from http://pubs.acs.org on September 24, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Information and Modeling is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Structural Basis of Fullerene Derivatives as Novel Potent Inhibitors of Protein Tyrosine Phosphatase 1B: Insight into the Inhibitory Mechanism through Molecular Modeling Studies

Mengdan Qiana, Yaming Shanb, Shanshan Guanb, Hao Zhang*a, Song Wang*a and Weiwei Han*bc

a

Institute of Theoretical Chemistry, Jilin University, Changchun 130023, China Key Laboratory for Molecular Enzymology and Engineering of Ministry of Education, National Engineering Laboratory for AIDS Vaccine, School of Life Science, Jilin University, Changchun 130012, China c Department of Computer Science, C.S. Bond Life Sciences Center, University of Missouri Columbia, Missouri 65211, USA b

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract: Protein tyrosine phosphatase 1B (PTP1B) has become an outstanding target for the treatment of diabetes and obesity. Recent, research has demonstrated that some fullerene derivatives serve as a new nanoscale-class of potent inhibitors of PTP1B, but the specific mechanism remains unclear. Several molecular modeling methods (molecular docking, molecular dynamics simulations, and molecular mechanic/generalized Born surface area calculations) were integrated to provide insight into the binding mode and inhibitory mechanism of the new class fullerene inhibitors. The results revealed that PTP1B with open WPD loop was more susceptible to the combination of fullerene inhibitor due to the comparable shape and size. When the WPD loop fluctuated to the open conformation, the inhibitor fell into the active pocket and induced conformational rotation of the WPD loop. This rotation was closely related to the reduction of the catalytic activity of PTP1B. In addition, it was suggested that compound 1, the same as compound 2, is a competitive inhibitor since it blocked the active site to prevent the binding of the substrate. The high binding affinity of fullerene-based compounds and the transition of WPD loop, caused by the specific structural property of hydrophobic fullerene core and the appending polar groups, make these fullerene derivatives efficient competitive inhibitors. The theoretical results provide useful clues for further investigation of the noval inhibitors of PTP1B at the nanoscale.

Keywords: Protein tyrosine phosphatase 1B (PTP1B), fullerene derivative inhibitors, molecular docking, molecular dynamics simulation, inhibitory mechanism .

ACS Paragon Plus Environment

Page 2 of 35

Page 3 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

1. INTRODUCTION More than 130 million people worldwide are afflicted with type II diabetes, and the incidence is expected to grow steadily over the next several years. The current treatments for type II diabetes include injectable insulin and oral hypoglycemic agents, such as metformin, sulfonlyureas, thiazolidinediones, and R-glucosidase inhibitors. A major goal of new therapies for type II diabetes is to potentiate the action of insulin. The protein tyrosine phosphatase 1B (PTP1B), known as the hydrolase of phosphotyrosine (pTyr)-containing proteins, functions as a major negative insulin signal transduction regulator by dephosphorylation of the insulin receptor.1-3 Inhibition of PTP1B can improve insulin sensitivity and result in the resistance to obesity.4-6 As a result, PTP1B has become an attractive target for the treatment of type II diabetes and obesity, and selective inhibitors of PTP1B could be of significant therapeutic utility.7-11 As a typical intercellular protein tyrosine phosphatase, PTP1B is composed of a single catalytic domain and various functional loops to regulate the activity of proteins.12 The overall structure of PTP1B is shown in Fig. 1a, and the important regions are marked in different colors, including the P-loop (His214 to Arg221), the substrate recognition loop (Tyr46, Val49 and Lys120), Q-loop (Gln262), the secondary binding site (Try20, Arg24, His25, Ala27, Phe52, Arg254, Met258, Gly259) and the WPD loop (Thr177 to Pro185).13-15 These functional regions, which are responsible for substrate catalysis and recognition, are located around the active site, forming an active pocket (Fig. 1b). The active pocket is positively charged because of the adjacent arginine and lysine residues; thus, the negatively charged phosphate moiety is easily attracted into the active site.16 The catalytic mechanism of PTPs has been thoroughly studied through numerous X-ray crystallographic studies, in which the featured conserved peptide chain Cys-(Xaa)5-Arg of the P-loop was identified to play key roles both in the stabilization and recognition of the bound phosphate substrate.17, 18 Before hydrolysis, Arg211 firstly forms a hex-coordinated intermediate with the phosphoryl substrate by hydrogen bonding to immobilize phosphate in the active pocket. Subsequently, residue Cys215 functions as a nucleophile to attack the phosphoryl moiety with the concerted assistance of general acid Asp181, leading to the cleavage of the P-O bond of the pTyr group.19, 20 Finally, the active Cys215 is regenerated and the intermediate breaks down under the attack of the water molecule activated by Gln262 as well as Asp181 (Fig. S1). Remarkably, PTP1B exhibits an equilibrium mixture of two distinct modes, the open and closed WPD conformations, in ligand-free solution.21 Fig. 1c is an overlay of the two PTP1B structures with different WPD loop modes. Residues Asp181 and Phe182 on the closed WPD loop are shifted down towards the active site compared to the open WPD loop (Fig. 1d). This conformational flexibility of the WPD loop is essential to substrate recruitment, solvent accessibility and product release. The WPD loop is of vital importance in phosphoester hydrolysis and the closure of

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the WPD loop enables the acid to approach the substrate, and allowing dephosphorylation to proceed successfully.22 Taken the characteristic property of the WPD loop, it is expected that active site with open WPD loop might be a potential target for the binding of new inhibitors.23, 24

Fig. 1 The overall structure of PTP1B with its important regions highlighted in different colors (a) and the key residues on these regions (b). The superimposition of PTP1Bs with open and closed WPD loop is shown in (c) and comparisons of residues in two different modes are depicted in (d). Carbon nanoparticles (CNPs), such as nanotubes and fullerenes, have been found to interact with and even influence protein assemblies.25 CNPs have attracted great attention due to their special physicochemical properties, such as large surface area and high penetration ability. Engineered or crafted fullerene derivatives with well-processed shapes and chemical properties have wide applications, ranging from medical therapy and drug delivery vectors to nanomaterial science.26, 27 Several fullerene derivatives have recently been synthesized and proposed as a novel class of PTP inhibitors. According to the experiment, the engineered fullerene derivatives, shown in Fig. 2 (supporting information), exhibit enhanced and potent inhibition compared with some other small molecule inhibitors due to their comparable size to accommodate the active pocket perfectly.28 Moreover, the fullerene derivatives exhibit lower toxicity and are safer in medicinal applications with decorated hydroxyl and carboxyl groups on the carbon cage. Compared with single-branched compound 2 (IC50 = 0.20 µM), compound 1 (IC50 = 0.08 µM), with polyhydroxyl groups, shows almost 2.5-fold inhibitory activity over PTP1B, indicating that the two inhibitors exhibit moderate selectivity over PTP1B.28 The experimental data suggest that compound 2 is likely to be a competitive inhibitor because its Lineweaver-Burk plot corresponds well with the competitive type. However, the inhibitory mechanism and interactions between the two inhibitors and PTP1B remain unknown. In this study, molecular docking, molecular dynamics simulation and molecular mechanic/generalized Born surface area calculation were employed to study the inhibitory mechanism responsible for the high potency of the fullerene inhibitors. It is promising that the investigation of

ACS Paragon Plus Environment

Page 4 of 35

Page 5 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

atomic interactions between the fullerene inhibitors and PTP1B can provide innovate insight into the design of a new class of inhibitors at the nanoscale.

Fig. 2 Chemical structures of two fullerene derivatives.

2. METHODS 2.1 Constructions of PTP1B and fullerene derivatives Since the X-ray crystallographic structure has validated that there exist two different modes of PTP1B (WPD loop open and closed), it is essential to take every structure into consideration. The two X-ray structures were obtained from RCSB protein data bank (www.rcsb.org, PDB code: 1NL9 and 2CM8 for open and closed mode respectively). Subsequently, the crystallographic waters, ions and the ligand were removed to get the pure protein. The 3D structures of two inhibitors were firstly built up by Gaussview 5.0. Then the geometry optimization of inhibitor was performed by Gaussian 09 with ab initio calculation method at the HF/6-31G (d) level. 2.2 Molecular docking calculation The AutoDock Tools were used for the setup of docking running and the molecular docking was simulated by AutoDock Vina software. AutoDock Vina is a new molecular docking program using the elaborate optimization algorithm and method to complete the docking as well as the scoring function.29 The definition of the binding site is described by a rectangular box with explicit coordinates. We adjusted the box to the center of primary and secondary binding sites. Both the open and closed PTP1B proteins were simulated to evaluate which conformation was more suitable for the fullerene binding. The three dimensional grids were set as 46× 40×36 and 40 × 40 × 24 points with a grid spacing of 1 Å for the open and closed mode, respectively. The center of grid was regarded as the mass center and the exhaustiveness was 20. Finally, the fullerene inhibitors were docked to the two PTP1B proteins by turn and the corresponding energy evaluations were also generated. The most popular and stable docking conformation was selected as the initial binding mode for the subsequent molecular dynamics simulation. 2.3 Molecular dynamics simulation The molecular dynamics simulation was performed using the GROMACS 4.5.2 software package30 on the basis of the predicted complexes generated by Autodock

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Vina. The GROMOS96-53A6 force field was applied to build the parameter files of the protein, while corresponding topology files of the ligands were generated by the Automated Topology Builder (ATB) server with the charge distribution calculated by the HF/STO-3G level quantum method.31, 32 The whole complex was solvated with the SPC water model in period boundary condition 10 Ǻ from all solute atoms.33 The proper amount of sodium ions was then added to the solvent in order to neutralize the charge of the whole system. Prior to the MD simulation, the energy of the system was minimized with the steepest descent algorithm to prevent steric clashes or improper geometry. After the minimization, a 500 ps positional restrained isochoric-isothermal (NVT) simulation and a 500 ps isothermal-isobaric (NPT) simulation were run sequentially to relax the system. Finally, the MD simulation was performed with a time step of 2 fs when the system was under constant temperature (300 K) and constant pressure (1.0 bar) conditions. The reference pressure and temperature were coupled with Parrinello-Rahman and V-rescale, respectively.34 During the simulation, the LINCs algorithm was employed to fix all bonds and the long-range electrostatics was computed by the Particle Mesh Ewald method with a grid spacing of 0.16 nm.35 2.4 MM/GBSA binding free energy calculation and the ligand-residue energy decompositions analysis In this study, the binding free energy was calculated with the procedure of molecular mechanic/generalized Born surface area (MM/GBSA) supported by the AMBER 12 software package.36-38 The calculation was performed on the basis of following equations: ∆Gbinding = Gcomplex− [ Gprotein + Gligand ] =∆H− T∆S (1) and ∆Gbinding = ∆EMM + ∆GGB + ∆GSA− T∆S (2), where Gbinding is the binding free energy and Gcomplex, Gprotein, and Gligand are the free energies of the complex, protein and ligand, respectively. ∆H is the enthalpy change, which includes three components: the molecular mechanics interaction energy between the protein and ligand (∆EMM), the polar solvent energy contribution (∆GGB) and the nonpolar solvent energy contribution (∆GSA) to the desolvation of ligand binding.39 The generalized Born approximation model was employed to compute the polar solvent free energy. The part of −T∆S represents the change in conformational entropy.40 Here, the contribution of the entropy term was omitted because the Normal-mode analysis for the entropy calculation had high computational demands but did not produce a sufficiently accurate result. A total of 100 snapshots were extracted from the MD trajectory to calculate the protein-ligand binding free energy.41 The free energies of each protein residue with the ligand were also calculated by AMBER 12. The binding interaction of each residue-ligand pair consists of three terms: the van der Waals contribution (∆Eval), the electrostatic contribution (∆Eele), and the solvation contribution (∆GGB + ∆GSA). ∆Gresidue-ligand = ∆Evdw + ∆Eele +∆GGB + ∆GSA (3) Same Snapshots extracted from the MD simulation were used to calculate all the energy components of the free energy decomposition.42 2.5 Principle component analysis and cross-correlation map analysis

ACS Paragon Plus Environment

Page 6 of 35

Page 7 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Principle component analysis (PCA) has developed to be an effective and useful tool to identify large-scale motion and the correlated movements of macromolecular biological systems.43 PCA of the covariance matrix is becoming popular for its wide application to protein systems by reducing or simplifying large and complicated data sets along the trajectory generated by molecular dynamics simulation.44-46 The structures during the MD simulation are extracted to construct the domain cross-correlation matrix map and the motions responsible for the most variance in atomic position are targeted when diagonalizing the covariance matrix of the atomic coordinates of the system. The covariance matrix over the collected structures can be calculated on the basis of the equation: Cij = 〈(χi − 〈χi〉)(χj − 〈χj〉)〉 (i, j = 1, 2, . . . , 3N) (1) where χi (χj) is the Cartesian coordinate of the ith (jth) atom and 〈χi〉 (〈χj〉) is the time average position of the ith (jth) atom computed throughout the simulated trajectory. The PCA calculation is implemented on the coordinates of the Cα atoms of the protein structure throughout the entire simulation time. The data can be visualized by means of correlation maps calculation based on the equation: Cij = 〈(χi − 〈χi〉)(χj − 〈χj〉)〉/(χi − 〈χi〉)1/2(χj − 〈χj〉)1/2 (i, j = 1, 2, . . . , 3N) (2) The covariance matrix map allows the identification of residue pairs with correlated and anticorrelated motion.47 i (j) is the ith (jth) residue and N is the total number of Cα atoms. The projections of the original structures onto the eigenvectors of the related principal components are represented as plots of the cross-correlation map.48, 49

3. RESULTS AND DISCUSSION 3.1 Determination of the optimal initial binding mode of fullerene inhibitors by docking analysis The molecular docking simulation was executed to provide a relatively reliable initial binding model of the inhibitor. Considering that there exist two modes of PTP1B, the inhibitors are respectively docked to the targeted active site of different receptors (open and closed mode) to identify the preferable binding conformation. The optimal docked conformation is selected from the 20 top poses of inhibitor generated from the Autodock Vina and the corresponding binding free energies are listed in Table 1. Generally, the inhibitors have lower binding free energy with open-WPD PTP1B than the corresponding closed WPD, indicating that inhibitors are more inclined to bind to the open active pocket which can be attributed to the comparable size to the fullerene derivatives. The optimal docked models of the two inhibitors with the different PTP1B are shown in Fig. 3. In terms of the open-PTP1B, both of the two inhibitors fall into the open active pocket and form hydrogen bonds with the polar residues Arg221, Lys116 and Gln266. The hydrophobic interaction from Tyr46 is also non-negligible in the stabilization of the fullerene inhibitors (Fig. 3a and 3c). The main difference is that inhibitor 1 is embedded onto the groove of the active pocket due to the bigger size, whereas inhibitor 2 can extend deep into the active site through the only modified hydrophilic branch. In contrast, as for the closed-mode PTP1B, both inhibitors interact with residues in the vicinity of the secondary binding site (Gln21, Arg24 and

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Gln262) rather than the active site. This binding mode is not as stable as that in the open mode because the closure of the WPD loop blocks the binding of inhibitor to the active site and thus they are just adsorbed on the surface of enzyme. Moreover, the closed loop provides a smaller space for the inhibitor binding and the possibility of clashes is bigger, which might be one of the reasons for the low score of docking evaluation. Both the binding affinity analysis and the binding mode inspection suggest that PTP1B with and open WPD loop is more favorable for the combination of fullerene derivatives.

Fig. 3 The optimal docking models of the three complexes with key residues shown in magenta lines while the ligands are shown in cyan sticks. (a) compound 1 is docked to open-PTP1B, (b) compound 1 is docked to closed-PTP1B, (c) compound 2 is docked to open-PTP1B. Table 1 The docking affinities of the two inhibitors in the different modes of PTP1B. Complex Open-PTP1B/compound1 Closed-PTP1B/compound1 Open-PTP1B/compound2 Closed-PTP1B/compound2

Docking Affinity(Kcal/mol) −7.0 −5.7 −11.7 −9.1

3.2 Setup of simulated systems and evaluation of the structural equilibrium and flexibility Now that the experimental data indicate that compound 2 might be the competitive

ACS Paragon Plus Environment

Page 8 of 35

Page 9 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

inhibitor, it should occupy the active site to prevent the binding of the substrate. Therefore, the predicted complex with compound 2 docked to the active pocket of the open-PTP1B is selected as the initial model for the subsequent simulation. But the inhibitory type and mechanism of compound 1 have not been identified, suggesting that the binding site is uncertain between the active site (competitive type) and some allosteric sites (noncompetitive type). Therefore, two sets of parallel simulations were attempted on the basis of the two docked models of compound 1 with both open and closed PTP1B to investigate the inhibitory type of inhibitor 1. Taken together, three 100 ns MD simulations were performed on the following three complexes: open-PTP1B/compound1, open-PTP1B/compound2 and closed-PTP1B/compound1. The root mean square deviation (RMSD) of the Cα atoms with respect to the initial positions is usually regarded as an indicator of systematic stability. As shown in Fig. 4a, the RMSD curves of the backbone-Cα with respect to the three complexes achieve acceptable equilibrium after 50 ns with the average value of 0.271 nm, 0.236 nm and 0.338 nm for open-PTP1B/com1, open-PTP1B/com2, and closed-PTP1B/com1, respectively (Fig. 4c). The RMSD value of closed-PTP1B undergoes a large fluctuation at approximately 50 ns in the presence of compound 1, and the corresponding maximum RMSD value frequency is 0.41 nm whereas it is just 0.30 nm in the open-PTP1B/inhibitor1 complex (Fig. 4b). All the data suggest that the docked complex with the closed WPD loop is not the optimal binding mode, and the protein experiences significant conformational adjustments. The regional conformational fluctuations can be assessed by calculating the Cα root mean square fluctuation (RMSF) of each residue from the MD trajectory. The profiles of RMSF are shown in Fig. 5. Generally, the overall motion of PTP1B with the inhibitor is very similar to the free protein, except for some crucial regions, such as the P-loop, Q-loop and the recognition loop. In fact, the variation of the RMSF curves is closely related to the interactions between the inhibitor and PTP1B. Located at the bottom of the active site, the P-loop (His214-Arg221) and Q-loop (Gln262) should have exhibited some flexibility, but these two loops show almost no fluctuations due to the influence of the inhibitors (Fig. 5c-e), which indicates that these loops form interactions with the ligand and thus are restricted to a certain extent. Additionally, the WPD loop show basically identical flexibility in the three complexes with the naked PTP1B, suggesting that the inhibitors interact little with the WPD loop. Notably, the recognition loop (Tyr46-Val49 and Lys116-Lys120) is highly stable under the interactions of inhibitor 1 (Fig. 5c and 5d), while it exhibits amplified fluctuation in the presence of inhibitor 2 (Fig. 5e). It is obvious that the recognition loop has a greater effect on the binding of inhibitor 1.

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Fig. 4 (a) The RMSD profiles of backbone Cα atoms of three complexes as function of time. (b) Relative occurrence frequency of RMSD value of average structure. (c) Average RMSD value of the three systems during the 100 ns MD simulation.

Fig. 5 (a) and (b) are the RMSF curves of the Cα atoms of the three complexes. (b), (c) and (d) are the visualization of backbone flexibility for open-PTP1B/inhibitor1, closed-PTP1B/inhibitor1 and open-PTP1B/inhibitor2, respectively. The motion degree corresponds to the tube with different colors and thickness. 3.3 Structural motion after binding the fullerene inhibitors The correlation matrix analysis of PCA can promote the understanding of protein regions that have intense relevant conformational changes and can shed light on the dynamic motion of PTP1B induced by ligand binding. Thus, the correlation matrix for all pairs of Cα atoms was calculated. The number of eigenvalues used for the construction of the covariance matrix of open-WPD complex was 849, while 891 were used for the closed-WPD complex. The cumulative variance of the first five eigenvectors is shown in Table S1. It is found that the first 5 eigenvectors represent at least 40% of the total motions, and the first two represent almost 25%. Here, the covariance matrix map is constructed on the basis of the first two eigenvectors to depict the linear correlations of each pair of Cα atoms during the dynamics. The covariance

ACS Paragon Plus Environment

Page 10 of 35

Page 11 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

matrix maps of the three complexes are illustrated in Fig. 6, in which the antiharmonic and large-scale motions are highlighted by diagonalizing the matrix. The positive regions, marked in red and yellow, indicate the strong correlated motions of residues, whereas the negative regions, which are colored by blue and dark, are associated with the anticorrelated movements. The diagonal is relatively highly correlated since it represents the variance of a residue with itself. In general, motion is thought to be in the normal range when the values fluctuate between −0.2 and 0.2. Notably, compound 2 leads to strong centralized self-correlated motion (Fig. 6b) whereas compound 1 is able to impair this motion and gives rise to a series of new correlated and anticorrelated motion of PTP1B (Fig. 6a and c). The matrix maps of the three systems indicate that significant motions consistent or inconsistent, mainly occurs among the regions of residues Lys58-Asn64, Arg113-Trp125 (recognition loop), Tyr175-Pro185 (WPD loop), and Phe196-Pro210. Particularly, the motion between the recognition loop and the WPD loop deserves more attention since it occurs near the active site. The binding of inhibitor 1 largely decreases the correlated motion between the recognition loop and the WPD loop (box a3 and c3 in Fig. 6), whereas these motions are still intense in the presence of inhibitor 2 (box b3 in Fig. 6), suggesting the decreased flexibility of the recognition loop with inhibitor 1. In addition, when combined with inhibitor 1, the loop Lys58-Asn64 has strong fluctuations itself (box a1 and c1 in Fig. 6) and at the same time it exhibits weak opposite motion with the recognition loop Arg113-Trp125 (box a2 and c2 in Fig. 6). Conversely, the correlated motion of the Lys58-Asn64 loop becomes weak in the presence of inhibitor 2 (box b1 and b2 in Fig. 6).

Fig.6 The cross-correlation matrix map for the open-PTP1B/inhibitor1 complex (a), the open-PTP1B/inhibitor2 complex (b) and the closed-PTP1B/inhibitor1 complex (c). Although the generally correlated movements between the important loops are illustrated according to the cross-correlation map analysis, the specific motion trend of each region remains ambiguous. Therefore, the extreme projections for PC1 and PC2 are visualized by presenting the mode directions for each residue with arrows (Fig. 7), which describes the dominant motions during the 100 ns MD simulation. Observation of PC1 and PC2 for the three

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

complexes reveals that the loop of the secondary binding domain moves in nearly the same direction as the recognition loop, thus maintaining the integrity and compactness of PTP1B. In addition, the motion of the recognition loop is much weaker with presence of inhibitor 1 (Fig. 7a and Fig. 7c) compared with inhibitor 2 (Fig. 7b), which corresponds well with the flexibility analysis above. Although the WPD loop is not in parallel with the overall motion, it shows obvious and stable trends to get away from the active site and point to the α-helix connected to the secondary binding site. This inconsistent motion of the WPD loop might be related to the depressed catalytic activity of PTP1B.

Fig. 7 The motions of PTP1B on the basis of the first two PCs, and the arrow represents the directions and amplitudes of the movements. (a) is for open-PTP1B/inhibitor1 complex, (b) is closed-PTP1B/inhibitor1 complex, and (c) is open-PTP1B/inhibitor2 complex. In summary, the binding of inhibitor 1 impairs the compactness of the protein and further results in more scattered correlated and anticorrelated motions among important regions, while inhibitor 2 can induce more strong and centralized self-correlated motion in PTP1B. Moreover, the motion between the recognition loop and the secondary binding loop in the parallel direction may be responsible for the integrity of PTP1B. However, the nonparallel movements of the WPD loop can depress the catalytic process, which will be discussed in detail below. 3.4 Comparative binding modes and inhibitory mechanism of the two inhibitors In this section, structural interactions between PTP1B and the two inhibitors will be

ACS Paragon Plus Environment

Page 12 of 35

Page 13 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

discussed according to the MD simulation results as well as MM/GBSA calculations to provide comprehensive insights into the key residues and dominant interactions that play roles in the inhibitory action. It is very interesting to find that inhibitor 1 eventually evolves into a similar position despite starting from distinct initial docking models for the open and closed forms of PTP1B. The last snapshots corresponding to the open- and closed-mode complexes are respectively shown in Fig. 8a and 8d. As seen, the inhibitor has similar polar interactions with key residues of the recognition loop, the P-loop as well as the Q-loop. In the open-mode complex, inhibitor 1 makes small adjustments within the active pocket and finally forms stable hydrogen bonding interactions with residue Arg221 of the P-loop, Gln266 of the Q-loop, and Tyr46 and Lys120 of the recognition loop, with occurrence ratio reaching 87.44%, 98.08%, 77.24% and 61.06%, respectively (Table S2). This result is within expectations since P-loop is so important in the catalytic reaction that it is universally treated as the binding target for inhibitors. Correspondingly, in addition to these interactions, inhibitor 1 in closed-mode PTP1B has predominant polar interactions with residues Asn111 (37.52%), Lys116 (58.04%) and Lys120 (48.52%) of the recognition loop, which is much different from its initial docking mode. The trajectory analysis suggests that inhibitor 1 undergoes a relatively large translocation from the secondary binding site to the active pocket. The statistics on the total number of hydrogen bonds between the inhibitor and closed-PTP1B during the simulation suggest that the polar interactions are not as stable as the open-mode complex until 60 ns (Fig. S2). In fact, the curve of hydrogen bond number is consistent well with the RMSD profile and both profiles indicate that both the inhibitor and the closed-PTP1B experience conformational adjustments to accommodate each other. Observation of the MD trajectory suggests that the negatively charged inhibitor 1 is anchored into the active cavity by the strong electrostatic interactions of the positively charged residues (Lys116 and Lys120) on the recognition loop along with the open of WPD loop. In this manner, the inhibitor firstly inserts into the binding groove and subsequently makes adjustments to accommodate the various electrostatic interactions of the P-loop and Q-loop. Also, MD simulations of open- and closed-mode PTP1B with inhibitor 1 located at the outside of the enzyme (about 2 nm away from the active site) show that the inhibitor gets close to the active centre under the electrostatic interactions from the recognition loop (supporting information). Thus, the recognition loop is of great importance in the inhibitor recruitment to the active site PTP1B when the WPD loop is in the open mode. The MM/GBSA calculation also validates the significant electrostatic contribution of residues Lys116, L ys120, Ser216 and Arg221 in the presence of inhibitor 1 (Fig. 8c and 8f). Apart from the polar interaction, hydrophobic interaction is another crucial factor responsible for the high binding affinity of fullerene inhibitor (Fig. 8b and 8e). Residues Ala217, Ile219, and Gly220 of the P-loop, Tyr46, Val49, Asn111, and Glu115 of the recognition loop, Gln262, Thr263, and Gln266 of

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the Q-loop and Trp179 of the WPD loop create a hydrophobic pocket to immobilize the inhibitor. In case of the open-PTP1B/inhibitor2 complex, the prominent contribution to ligand stabilization is hydrophobic interaction, which is different to inhibitor 1. Both the fullerene cage and the pyrrolidine ring make numerous hydrophobic contacts with the side chains of Val49, Asn111, Cys215, Ser216, Ala217, Ile219, Gln262 and Thr263 (Fig. 8h). Specifically, the aromatic chains of Phe182, Trp179 and Tyr46 form stacking interactions with the inhibitor. The hydrophobic characteristics of inhibitor 2 make van der Waals contact the leading reason for the high binding affinity, as verified by the remarkable van der Waals contributions of the decomposition of the binding energy on each residue (Fig. 8i). By comparison, polar interaction, though weak, is induced by residues Ser216 (39.57%), Ala217 (44.83%), Arg219 (34.48%) and Arg221 (5.45%) of the P-loop. The hydrophilic branch of inhibitor 2 penetrates into the active pocket and establishes several hydrogen bonds with these polar residues (Fig. 8g). And Lys120 and Tyr46 of the recognition loop also make small contributions to the stability of inhibitor.

Fig. 8 The binding modes of the two inhibitors to the different modes PTP1B and the important residues that contribute to the stability of inhibitors. The first column shows key residues that make hydrogen bonding with the inhibitor, the second column is hydrophobic interaction, and the third column is the decomposition energy of each important residue. The first row is the open-PTP1B/inhibitor1 complex, the second row is closed-PTP1B/inhibitor1 complex and the third row is open-PTP1B/inhibitor2

ACS Paragon Plus Environment

Page 14 of 35

Page 15 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

complex. In order to measure the binding affinity of the inhibitor, the total binding free energies of the three complexes are calculated by MM/GBSA. The results are listed in Table 2. The total ∆Gbinding of inhibitor 1 is −44.63 kcal mol−1 and −42.18 kcal mol−1 for the open and closed PTP1B respectively (IC50=0.08), slightly lower than the inhibitor 2 complex at −40.74 kcal mol−1 (IC50=0.20). The order of magnitude of the three complexes is identical, indicating that two compounds exhibit modest inhibitory selectivity, which is in agreement with the experiment. Even so, the binding energy of compound 1 is a bit lower than compound 2 and this is also consistent with the IC50 value in the experiment. It is noteworthy that the van der Waals contributions are very similar in the binding energy of three complexes due to the hydrophobic contacts of the fullerene cage with the nonpolar residues around the active site. The electrostatic contribution is significant for inhibitor 1 with the electrostatic energy of −119.15 kcal mol−1 and −90.23 kcal mol−1 for the open and closed mode PTP1B respectively, whereas it is only −8.47 kcal mol−1 for the inhibitor 2 complex. However, the solvent energy in PTP1B/inhibitor1 complex increases correspondingly (134.36 kcal mol−1 and 108.22 kcal mol−1 for the open and closed modes respectively), which largely cuts down the total contribution. To sum up, the high binding affinity of inhibitors is attributed to the significant val der Waals contribution between fullerene cage and the nonpolar surface of amino acid residues at active site as well as the hydrogen bonding interactions of the decorated groups with the active site. The high binding affinity should be one primary reason for the efficient inhibitory potency of fullerene-based inhibitors. For another, the active is precisely blocked by the fullerene cage, which hampers the entrance of substrate. Similar interactions are observed at the fullerene-based inhibitory mechanism in terms of other enzymes,50, 51 indicating the similar mechanism of fullerene derivatives caused by the special structural features. Remarkably, conformational transitions of the WPD loop are observed during the simulation of the three complexes corresponding to the PCA analysis above. The rotation of the WPD loop can be quantified by the RMSD between the target structure and the initial structure of this loop (Thr177 to Pro185). As shown in Fig. 9, the movement tendency of WPD loop in the open-PTP1B/inhibitor1 and open-PTP1B/inhibitor2 complexes is similar (Fig. 9a and Fig. 9e), and both experience transient rotation from 30 ns to 40 ns before make irreversible deflection from the initial position. In comparison, the transition of the WPD loop in the closed-mode complex does not occur until 80 ns, corresponding to the rising stage of the RMSD curve (Fig. 9c). Trajectory analysis indicates that the closed WPD loop firstly fluctuates to the open state to allow the binding of inhibitors; then, the ligand adjusts to accommodate the active pocket, which impels the rotation of open-mode WPD loop. The transition of the WPD loop can be clearly observed by overlapping the developed complex with the initial PTP1B (Fig. 9b, Fig. 9d and Fig. 9f). The binding of inhibitors causes some of aromatic residues on the open-mode WPD loop, such as Pro180 and Phe182, to be exposed to the solvent. This

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

unfavorable state can be released by the deflection of aromatic residues to form stacking interactions with the fullerene backbone or turn to the hydrophobic cavity formed by Trp179, Phe191 and Phe269. The rotation of Pro180 further triggers the movement of the opposite-site residue Pro185, leading to the deviation of the whole WPD loop. However, residue Trp179 remains almost constant during the whole simulation since it forms strong cation-π interactions with Arg221. The calculation of the cation-π interactions of PTP1B shows that the residue pair Trp179-Arg221 makes a significant cation-π contribution, where Eelec is −4.5 kcal mol−1 (Table. S4) in the three developed complexes, and only −2 ~ −3 kcal mol−1 in the crystal structures (Table. S3). The overall effect of the rotation is that the hydrophobic residues Pro180, Phe182 and Pro185 are buried inside the protein while polar residues Asp181 and Glu186 are exposed to the solvent. Moreover, catalytic acid Asp181 shifts far away from Cys215 due to the rotation, which might also exert a negative influence on the catalytic activity of PTP1B. The distance between Asp181 and Cys215 is shown in Fig. S3. Generally speaking, distance between Aps181 and Cys215 increases by at least 0.3 nm as a consequence of the rotation of the WPD loop in the three systems.

Fig. 9 Root mean square deviation profile of the WPD loop in the three simulated systems and the superposition of the target conformation with the initial conformation of PTP1B highlighting the different positions of key residues on the WPD loop. (a) and (b) for the open-PTP1B/inhibitor1 complex; (c) and (d) for the open-PTP1B/inhibitor2 complex; (e) and (f) for the closed-PTP1B/inhibitor1 complex. Table 2 The binding free energy andenergy components calculated by MM/GBSA for the three PTP1B/inhibitor complexes.

ACS Paragon Plus Environment

Page 16 of 35

Page 17 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Energy components ∆Eele ∆Evdw ∆EMM ∆Gpolar ∆Gnonpolar ∆Gsol ∆Gele ∆Gbinding

Open-PTP1B/com1 Open-PTP1B/com2 −119.15 −59.85 −178.99 −6.47 140.83 134.36 21.68 −44.63

−8.47 −56.73 −65.52 −5.63 30.41 24.78 21.95 −40.74

Closed-PTP1B/com1 −90.23 −60.17 −150.40 −6.49 114.71 108.22 24.48 −42.18

The conformational transitions of the WPD loop can be reflected from different regions of the free energy landscape (FEL), which is depicted on the basis of the projection of the first two principal components of the C-α trajectory. As shown in Fig. 10, the lowest energy region of the FEL is achieved after 40 ns for the open-mode complexes while the closed-PTP1B/inhibitor1 complex does not become stable until 90 ns, indicating that the protein is in a metastable state during the initial stage. Visual inspection of the structures at low energy reveals that the WPD loop rotates to different degrees compared with the initial structure. The rotation is finally realized by the medium of the open WPD conformation, as depicted in Fig. 10a (13.5 ns), Fig. 10b (10.5 ns) and Fig. 10c (9.5 ns). In addition, the low-energy zone of open-mode PTP1B is larger than that of the closed-mode PTP1B, which suggests that the closed-mode protein undergoes a relatively long transition stage to reach equilibrium. The FEL analysis verifies the spontaneous rotation of the WPD loop under the influence of fullerene inhibitors.

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Fig. 10 Free energy landscape analysis and representative complex structures extracted from the MD trajectory. (a) is open-PTP1B/inhibitor1 complex, (b) is closed-PTP1B/inhibitor1 complex, and (c) is open-PTP1B/inhibitor2 complex. The corresponding structures are represented as cartoon (light pink) and the WPD loop is highlighted in dark blue. The transition of the WPD loop seems to be responsible for the high potency of fullerene inhibitors since it has been reported previously that aryl diketoacid derivatives can function as potent PTP1B noncompetitive inhibitors by targeting the inactive, WPD loop open conformation.23 The substrate cannot be catalyzed even though it binds the active site, because the inhibitors LZP25 and LZP6 block the WPD loop closure. In our work, fullerene derivatives exhibit a similar and might even enhanced inhibitory mechanism compared to the aryl diketoacid derivatives. On the one hand, the large fullerene derivatives can fully occupy the active pocket to prevent the binding of the substrate, which is the factor that needs to be improved for the aryl diketoacid derivatives. Therefore, fullerene derivatives are more likely to be

ACS Paragon Plus Environment

Page 18 of 35

Page 19 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

competitive inhibitors than noncompetitive inhibitors, consistent with the experiment. On the other hand, the fullerene backbone has a greater influence on the structure of side loops and triggers rotation instead of maintenance of the open-mode WPD loop. The rotated WPD loop effectively restrains the formation of the catalytically competent form of the enzyme. The inhibitory potency is closely associated with the high binding affinity as well as the transition of the WPD loop caused by fullerene derivatives.

4. CONCLUSION In the present work, the structural basis of the inhibitory mechanism of the different novel fullerene inhibitors on PTP1B was investigated by MD simulations and MM/GBSA free energy calculations. The docking analysis indicated that the open WPD conformation was preferable for the binding of fullerene inhibitors. The MD simulation of closed-mode PTP1B with inhibitor 1 showed that the inhibitor was finally induced into the active pocket by the strong electrostatic interactions from some positively charged residues of the P-loop and the recognition loop. Both inhibitors showed high binding affinity at the active site even though their dominant driving force might be different. Inhibitor 1 was mainly immobilized by the polar interactions from the P-loop and the recognition loop, resulting in impaired compactness of the protein and scattered correlated motion among important regions. In comparison, inhibitor 2 was mainly stabilized by hydrophobic interaction in addition to several polar interactions generated from the P-loop. And the recognition loop provided little assistance in the binding of inhibitor 2. The high binding affinities of fullerene-based inhibitors resulted from the strong hydrogen bonding of decorated hydrophilic groups with the active site, and van der Waals interactions between nonpolar residues around the active center and the fullerene core. In addition, the WPD loops of the three systems made similar rotations during the simulation. This conformational transition induced by the inhibitors was favored the reduction of the catalytic activity of PTP1B since the rotation prevented closure of the WPD loop. The fullerene inhibitors also blocked the binding of substrate because of its bulky size and thus were thought of as the competitive inhibitors. Above all, it comes to conclusion that both the high binding affinity and the transition of the WPD loop should be responsible for outstanding inhibitory activity of the fullerene inhibitors. The thorough investigation of the inhibitory modes and interactions of fullerene inhibitors can provide guidance to the design of promising, potent nanoparticle inhibitors. ASSOCIATED CONTENT Supporting Information Fig. S1-S3 and Table S1-S4 are shown in the supporting information. AUTHOR INFORMATION Corresponding Author *W. H.: E-mail: [email protected] *S. W.: E-mail: [email protected], Fax: 086-431-88498026 *H. Z.: E-mail: [email protected]

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACKNOWLEDGEMENTS This work was financial supported by the National Natural Science Foundation of China (Grant No. 21443008), Science and Technology Research of Department of Education of Jilin Province (Grant No. 2015468, 2016406 and 2016407), and the Major Scientific Research Projects of Jilin Province (Grant No. 20140203025NY). REFERENCES (1) Koren, S.; Fantus, I. G. Inhibition of the Protein Tyrosine Phosphatase PTPIB: Potential Therapy For Obesity, Insulin Resistance and Type-2 Diabetes Mellitus. Best. Pract. Res Clin Endocrinol Metab. 2007, 21, 621-640. (2) Elchebly, M.; Payette, P.; Michaliszyn, E.; Cromlish, W.; Collins, S.; Loy, A. L.; Normandin, D.; Cheng, A.; Himms-Hagen, J.; Chan, C. C.; Ramachandran, C.; Gresser, M. J.; Tremblay, M. L.; Kennedy, B. P. Increased Insulin Sensitivity and Obesity Resistance in Mice Lacking the Protein Tyrosine Phosphatase-1B Gene. Science 1999, 283, 1544-1548. (3) Bence, K. K.; Delibegovic, M.; Xue, B.; Gorgun, C. Z.; Hotamisligil, G. S.; Neel, B. G.; Kahn, B. B. Neuronal PTP1B Regulates Body Weight, Adiposity and Leptin Action. Nat. Med. 2006, 12, 917-24. (4) Zhang, S.; Zhang, Z. Y. PTP1B as a drug target: Recent Developments in PTP1B Inhibitor Discovery. Drug Discov. Today. 2007, 12, 373-81. (5) Mahadev, K.; Ruffolo, S.; Wu, X. D.; Ouedraogo, R.; Kennedy, B. P.; Goldstein, B. J. Protein-tyrosine Phosphatase 1B (PTP1B), A Target of Cellular Oxidative Inhibition Following Insulin Stimulation, Interacts with the NADPH Oxidase Homolog Nox4. Diabetes. 2006, 55, A301-A301. (6) Zabolotny, J. M.; Bence-Hanulec, K. K.; Stricker-Krongrad, A.; Haj, F.; Wang, Y.; Minokoshi, Y.; Kim, Y. B.; Elmquist, J. K.; Tartaglia, L. A.; Kahn, B. B.; Neel, B. G. PTP1B Regulates Leptin Signal Transduction in Vivo. Dev. Cell. 2002, 2, 489-95. (7) Pei, Z. H.; Liu, G.; Lubben, T. H.; Szczepankiewicz, B. G. Inhibition of Protein Tyrosine Phosphatase 1B as a Potential Treatment of Diabetes and Obesity. Curr. Pharm. Des. 2004, 10, 3481-3504. (8) Zhang, Z. Y. Protein Tyrosine Phosphatases: Structure and Function, Substrate Specificity, and Inhibitor Development. Annu. Rev. Pharmacool. Toxicol. 2002, 42, 209-234. (9) Liu, G.; Xin, Z. L.; Liang, H.; Abad-Zapatero, C.; Hajduk, P. J.; Janowick, D. A.; Szczepankiewicz, B. G.; Pei, Z. H.; Hutchins, C. W.; Ballaron, S. J.; Stashko, M. A.; Lubben, T. H.; Berg, C. E.; Rondinone, C. M.; Trevillyan, J. M.; Jirousek, M. R. Selective Protein Tyrosine Phosphatase 1B Inhibitors: Targeting the Second Phosphotyrosine Binding Site with Non-carboxylic Acid-containing Ligands. J. Med. Chem. 2003, 46, 3437-3440. (10) Liu, G. Recent Advances in Protein-Tyrosine-Phosphatase 1B (PTP1B) Inhibitors for the Treatment of Type 2 Diabetes and Obesity. Drugs of the Future. 2004, 29, 1245-1259. (11) Baskaran, S. K.; Goswami, N.; Selvaraj, S.; Muthusamy, V. S.; Lakshmi, B. S. Molecular Dynamics Approach to Probe the Allosteric Inhibition of PTP1B by Chlorogenic and Cichoric Acid. J. Chem. Inf. Model. 2012, 52, 2004-2012.

ACS Paragon Plus Environment

Page 20 of 35

Page 21 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

(12) Zhang, Z. Y.; Palfey, B. A.; Wu, L.; Zhao, Y. Catalytic Function of the Conserved Hydroxyl Group in the Protein Tyrosine Phosphatase Signature Motif. Biochemistry 1995, 34, 16389-96. (13) Montalibet, J.; Skorey, K.; McKay, D.; Scapin, G.; Asante-Appiah, E.; Kennedy, B. P. Residues Distant from the Active Site Influence Protein-tyrosine Phosphatase 1B Inhibitor Binding. J. Biol. Chem. 2006, 281, 5258-66. (14) Andersen, J. N.; Mortensen, O. H.; Peters, G. H.; Drake, P. G.; Iversen, L. F.; Olsen, O. H.; Jansen, P. G.; Andersen, H. S.; Tonks, N. K.; Moller, N. P. Structural and Evolutionary Relationships among Protein Tyrosine Phosphatase Domains. Mol. Cell. Biol. 2001, 21, 7117-36. (15) Wang, J. F.; Gong, K.; Wei, D. Q.; Li, Y. X.; Chou, K. C. Molecular Dynamics Studies on the Interactions of PTP1B with Inhibitors: From the First Phosphate-binding Site to the Second One. Protein Eng. Des. Sel. 2009, 22, 349-55. (16) Hoff, R. H.; Wu, L.; Zhou, B.; Zhang, Z. Y.; Hengge, A. C. Does Positive Charge at the Active Sites of Phosphatases Cause a Change in Mechanism? The Dffect of the Conserved arginine on the Transition State for Phosphoryl Transfer in the Protein-tyrosine Phosphatase from Yersinia. J. Am. Chem. Soc. 1999, 121, 9514-9521. (17) Zhang, Z. Y.; Wang, Y.; Wu, L.; Fauman, E. B.; Stuckey, J. A.; Schubert, H. L.; Saper, M. A.; Dixon, J. E. The Cys(X)5Arg Catalytic Motif in Phosphoester Hydrolysis. Biochemistry 1994, 33, 15266-70. (18) Fang, L.; Zhang, H.; Cui, W.; Ji, M. J. Studies of the Mechanism of Selectivity of Protein Tyrosine Phosphatase 113 (PTP1B) Bidentate Inhibitors Using Molecular Dynamics Simulations and Free Energy Calculations. J. Chem. Inf. Model. 2008, 48, 2030-2041. (19) Barford, D.; Das, A. K.; Egloff, M. P. The Structure and Mechanism of Protein Phosphatases: Insights into Catalysis and Regulation. Annu. Rev. Biophys. Biomol. Struct. 1998, 27, 133-64. (20) Sanchini, S.; Perruccio, F.; Piizzi, G. Rational Design, Synthesis and Biological Evaluation of Modular Fluorogenic Substrates with High Affinity and Selectivity for PTP1B. ChemBioChem. 2014, 15, 961-76. (21) Khajehpour, M.; Wu, L.; Liu, S.; Zhadin, N.; Zhang, Z. Y.; Callender, R. Loop Dynamics and Ligand Binding Kinetics in the Reaction Catalyzed by the Yersinia Protein Tyrosine Phosphatase. Biochemistry 2007, 46, 4370-8. (22) Xiao, P.; Wang, X.; Wang, H. M.; Fu, X. L.; Cui, F. A.; Yu, X.; Wen, S. S.; Bi, W. X.; Sun, J. P. The second-sphere residue T263 is important for the function and catalytic activity of PTP1B via interaction with the WPD-loop. Int. J. Biochem. Cell Biol. 2014, 57, 84-95. (23) Liu, S. J.; Zeng, L. F.; Wu, L.; Yu, X.; Xue, T.; Gunawan, A. M.; Long, Y. Q.; Zhang, Z. Y. Targeting Inactive Enzyme Conformation: Aryl Diketoacid Derivatives as a New Class of PTP1B Inhibitors. J. Am. Chem. Soc. 2008, 130, 17075-17084. (24) Kumar, R.; Shinde, R. N.; Ajay, D.; Sobhia, M. E. Probing Interaction Requirements in PTP1B Inhibitors: A Comparative Molecular Dynamics Study. J. Chem. Inf. Model. 2010, 50, 1147-1158.

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(25) Zhou, X.; Xi, W.; Luo, Y.; Cao, S.; Wei, G. Interactions of a Water-soluble Fullerene Derivative with Amyloid-beta Protofibrils: Dynamics, Binding Mechanism, and the Resulting Salt-bridge Disruption. J. Phys. Chem. B. 2014, 118, 6733-41. (26) Leonis, G.; Avramopoulos, A.; Papavasileiou, K. D.; Reis, H.; Steinbrecher, T.; Papadopoulos, M. G. A Comprehensive Computational Study of the Interaction between Human Serum Albumin and Fullerenes. J. Phys. Chem. B. 2015, 119, 14971-85. (27) Calvaresi, M.; Bottoni, A.; Zerbetto, F. Thermodynamics of Binding Between Proteins and Carbon Nanoparticles: The Case of C60@Lysozyme. J. Phys. Chem. C. 2015, 119, 28077-28082. (28) Kobzar, O. L.; Trush, V. V.; Tanchuk, V. Y.; Zhilenkov, A. V.; Troshin, P. A.; Vovk, A. I. Fullerene Derivatives as a New Class of Inhibitors of Protein Tyrosine Phosphatases. Bioorg. Med. Chem. Lett. 2014, 24, 3175-9. (29) Trott, O.; Olson, A. J. AutoDock Vina: Improving the Speed and Accuracy of Docking with a New Scoring Function, Efficient Optimization, and Multithreading. J. Comput. Chem. 2010, 31, 455-61. (30) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-balanced, and Scalable Molecular Simulation. J. Chem. Theory. Comput. 2008, 4, 435-447. (31) Canzar, S.; El-Kebir, M.; Pool, R.; Elbassioni, K.; Malde, A. K.; Mark, A. E.; Geerke, D. P.; Stougie, L.; Klau, G. W. Charge Group Partitioning in Biomolecular Simulation. J. Comput. Biol. 2013, 20, 188-198. (32) Malde, A. K.; Zuo, L.; Breeze, M.; Stroet, M.; Poger, D.; Nair, P. C.; Oostenbrink, C.; Mark, A. E. An Automated Force Field Topology Builder (ATB) and Repository: Version 1.0. J. Chem. Theory. Comput. 2011, 7, 4026-4037. (33) Hess, B.; van der Vegt, N. F. A. Hydration Thermodynamic Properties of Amino Acid Analogues: A Systematic Comparison of Biomolecular Force Fields and Water Models. J. Phys. Chem. B. 2006, 110, 17616-17626. (34) Berendsen, H. J.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684-3690. (35) Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: An N⋅ log (N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089-10092. (36) Rastelli, G.; Del Rio, A.; Degliesposti, G.; Sgobba, M. Fast and Accurate Predictions of Binding Free Energies Using MM-PBSA and MM-GBSA. J. Comput. Chem. 2010, 31, 797-810. (37) Case, D. A.; Cheatham, T. E.; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J. The Amber Biomolecular Simulation Programs. J. Comput. Chem. 2005, 26, 1668-1688. (38) Gohlke, H.; Kiel, C.; Case, D. A. Insights into Protein-protein Binding by Binding Free Energy Calculation and Free Energy Decomposition for the Ras-Raf and Ras-RalGDS Complexes. J. Mol. Biol. 2003, 330, 891-913. (39) Vorontsov, I. I.; Miyashita, O. Crystal Molecular Dynamics Simulations to Speed Up MM/PB(GB)SA Evaluation of Binding Free Energies of Di-Mannose

ACS Paragon Plus Environment

Page 22 of 35

Page 23 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Deoxy Analogs with P51G-m4-Cyanovirin-N. J. Comput. Chem. 2011, 32, 1043-1053. (40) Hou, T.; Chen, K.; McLaughlin, W. A.; Lu, B.; Wang, W. Computational Analysis and Prediction of the Binding Motif and Protein Interacting Partners of the Abl SH3 Domain. PLoS. Comput. Biol. 2006, 2, e1. (41) Niu, X.; Gao, X.; Wang, H.; Wang, X.; Wang, S. Insight into the Dynamic Interaction between Different Flavonoids and Bovine Serum Albumin Using Molecular Dynamics Simulations and Free Energy Calculations. J. Mol. Model. 2013, 19, 1039-47. (42) Niu, X.; Wang, X.; Wang, H.; Gao, X.; Wang, Y.; Wang, S. Insights into the Specific Binding Site of Adenosine to the Stx2, the Protein Toxin from Escherichia Coli O157:H7 using molecular dynamics simulations and free energy calculations. Mol. Simul. 2013, 39, 199-205. (43) Batista, P. R.; Costa, M. G.; Pascutti, P. G.; Bisch, P. M.; de Souza, W. High Temperatures Enhance Cooperative Motions between CBM and Catalytic Domains of a Thermostable Cellulase: Mechanism Insights from Essential Dynamics. Phys. Chem. Chem. Phys. 2011, 13, 13709-20. (44) Sarmiento, M.; Puius, Y. A.; Vetter, S. W.; Keng, Y. F.; Wu, L.; Zhao, Y.; Lawrence, D. S.; Almo, S. C.; Zhang, Z. Y. Structural Basis of Plasticity in Protein Tyrosine Phosphatase 1B Substrate Recognition. Biochemistry 2000, 39, 8171-8179. (45) Laatikainen, R.; Saarela, J.; Tuppurainen, K.; Hassinen, T. Internal Motions of Native Lysozyme are More Organized than Those of Mutants: A Principal Component Analysis of Molecular Dynamics Data. Biophys. Chem. 1998, 73, 1-5. (46) Spiliotopoulos, D.; Spitaleri, A.; Musco, G. Exploring PHD Fingers and H3K4me0 Interactions with Molecular Dynamics Simulations and Binding Free Energy Calculations: AIRE-PHD1, a Comparative Study. PloS One. 2012, 7, e46902. (47) Zoete, V.; Michielin, O.; Karplus, M. Relation between Sequence and Structure of HIV-1 Protease Inhibitor Complexes: A Model System for the Analysis of Protein Flexibility. J. Mol. Biol. 2002, 315, 21-52. (48) Zhou, Y.; Zhang, N.; Chen, W.; Zhao, L.; Zhong, R. Underlying Mechanisms of Cyclic Peptide Inhibitors Interrupting the Interaction of CK2alpha/CK2beta: Comparative Molecular Dynamics Simulation Studies. Phys. Chem. Chem. Phys. 2016, 18, 9202-10. (49) Yamada, H.; Mori, S.; Miyakawa, T.; Morikawa, R.; Katagiri, F.; Hozumi, K.; Kikkawa, Y.; Nomizu, M.; Takasu, M. Structural Study of Cell Attachment Peptide Derived from Laminin by Molecular Dynamics Simulation. PloS One. 2016, 11, e0149474. (50) Durdagi, S.; Supuran, C. T.; Strom, T. A.; Doostdar, N.; Kumar, M. K.; Barron, A. R.; Mavromoustakos, T.; Papadopoulos, M. G. In Silico Drug Screening Approach for the Design of Magic Bullets: A Successful Example with Anti-HIV Fullerene Derivatized Amino Acids. J. Chem. Inf. Model. 2009, 49, 1139-43. (51) Innocenti, A.; Durdagi, S.; Doostdar, N.; Strom, T. A.; Barron, A. R.; Supuran, C. T. Nanoscale Enzyme Inhibitors: Fullerenes Inhibit Carbonic Anhydrase by Occluding the Active Site Entrance. Bioorg. Med. Chem. 2010, 18, 2822-8.

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

TOC Graphic

ACS Paragon Plus Environment

Page 24 of 35

Page 25 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Fig. 1 The overall structure of PTP1B with the important regions highlighted in different colors (a) and the key residues on these regions (b). The superimposition of PTP1B with the open and closed WPD loops is shown in (c) and a comparison of the residues in two different modes is depicted in (d). 61x47mm (600 x 600 DPI)

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Fig. 2 Chemical structures of the two fullerene derivatives. 28x10mm (600 x 600 DPI)

ACS Paragon Plus Environment

Page 26 of 35

Page 27 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Fig. 3 The optimal docking models of the three complexes with key residues shown as magenta lines and ligands shown as cyan sticks. (a) Compound 1 docked to open-PTP1B, (b) compound 1 docked to closedPTP1B, (c) compound 2 docked to open-PTP1B. 102x132mm (300 x 300 DPI)

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Fig. 4 (a) The RMSD profiles of the backbone Cα atoms of the three complexes as a function of time. (b) Relative occurrence frequency of the RMSD of the average structure. (c) Average RMSD value of the three systems during the 100 ns MD simulation. 38x8mm (600 x 600 DPI)

ACS Paragon Plus Environment

Page 28 of 35

Page 29 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Fig. 5 (a) and (b) RMSF curves of the Cα atoms of the three complexes. (b), (c) and (d) Visualization of the backbone flexibility for open-PTP1B/inhibitor1, closed-PTP1B/inhibitor1 and open-PTP1B/inhibitor2, respectively. The degree of motion corresponds to the tube with different colors and thickness. 170x95mm (300 x 300 DPI)

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Fig. 6 The cross-correlation matrix map for the open-PTP1B/inhibitor1 complex (a), the openPTP1B/inhibitor2 complex (b) and the closed-PTP1B/inhibitor1 complex (c). 59x19mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 30 of 35

Page 31 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Fig. 7 The motion of PTP1B on the basis of the first two PCs; the arrow represents the direction and amplitude of the movement. (a) open-PTP1B/inhibitor1 complex, (b) closed-PTP1B/inhibitor1 complex, and (c) open-PTP1B/inhibitor2 complex. 103x134mm (300 x 300 DPI)

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Fig. 8 The binding modes of the two inhibitors to the different modes PTP1B and the important residues that contribute to the stability of inhibitors. The first column shows key residues that make hydrogen bonding with the inhibitor, the second column is hydrophobic interaction, and the third column is the decomposition energy of each important residue. The first row is the open-PTP1B/inhibitor1 complex, the second row is closed-PTP1B/inhibitor1 complex and the third row is open-PTP1B/inhibitor2 complex. 139x107mm (600 x 600 DPI)

ACS Paragon Plus Environment

Page 32 of 35

Page 33 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Fig. 9 Root mean square deviation profile of the WPD loop in the three simulated systems and the superposition of the target conformation with the initial conformation of PTP1B, highlighting the different positions of key residues on the WPD loop. (a) and (b) open-PTP1B/inhibitor1 complex; (c) and (d) openPTP1B/inhibitor2 complex; (e) and (f) closed-PTP1B/inhibitor1 complex. 87x96mm (300 x 300 DPI)

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Fig. 10 Free energy landscape analysis and representative complex structures extracted from the MD trajectory. (a) open-PTP1B/inhibitor1 complex, (b) closed-PTP1B/inhibitor1 complex, and (c) openPTP1B/inhibitor2 complex. The corresponding structures are represented as cartoons (light pink), and the WPD loop is highlighted in dark blue.

154x299mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 34 of 35

Page 35 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

34x24mm (300 x 300 DPI)

ACS Paragon Plus Environment