Design, synthesis, biological evaluation and molecular dynamics simulation studies of (R)-5-methylthiazolidin-4-One derivatives as megakaryocyte protein tyrosine phosphatase 2 (PTP-MEG2) inhibitors for the treatment of type 2 diabetes

Article abstract of DOI:10.1080/07391102.2019.1654410

PTP-MEG2 plays a significant role in insulin production and is able to enhance insulin signaling and improve insulin sensitivity. So, PTP-MEG2 inhibitors are closely associated with type 2 diabetes therapy. A series of novel (R)-5-methylthiazolidin-4-one

Full text of DOI:10.1080/07391102.2019.1654410

Journal of Biomolecular Structure and Dynamics
ISSN: 0739-1102 (Print) 1538-0254 (Online) Journal homepage: https://www.tandfonline.com/loi/tbsd20
Design, synthesis, biological evaluation and
molecular dynamics simulation studies of
(R)-5-methylthiazolidin-4-one derivatives as
megakaryocyte protein tyrosine phosphatase 2
(PTP-MEG2) inhibitors for the treatment of type 2
diabetes
Jingwei Wu, Yingzhan Sun, Hui Zhou, Ying Ma & Runling Wang
To cite this article: Jingwei Wu, Yingzhan Sun, Hui Zhou, Ying Ma & Runling Wang (2019):
Design, synthesis, biological evaluation and molecular dynamics simulation studies of (R)-5-
methylthiazolidin-4-one derivatives as megakaryocyte protein tyrosine phosphatase 2 (PTP-MEG2)
inhibitors for the treatment of type 2 diabetes, Journal of Biomolecular Structure and Dynamics,
DOI: 10.1080/07391102.2019.1654410
To link to this article: https://doi.org/10.1080/07391102.2019.1654410
Accepted author version posted online: 11
Aug 2019.
Submit your article to this journal
View Crossmark data
Full Terms & Conditions of access and use can be found at
https://www.tandfonline.com/action/journalInformation?journalCode=tbsd20

Design, synthesis, biological evaluation and
molecular dynamics simulation studies of (R)-5-
methylthiazolidin-4-one derivatives as
megakaryocyte protein tyrosine phosphatase 2 (PTP-
MEG2) inhibitors for the treatment of type 2
diabetes
Jingwei Wu a, ‡, Yingzhan Sun a, ‡, Hui Zhou a, Ying Ma a, *, Runling Wang a, *
aTianjin Key Laboratory on Technologies Enabling Development of Clinical Therapeutics and
Diagnostics (Theranostics), School of Pharmacy, Tianjin Medical University, Tianjin 300070,
China
‡ These authors contributed equally to this work
*Corresponding authors: Ying Ma, E-mail: maying@tmu.edu.cn
1

Runling Wang , E-mail:wangrunling@tmu.edu.cn
Keywords: PTP-MEG2; type 2 diabetes; De novo design; synthesis; molecular dynamics
simulation
Abstract
PTP-MEG2 plays a significant role in insulin production and is able to enhance insulin signaling
and improve insulin sensitivity. So, PTP-MEG2 inhibitors are closely associated with type 2
diabetes therapy. A series of novel (R)-5-methylthiazolidin-4-one derivatives were designed and
synthesized, and their PTP-MEG2 inhibitory activities (IC50) were determined. Among the
desired compounds, 1h shares the highest inhibitory activity (IC50 = 1.34 μM) against PTP-
MEG2. Additionally, various post-dynamic analyses confirmed that when compound 1h bound to
the PTP-MEG2, the protein conformations became unstable and the function of the pTyr
recognition loop (Asn331-Cys338) would be disturbed. And thus, the ideal conformations needed
for the catalytic activity was difficult to be maintained. In brief, these might be how the
compound 1h worked. Furthermore, we also found that the key residues Arg332 would play a
critical role in disturbing the residue interactions.
List of abbreviations
DCCM (dynamic cross-correlation mapping); DMF (N,N-dimethylformamide); DSSP (definition
of secondary structure of proteins); FOXO (forkhead transcription factors); MD (molecular
dynamics); PCA (principal component analysis); PDB (protein data bank); PTKs (protein
tyrosine kinases); PTPs (protein tyrosine phosphatases); PTP-MEG2 (megakaryocyte protein
2

tyrosine phosphatase 2); RIN (residue interaction network); RING (Residue Interaction Network
Generator); RMSD (root means square deviation); RMSF (root mean square fluctuation);
1 Introduction
Protein tyrosine phosphatases (PTPs) comprise a superfamily of signaling enzymes which are
important in lots of cellular processes, such as cell proliferation, differentiation, migration,
adhesion and apoptosis(Hendriks et al. 2013; Lee, Yi, Lawan, Min & Bennett 2015). PTPs are
catalyzed by protein tyrosine kinases (PTKs), which are vital for proper level of cellular protein
tyrosyl phosphorylation. Alterations in PTP genes contribute to human diseases, including cancer
or metabolic disorders(He, Zeng, He, Zhang & Zhang 2013; Hendriks, Bourgonje, Leenders &
Pulido 2018). Therefore, PTPs are promising drug targets(Zhang 2001).
The non-receptor PTP-MEG2 (also namedPTPN9)(Wang et al. 2019) is the unique mammalian
representative of PTPs, which is highly expressed at highest levels in leukocytes, brain, exocrine
and endocrine cells(Qi et al. 2002; Huynh et al. 2003). It has been proved that PTP-MEG2 is a
regulator of insulin-dependent FOXO1 subcellular localization and controls insulin signaling and
insulin production (Cho et al. 2006; Gurzov, Stanley, Brodnicki & Thomas 2015; Chatterjee,
Khunti & Davies 2017). Indeed, treatment with PTP-MEG2 inhibitors can enhance insulin action
both in vivo and in vitro, indicating that PTP-MEG2 may be an effective therapeutic target of
type 2 diabetes(Cho et al. 2006).
The PTP catalytic site (i.e., pTyr) alone is not sufficient for high-affinity binding while the
adjacent residues are also conductive to PTP substrate recognition. Reports have showed that
there are neighbouring sub-pockets, which can be targeted for PTP-MEG2 inhibitors(Zhang et al.
2009). Inspired by these discoveries, we design a series of PTP-MEG2 inhibitors by means of
3

fragment-based de novo design involving virtual screening, molecular docking and molecular
dynamics simulation., In order to discover some potential PTP-MEG2 inhibitors for the treatment
of type 2 diabetes, a series of thiazolidin-4-one derivatives were synthesized and their PTP-
MEG2 inhibitory activities (IC50) were determined. In the meantime, molecular dynamics
simulations could provide the fundamental information concerning individual particle motions
over time. Therefore, the application of the 100 ns MD molecular dynamics could explore the
effect of 1h on the PTP-MEG2, and various post-analyses were carried out to analyze the
mechanisms of inhibition at the molecular level. The root means square deviation (RMSD) firstly
verified the PTP-MEG2 and PTP-MEG2/1h systems were stable and available. The root mean
square fluctuation (RMSF) and radius of gyration (Rg) analyses showed that the protein
conformation became unstable and the fluctuation of the the pTyr recognition loop (Asn331-
Cys338) was decreased in the PTP-MEG2/1h system. The Dynamics cross-correlation matrix
analysis (DCCM) results and principal component analysis (PCA) had similar results with the
above Rg and RMSF analyses, which further verified that when inhibitor 1h bound to the protein,
the PTP-MEG2 owned unstable conformation and the pTyr recognition loop had low flexibility.
The study on the residue interaction network (RIN) revealed that the key residues Arg332, played
an important role in breaking the residue interactions in the PTP-MEG2/1h system. In brief, the
inhibitor owned large inhibitor effect for the PTP-MEG2 in this study.
2 Material and Methods
Our calculations were carried out on Dell Precision TM T5500 computer with Discovery Studio
3.5 software package (http://accelrys.com/).
2.1 Compound design
4

2.1.1 Virtual screening
Virtual screening has become an effective method for the lead discovery and optimization.
(Shoichet 2004; Li et al. 2009; Liu et al. 2011). In the present study, docking-based virtual
screening was performed using LibDock (Accelrys Discovery Studio package), because it has
been applied to the GSK validation data set(Diller & Merz 2001).
The receptor crystal structure of catalytic domain of PTP-MEG2 (PDB ID:4GE6) (Zhang et al.
2009) was downloaded from the protein data bank (PDB) (Berman et al. 2002)and was used for
molecular docking. The ligands utilized for virtual screening were retrieved from Drugs-Now
database and Drugs-Like in ZINC database(Irwin & Shoichet 2005). The preparation for all
compounds and the protein were performed by the Prepare Ligands and Prepare Protein
protocols, respectively. The binding sites sphere was defined using “Find Site From Current
Selection” based on the key residues in the catalytic pocket (The pTyr recognition loop, residues
331–338) and the sub-pocket (residues 514–521, residues 558–564)(Zhang et al. 2009). All
investigated ligands from ZINC database were docked to the binding pocket with the method of
LibDock. According to the docking scores, the best candidate was identified as potential scaffold
used for the following study of De Novo Design.
2.1.2 De Novo Design
The technology of “De novo design” has been often used in the design of drug candidates for lots
of pharmaceutically relevant targets (Loving, Alberts & Sherman 2010). The flowchart showing
the process of De novo design was illustrated in Figure 1.The protocol of “De Novo Library
Generation” was used to create fragment library, which was used in the following step “De Novo
Evolution”. In view of pharmacological activity, the fragments were retrieved from drugs-now
5

database. The protocol of De Novo Evolution could develop the new derivatives in the binding
pocket environment of the receptor based on the scaffold native structure. The Ludi
algorithm(Böhm 1992) was used to fit appropriate fragments to scaffold obtained by virtual
screening, and then produced a collection of the new derivatives with higher scores(Law et al.
2009). The evolution mode was set to full evolution, which allows the scaffold to link up to a
maximum of three fragments. The maximum number of survivors, population size and
generations were set to 2, 10 and 2, respectively(Wu et al. 2014). Finally, top-ranked “LibDock
SCORE” molecules were selected for the following chemical synthesis.
Figure 1. The de novo design flowchart. The database and the methods used for de novo design
were depicted as blue cylinder and yellow arrows, respectively.
2.2 Chemistry
Melting points were measured with an RY-2 microscopic melting point apparatus and are
uncorrected. 1H NMR and 13C NMR spectra were recorded on a Bruker AV400 NMR
spectrometer with CDCl3 or DMSO-d6 as solvent and known chemical shifts of residual proton
6

signals of deuterated solvents ( 1H NMR: 7.26 ppm for CDCl3) or carbon signals of deuterated
solvents (13C NMR: 39.52 ppm for DMSO-d6) as internal standard. MS spectra were obtained
with an Agilent Q-TOF 6510 mass spectrometer with electro-spray ionization technique (in
ESI+/ESI- mode) using direct injection method. The synthetic route to the desired products 1a-1h
is illustrated in Schemes 1. For synthetic details, see the supporting information.
Scheme 1 Synthetic routes to desired products 1a-1h
Reagents and conditions: a. 1,2-dibromoethane (for 3a-3d) or 1,4-dibromobutane (for 3e-3h), aq.
K2CO3, EtOH/H2O (v/v = 1:1), 80 °C; b. (R,Z)-ethyl 2-(5-methyl-4-oxothiazolidin-2-
ylidene)acetate, KI, K2CO3, DMF, r.t.
2.3 Enzymatic assays
The activity of PTP-MEG2 was assessed at 37 °C in 96-well plates with p-nitrophenyl phosphate
(pNPP) as substrate. Purifed recombinant PTP-MEG2 (0.1 μg) in 60 μL buffer with 60 mM
calcium citrate, 2 mM EDTA, and 2 mM dithiothreitol and each sample were added to each well
of a 96-well plate. After preincubation for 2 min at 37°C, 50 μL of buffer with 40 μL pNPP was
added and incubated. After incubation for 30 min, the reaction was quenched by adding 10μL of
7

a 0.2 N NaOH solution, and the absorbance at 405 nm was measured. IC50 values were
determined as the concentration of compounds that gave 50% of the control enzyme activity.
2.4 Molecular Dynamics Simulations
2.4.1 System preparation
The X-ray crystal structure of PTP-MEG2 (PDB ID: 4GE6) was derived from the Protein Data
Bank (PDB) and prepared by the “prepare protein” in DS v3.5(Ma, Wang, Xu, Wang & Chou
2012). The process of the “prepare protein” included the following steps: firstly, the H2O in the
protein was removed; secondly, the intramolecular bond orders was assigned; thirdly, the
hydrogen atoms was added to the structure; and then, the disulfides and metals in protein was
treated; finally, Modeler (Chimera) was used to add the missing residues(Cheng et al. 2012;
Zhang, Wang, Xu, Wang & Wang 2012). The structure of the PTP-MEG2 was crystallized as a
dimer; only chain A were considered for the simulations, therefore, the chain B was deleted. And
then, the molecular dynamic simulations were carried out.
2.4.2 Molecular Dynamics Simulations
Conventional MD simulations in the study were performed to assess the effect of complex 1h.
Both PTP-MEG2 and PTP-MEG2/1h systems were performed to run 100 ns MD
simulations(Fakhar et al. 2017; Li et al. 2018). With the help of the single-point charge (SPC)
water model(MacKerell et al. 1998) and periodic boundary conditions, the two proteins were
solvated in dodecahedral 1.0 nm, respectively. The PTP-MEG2 system was solvated in about
3547 SPC water molecules and the PTP-MEG2/1h system was solvated in about 3552 SPC water
molecules. To neutralize the positive charge in the two systems, three chlorine (3 Cl- ions)
8

counter ions were added to replace the solvent molecules randomly. In the process of the
minimization, there were 50,000 steps in the simulations to perform steepest descent energy
minimization until the maximum force was less than 10.0 KJ/mol(Berhanu & Masunov 2014).
Afterwards, the temperature and pressure of the two systems were stable at 310 K and 1 bar
under the 100 ps NVT and NPT simulation conditions, respectively(Foster & West 2017).
Meanwhile, the particle-mesh Ewald algorithm was utilized to manage the long-range
electrostatic interactions, and bond lengths was constrained at their equilibrium values by the
LINCS algorithm(Li et al. 2018). The two systems were performed 100 ns molecular dynamics
simulations, respectively. The trajectories were stored in the traj.trr file for the two models at
every 20 ps, and structural analysis was done at every picosecond. Post dynamics analysis were
carried out to compare the two models to find out the effect of the compound 1h on the PTP-
MEG2, including the RMSD, Rg, RMSF, PCA, DCCM and RIN analyses.
2.4.3 Principal component analysis
PCA has been utilized to give a deep insight of the dynamics of a protein system(Sittel, Jain &
Stock 2014). In this study, PCA was carried out on Cα atoms on 1000 snapshots extracted from
8~100 ns MD trajectories. The conformational changes of proteins and the atomic displacement
in a collective manner can be defined using PCA(Buslaev, Gordeliy, Grudinin & Gushchin
2016). The correlated motions of whole protein can be represented by the eigenvectors and
eigenvalues. Therefore, the PC1, PC2 and PC3 (first three principal components) were computed,
which contributed more significantly in PCA analysis. The PCA method could be calculated by
the following formula(Balmith & Soliman 2017):
〈
〈 〉 〉( )
퐶푖푗 = (푋푖 − 푋 ) 푖, 푗 = 1,2,3 … ,3푁
푗
9

Where, the i/j-th alpha carbon atoms are represented by the Xi/j, the time average is represented
by the 〈Xi/j〉 and N is considered as the number of Cα atoms. The tool of the Bio3D library in the
R and ProDy software are utilized to carry out the analysis.
2.4.4 Dynamics cross-correlation matrix analysis
To gather a deeper insight of the dynamics of PTP-MEG2 and PTP-MEG2/1h structures, the
DCCM was used to determine cross-correlated displacements of backbone alpha carbon atoms in
the molecular dynamic trajectories for the two systems, and was analyzed by the following
formula(Lukman, Lane & Verma 2013):
(∆푟 × ∆푟 )
푖
푗
퐶푖푗 =
⁄
1 2
2
〈
2〉〈
〉
( ∆푟 ∆푟 )
푖
푗
Here, i/j corresponds to i/j-th residues and ∆푟 represents the displacement of i/j-th atoms from
푖/푗
the mean, respectively. This analysis was also performed by the Bio3D library in the R and
ProDy software(Aier, Varadwaj & Raj 2016).
2.4.5 Residue interaction networks analysis
RIN analysis could show the detailed interactions between residues in protein(Niu et al. 2017),
which was an alternative way of representing protein structures and has been utilized for
analyzing the structural changes of protein interaction network with inhibitor, protein folding and
catalytic activity. In the RINs, the edges represented the interactions between pairs of residues
and the nodes represented residues of a protein. According to the type of interaction, the edges
were labeled as van der Waal, salt bridges and hydrogen bonds interactions, respectively. The
Cytoscape network file for the PTP-MEG2 and PTP-MEG2/ inhibitor systems could be
10

constructed by the RING v2.0.1(Piovesan, Minervini & Tosatto 2016). The Cytoscape
3.6.1(Shannon et al. 2003) was utilized to constructed the RINs, using the Cytoscape network
file. The tool of RINalyzer in Cytoscape 3.6.1 was used to observe the network. The closeness
centrality, shortest path betweenness and the degree of residues could be calculated by the tool of
RINalyzer and NetworkAnalyzer in the Cytoscape.
3. Results and discussion
3.1 Compounds design
Table 1. Structure and PTP-MEG2 inhibitory activity of compounds 1a-1h and ZINC5852926
De Novo Evolution
ZINC5852926
1a-1h
LibDock
Compounds
R
n
MW
IC50 (μM)
Score
-
-
201.243
393.45
393.45
366.39
101.732
129.736
135.991
131.773
>50
ZINC5852926
p-COOEt
m-COOEt
p-NO2
1
1
1
8.41±0.51
7.82±0.46
13.33±0.78
1a
1b
1c
11

p-CH2COOEt
p-COOEt
1
3
3
3
3
407.48
406.47
406.47
394.44
435.53
137.817
128.492
130.559
132.410
136.556
3.62±0.21
2.77±0.16
5.03±0.31
7.54±0.45
1.34±0.08
1d
1e
1f
m-COOEt
p-NO2
1g
1h
p-CH2COOEt
12

Figure 2. Structures and binding interactions of 4GE6 in complex with ZINC5852926 and
compound 1h. (A) 2D diagram of the interaction of ZINC5852926- PTP-MEG2. Residues
involved in hydrogen-bond, charge or polar interactions are represented by pink boxes. Residues
involved in van der waals interactions are represented by green boxes. (B) 2D diagram of the
interaction of compound 1h-PTP-MEG2. Hydrogen-bond interactions with amino acid main
13

chains and side chains are represented by a green and blue dashed arrow directed towards the
electron donor. (C) The docked pose of ZINC5852926 in the catalytic pocket of 4GE6. (D) The
docked pose of compound 1h in the catalytic pocket and the adjacent sub-pocket of PTP-MEG2.
(E) The aligned pose of ZINC5852926 (yellow) and compound 1h (green) in the active sites of
4GE6. (F) Surface representation of 4GE6 in complex with ZINC5852926 (yellow) and
compound 1h (green). Hydrogen bond interactions were depicted as yellow-dashed lines. The
Figures were created using Discovery Studio 3.5 software package (http://accelrys.com/) and
PyMol (http:/www.pymol.org).
ZINC5852926, as the best scaffold, was discovered in the virtual screening procedure (Figure 1).
The following de novo design led to the discovery of 8 novel compounds (1a-1h) by the protocol
of “De Novo Evolution”. The properties and doking scores of 8 compounds (1a-1h) as well as
the scaffold (ZINC5852926) were summarized in Table 1. The LibDock score is a measure of
the strength of binding between the ligand and the receptor. The higher score indicates stronger
binding interaction between them(Diller & Merz Jr 2001). Examination of the LibDock score
(Table 1) revealed that the scores for those 8 compounds have exceeded that of the scaffold
(ZINC5852926), signifying that these desired compounds (1a-1h) would bind to the PTP-MEG2
more strongly than the scaffold (ZINC5852926).
The 2D diagrams of ZINC5852926-PTP-MEG2 as well as compound 1h- PTP-MEG2 were
shown in Figure 2A and Figure 2B, respectively. Comparing with the Figure 2A, we could find
more interactions formed between compound 1h- PTP-MEG2 in Figure 2B, green boxes such as
Arg312, Asn314, Pro315, Val316, Phe319, Leu339, Thr552, Ala555, Ser557 and Gln559
represented van der waals interactions, while pink boxes like Ile310, Arg311, Glu313, Asp335,
Arg551 and Phe556 were involved in hydrogen bond interaction, polar or charge interaction.
14

Moreover, the oxygen atoms formed one potential H-bond with the main chain of Phe556 (green
dashed arrow in Figure 2B) and three potential H-bond with the side chains of Arg311 and
Arg551 (blue dashed arrow in Figure 2B). The docked conformations of ZINC5852926 and 1h
to PTP-MEG2 were showed in Figure 2C, Figure 2D and Figure 2E, respectively. Obviously,
compounds 1h (green, Figure 2E) occupied both of the catalytic pocket and the adjacent sub-
pocket while ZINC5852926 (yellow, Figure 2E) could only bind the catalytic pocket. In
addition, the binding surfaces of 4GE6 in complex with ZINC5852926 (yellow) and compound
1h (green) were shown in Figure 2F, compound 1h formed the important H-bond (depicted as
yellow-dashed lines) in the sub-pocket. This may be able to explain why the compound 1h
inhibited the PTP-MEG2 with a lower IC50 than ZINC5852926.
3.2 Chemistry
Firstly, the synthetic route to desired products 1a-1h was shown in Scheme 1. O-alkylation of
phenol derivatives (2a-2d) with 1,2-dibromoethane (for 3a-3d) or 1,4-dibromobutane (for 3e-3h)
in the presence of K2CO3 in EtOH/H2O at 80°C smoothly afforded 3a-3h according to known
procedures(Zhang et al. 2017). Secondly, N-alkylation of (R,Z)-ethyl 2-(5-methyl-4-
oxothiazolidin-2-ylidene)acetate with 3a-3h in the presence of KI as catalyst and K2CO3 as base
in DMF at room temperature produced the final desired products 1a-1h(Wu et al. 2019).
3.3 In vitro PTP-MEG2 inhibitory activity
The results of in vitro inhibitory assay of 8 synthesized compounds (1a-1h) as well as scaffold
ZINC5852926 against PTP-MEG2 are summarized in Table 1. It is very clear that all the desired
compounds (1a-1h) showed significantly improved PTP-MEG2 inhibitory activity ( IC50 values
were range from 1.34 μM to 13.33 μM) as compared with scaffold ZINC5852926 (IC50 >50 μM).
15

Besides, it was found that the PTP-MEG2 inhibitory activity slightly increased as
the length of alkyl groups increased (1a < 1e, 1b < 1f, 1c < 1g, 1d < 1h). Moreover, it should be
noted that the best and the worst R group in the benzene ring was p-CH2COOEt and p-NO2,
respectively (1d > 1b> 1a > 1c, 1h > 1e >1f >1g). In summary, in vitro inhibitory assay of 8
synthesized compounds led to the discovery of a potent novel compound 1h (IC50 = 1.34 μM),
which had more inhibitory activity than the scaffold ZINC5852926 .
3.4 Post-dynamics analyses
3.4.1 Root mean square deviation (RMSD)
It could be seen from the Figure 3 that both the two systems reached equilibrium and
stabilization from 8 ns simulation. This result illustrated that the 100 ns molecular dynamics
simulations for the two systems were credible and available.
Figure 3 Cα backbone RMSD for PTP-MEG2 system (blue line) and PTP-MEG2/1h system
(green line) over a 100 ns MD simulation.
16

As shown in the Figure 3, the results indicated more flexible nature of residues in PTP-MEG2/1h
system than that in PTP-MEG2 system. The results showed that the fluctuation of Cα in the PTP-
MEG2/1h system was increased, which could obviously show that the stability of conformations
of PTP-MEG2/1h would be reduced.
3.4.2 Radius of gyration (Rg) and Root mean square fluctuation (RMSF) analyses.
Figure 4 Rg comparison across the 100 ns MD simulation of the PTP-MEG2 system (blue line)
and PTP-MEG2/1h system (green line).
As shown in the Figure 4, to provide further insight into the stability of biological molecules, as
well as to determine the compactness of the tertiary structure of a protein throughout the MD
simulation, Rg was calculated(Ndagi, Mhlongo & Soliman 2017). It could be seen form the
Figure that the whole trend of the fluctuations of the two systems was similar compactness.
However, we could still discover that the PTP-MEG2 system showed a lower average Rg of 1.86
Å than PTP-MEG2/1h system with an average Rg of 1.88 Å. The increase might reflect highly
unstable nature of PTP-MEG2/1h system as compared to the PTP-MEG2 system. This was
17

agreed with the above RMSF and RMSD assumption that the compound 1h decreased the
interactions between the amino acids and the amino acids resulting in the changes of the
conformational stability in the PTP-MEG2/1h.
The RMSF is a measure of the average atomic mobility of backbone atoms (N, Cαand C) during
MD simulation. As shown in the Figure 5, RMSF was monitored throughout the 100 ns
molecular dynamics simulation. The function of a protein does not only depend on its structure
but also on its dynamics. The function of a protein can be altered through interactions with the
active site, or indirectly by disrupting motions essential to its function; more specifically, the
essential motions induced by the conformational changes. It was obviously found that the
fluctuations of three regions (Leu328-Asp335, Tyr366-Asn370 and Val495-Ser502) were
reduced in the PTP-MEG2/1h system, with the average value 0.229 nm, 0.180 nm and 0.222 nm,
respectively. In the PTP-MEG2 system, the average values of these regions were 0.314 nm, 0.272
nm and 0.337 nm, respectively. The differences of these regions between the PTP-MEG2 system
and the PTP-MEG2/1h system were 0.085 nm, 0.092 nm and 0.115 nm. The results suggested
that the compound 1h had effect on the three regions with decreasing the interaction between the
residues and residues. Interestingly, the conformations of residues of one region (Leu328-
Asp335) partly located in the pTyr recognition loop (Asn331-Cys338)(Barr et al. 2009; Zhang et
al. 2012) might be changed, such as Asn331, Arg332, Tyr333, Gly334 and Asp335. We inferred
that the compound 1h would influence the function of this loop region, which could be one of the
reasons for its inhibitory effect. It also could be seen that the whole fluctuations of PTP-
MEG2/1h system was larger than that in the PTP-MEG2 system, indicating that the average
atomic mobility of backbone atoms (N, Cα and C) during MD simulation was increased in the
PTP-MEG2/1h system. Therefore, the large fluctuations of whole system would change the
18

stability of conformations of PTP-MEG2/1h. The conformations needed for the catalytic activity
could be no longer maintained. The results would be another reason for 1h inhibitory effect.
To sum up, when the compound 1h bound to the PTP-MEG2, the conformations of the protein
became unstable, and thus, we speculated that the conformations needed for the catalytic activity
might be no longer maintained due to the increased fluctuations of structure of the PTP-
MEG2/1h. Meanwhile, the fluctuation of the the pTyr recognition loop (Asn331-Cys338) was
decreased in the PTP-MEG2/1h system, inferring that the conformations of residues in this
region might be changed and the interactions between residues in this region and residues in other
regions might also be changed, therefore, it could be deduced that the catalytic activity might be
directly influenced by the inhibitor 1h. These results might be the way that inhibitor 1h inhibited
the PTP-MEG2.
Figure 5 RMSF of the PTP-MEG2 system (blue line) and PTP-MEG2/1h system (green line)
over 100 ns MD simulation. The red boxes represent the large differences between the two
systems.
19

3.4.3 Principal component analysis (PCA)
To gain insight into the motions associated with conformation, we examined the conformational
behavior of PTP-MEG2 system and PTP-MEG2/1h system by projecting them along the first
three principal component or eigenvector (PC1, PC2 and PC3) directions.
Figure 6. PCA scatter plots projecting principal components (PC1, PC2 and PC3) in the space of
(A) PTP-MEG2 system and (B) PTP-MEG2/1h system.
As shown in the Figure 6, eigenvectors computed from the MD trajectories for the two systems
were quite varied, displaying significant differences in motion across the three principal
components, which clearly indicated the differences in protein motion and the conformational
landscape between the PTP-MEG2 and the PTP-MEG2/1h. In the PTP-MEG2/1h system (Figure
6B), there were the highest correlated motion along PC1, PC2 and PC3 with a percentage of
41%, 8.3% and 4.3%,respectively, whereas in the PTP-MEG2 system (Figure 6A), the values of
the three PCs were 21.9%, 11.6% and 6.7%, respectively, suggesting that the correlated motion
was much more in the PTP-MEG2/1h system than in PTP-MEG2 system. Meanwhile, we could
20

also see from the scatter plot that there were much more the light red or blue dots in the PTP-
MEG2/1h system than in the PTP-MEG2 system, indicating that the unstable conformations were
increased in the PTP-MEG2/1h system due to the compound 1h. The results showed that the
PTP-MEG2/1h system was more flexible than the PTP-MEG2 system, which was agreed with
the above RMSF, RMSD and Rg analyses.
3.4.4 Dynamics cross-correlation matrix (DCCM) analysis
To get a further insight of the confrontational changes of PTP-MEG2 upon the compound 1h,
DCCM analysis was carried out to determine the presence of correlated motions. As shown in the
Figure 7, the DCCM of the two systems was compared. The cyan and light cyan represented the
positive correlated motions of residues, whereas the pink and light pink represented the anti-
correlated movements of residues.
Figure 7. Cross-correlation matrices of the fluctuations of C-α atoms in PTP-MEG2 system (A)
and PTP-MEG2/1h (B).
21

DCCM analysis showed that binding of 1h changed the structure conformation of PTP-MEG2 as
reflected in the correlated motions and dynamics. Firstly, it could be seen from the Figure7 that
the the positive correlated motions and the anti-correlated motions of residues were both
increased in the PTP-MEG2/1h system, which indicated that fluctuations of the residues in the
system was enhanced. Therefore, the stability of conformations of the PTP-MEG2/1h was
decreased, which was agreed with the above analyses. Secondly, it could be also shown that the
correlations of three regions marked by the blue boxes changed a lot when the compound 1h
bound to the protein. The positive correlated motions of residues between the region Leu328-
Cys338 and Asp335-Asp340 were subdued, interestingly, we could infer that the fluctuations of
residues in the pTyr recognition loop, such as Asn331, Arg332, Tyr333, Gly334，Asp335,
Val336, Pro337 and Cys338 would be decreased due to the decreased correlations between this
region and other region. Similarly, both the anti-correlated movements of residues between
His352-Phe362 and Gly377-Asp385 and that between Val495-Pro507 and Ile530-Thr540 were
also diminished, indicating that the correlations of these regions were reduced in the PTP-
MEG2/1h system.
Thus, the fluctuations of these regions were decreased in this system, which was agreement with
the RMSF analysis. In brief, DCCM analysis showed that the protein conformations became
unstable and the fluctuation of pTyr recognition loop (Asn331-Cys338) was inhibited owing to
the decreased correlations between this region and other region，which was agreement with the
above analyses.
3.4.5 Residue interaction network analysis
22

To further explore the effect of the compound 1h on the residues interactions in PTP-MEG2, the
residue interaction network (RIN) analysis was performed. Firstly, the closeness centrality and
shortest path betweenness of the residues were calculated by the tool of RINalyzer and
NetworkAnalyzer in the Cytoscape. And thus, the notable residues in protein could be screened
out, which could be utilized to reflect whether the protein system was stable or not. It could be
seen from the table S1 that there were 10 notable residues, such as Cys567, Met550, Ser548,
Asp529, Leu528, Phe525, Met400, Val399, Asp385 and Asn359 in the PTP-MEG2 system,
whereas in the PTP-MEG2/1h system, the number of the notable residues was 7, and they were
Met550, Phe525, Met400, Val399, Asp385, Asn359 and Arg332, respectively. We could find
that the number of notable residues was decreased, which might cause the protein conformation
unstable. Meanwhile, the Arg332 was found as notable residue, with the closeness centrality
value 1.0 and shortest path betweenness value 0.1327. Secondly, the degree of the residues was
calculated to further study the key residue. The three regions (Leu328-Asp335, Tyr366-Asn370
and Val495-Ser502) which showed large differences in above analysis were carried out to
analyze. The table S2 showed that the Arg332 owned biggest changes, and the difference of it
between PTP-MEG2/1h and PTP-MEG2 was 3. In brief, the Arg332 might be the key residues
which would play an important role in disturbing the residue interaction. Furthermore, the
detailed changes of the correlations in the above analyses were shown in Figure 8. The red doted
lines and the red rectangles represented the interactions and residues in the PTP-MEG2/1h,
respectively, whereas the green doted lines represented the interactions in the PTP-MEG2
system. The black doted lines represented the interactions in both systems. Thus, we could find
the Arg332 owned three interactions with Glu329, Leu379 and Leu328 in the PTP-MEG2/1h
system, whereas in the PTP-MEG2 system it had only one interaction with Glu329. To sum up,
we speculated that the Arg332 would be the key residues.
23

Figure 8 The residue interaction network analysis for the PTP-MEG2 and PTP-MEG2/1h
systems. The red doted lines and the red rectangles represented the interactions and residues in
the PTP-MEG2/1h, respectively, whereas the green doted lines represented the interactions in the
PTP-MEG2 system. The black doted lines represented the interactions in both systems.
4 Conclusions
In conclusion, a series of novel (R)-5-methylthiazolidin-4-one derivatives targeting PTP-MEG2
were designed with the aid of de novo design method. The designed compounds were
synthesized, and their activities against PTP-MEG2 were evaluated. The results showed that the
(R)-5-methylthiazolidin-4-one derivatives inhibited PTP-MEG2 within the micromolar range
(1.34 μM ~ 13.33 μM). Compound 1h, displaying the best inhibitory activity with an IC50 value
of 1.34 μM, was further validated by the method of molecular dynamics simulations. In this
study, the various post-dynamics analyses were used to explore the effect of the compound 1h on
the PTP-MEG2, such as the RMSD, Rg, RMSF, PCA, DCCM and RIN analyses. The results
suggested that when the compound 1h bound to the PTP-MEG2, the protein conformations
became unstable, indicating that the conformations needed for the catalytic activity was no longer
maintained. Meanwhile, the fluctuations of the pTyr recognition loop (Asn331-Cys338) was
increased in the PTP-MEG2/1h system, indicating that the conformations of this region might
24

also be changed, and thus, the function of catalytic activity site might be inhibited by the 1h. This
is exactly how the compound 1h works, disturbing not only the pTyr recognition loop
conformations, but also the protein conformations. Furthermore, we also found the key residues
Arg332 would play an important role in disturbing the residue interactions.
Supporting Information.
The synthetic details of the desired products 1a-1h and Tables (Table S1 and Table S2) for
residue interaction network analysis have been provided.
Acknowledgment
This work was supported by the National Natural Science Foundation of China (91773569),
Natural Science Foundation of Tianjin (16JCZDJC32500 and 18JCQNJC13700) and Science &
Technology Development Fund of Tianjin Education Commission for Higher Education
(2017KJ229). Furthermore, we are grateful for Tianjin Key Laboratory on Technologies Enabling
Development of Clinical Therapeutics and Diagnostics (Theranostics).
Conflict of Interest: The authors declare that they have no conflict of interest.
Reference
Aier, I., Varadwaj, P. K. & Raj, U. (2016). Structural insights into conformational stability of
both wild-type and mutant EZH2 receptor. Scientific Reports 6.
Böhm, H.-J. (1992). The computer program LUDI: a new method for the de novo design of
enzyme inhibitors. Journal of Computer-Aided Molecular Design 6: 61-78.
Balmith, M. & Soliman, M. E. (2017). Non-active site mutations disturb the loop dynamics,
dimerization, viral budding and egress of VP40 of the Ebola virus. Mol Biosyst 13: 585-
597.
Barr, A. J., Ugochukwu, E., Lee, W. H., King, O. N., Filippakopoulos, P., Alfano, I., Savitsky, P.,
Burgess-Brown, N. A., Muller, S. & Knapp, S. (2009). Large-scale structural analysis of
the classical human protein tyrosine phosphatome. Cell 136: 352-363.
Berhanu, W. M. & Masunov, A. E. (2014). Full length amylin oligomer aggregation: insights
from molecular dynamics simulations and implications for design of aggregation
inhibitors. J Biomol Struct Dyn 32: 1651-1669.
25

Berman, H. M., Battistuz, T., Bhat, T., Bluhm, W. F., Bourne, P. E., Burkhardt, K., Feng, Z.,
Gilliland, G. L., Iype, L. & Jain, S. (2002). The protein data bank. Acta Crystallographica
Section D: Biological Crystallography 58: 899-907.
Buslaev, P., Gordeliy, V., Grudinin, S. & Gushchin, I. (2016). Principal Component Analysis of
Lipid Molecule Conformational Changes in Molecular Dynamics Simulations. J Chem
Theory Comput 12: 1019-1028.
Chatterjee, S., Khunti, K. & Davies, M. J. (2017). Type 2 diabetes. The Lancet 389: 2239-2251.
Cheng, X. C., Wang, R. L., Dong, Z. K., Li, J., Li, Y. Y. & Li, R. R. (2012). Design, synthesis
and evaluation of novel metalloproteinase inhibitors based on L-tyrosine scaffold. Bioorg
Med Chem 20: 5738-5744.
Cho, C. Y., Koo, S.-H., Wang, Y., Callaway, S., Hedrick, S., Mak, P. A., Orth, A. P., Peters, E.
C., Saez, E. & Montminy, M. (2006). Identification of the tyrosine phosphatase PTP-
MEG2 as an antagonist of hepatic insulin signaling. Cell Metabolism 3: 367-378.
Diller, D. J. & Merz Jr, K. M. (2001). High throughput docking for library design and library
prioritization. Proteins: Structure, Function, and Bioinformatics 43: 113-124.
Diller, D. J. & Merz, K. M. (2001). High throughput docking for library design and library
prioritization. Proteins: Structure, Function, and Bioinformatics 43: 113-124.
Fakhar, Z., Govender, T., Maguire, G. E. M., Lamichhane, G., Walker, R. C., Kruger, H. G. &
Honarparvar, B. (2017). Differential flap dynamics in l,d-transpeptidase2 from
mycobacterium tuberculosis revealed by molecular dynamics. Mol Biosyst 13: 1223-1234.
Foster, C. A. & West, A. H. (2017). Use of restrained molecular dynamics to predict the
conformations of phosphorylated receiver domains in two-component signaling systems.
Proteins 85: 155-176.
Gurzov, E. N., Stanley, W. J., Brodnicki, T. C. & Thomas, H. E. (2015). Protein tyrosine
phosphatases: molecular switches in metabolism and diabetes. Trends in Endocrinology
and Metabolism 26: 30-39.
He, R., Zeng, L. F., He, Y., Zhang, S. & Zhang, Z. Y. (2013). Small molecule tools for functional
interrogation of protein tyrosine phosphatases. The FEBS journal 280: 731-750.
Hendriks, W., Bourgonje, A., Leenders, W. & Pulido, R. (2018). Proteinaceous regulators and
inhibitors of protein tyrosine phosphatases. Molecules 23: 395.
Hendriks, W. J., Elson, A., Harroch, S., Pulido, R., Stoker, A. & den Hertog, J. (2013). Protein
tyrosine phosphatases in health and disease. The FEBS journal 280: 708-730.
Huynh, H., Wang, X., Li, W., Bottini, N., Williams, S., Nika, K., Ishihara, H., Godzik, A. &
Mustelin, T. (2003). Homotypic secretory vesicle fusion induced by the protein tyrosine
phosphatase MEG2 depends on polyphosphoinositides in T cells. The Journal of
Immunology 171: 6661-6671.
Irwin, J. J. & Shoichet, B. K. (2005). ZINC− a free database of commercially available
compounds for virtual screening. Journal of Chemical Information and Modeling 45: 177-
182.
Law, R., Barker, O., Barker, J. J., Hesterkamp, T., Godemann, R., Andersen, O., Fryatt, T.,
Courtney, S., Hallett, D. & Whittaker, M. (2009). The multiple roles of computational
chemistry in fragment-based drug design. Journal of computer-aided molecular design
23: 459-473.
Lee, H., Yi, J.-S., Lawan, A., Min, K. & Bennett, A. M. (2015). Mining the function of protein
tyrosine phosphatases in health and disease. Seminars in Cell & Developmental Biology,
Elsevier.
26

Li, H., Huang, J., Chen, L., Liu, X., Chen, T., Zhu, J., Lu, W., Shen, X., Li, J. & Hilgenfeld, R.
(2009). Identification of novel falcipain-2 inhibitors as potential antimalarial agents
through structure-based virtual screening. Journal of medicinal chemistry 52: 4936-4940.
Li, H. L., Ma, Y., Zheng, C. J., Jin, W. Y., Liu, W. S. & Wang, R. L. (2018). Exploring the effect
of D61G mutation on SHP2 cause gain of function activity by a molecular dynamics
study. J Biomol Struct Dyn 36: 3856-3868.
Liu, X., Dong, G., Zhang, J., Qi, J., Zheng, C., Zhou, Y., Zhu, J., Sheng, C. & , J. (2011).
Discovery of novel human acrosin inhibitors by virtual screening. Journal of computer-
aided molecular design 25: 977-985.
Loving, K., Alberts, I. & Sherman, W. (2010). Computational approaches for fragment-based and
de novo design. Current Topics in Medicinal Chemistry 10: 14-32.
Lukman, S., Lane, D. P. & Verma, C. S. (2013). Mapping the structural and dynamical features
of multiple p53 DNA binding domains: insights into loop 1 intrinsic dynamics. PLoS One
8: e80221.
Ma, Y., Wang, S. Q., Xu, W. R., Wang, R. L. & Chou, K. C. (2012). Design novel dual agonists
for treating type-2 diabetes by targeting peroxisome proliferator-activated receptors with
core hopping approach. PLoS One 7: e38546.
MacKerell, A. D., Bashford, D., Bellott, M., Dunbrack, R. L., Evanseck, J. D., Field, M. J.,
Fischer, S., Gao, J., Guo, H., Ha, S., Joseph-McCarthy, D., Kuchnir, L., Kuczera, K., Lau,
F. T., Mattos, C., Michnick, S., Ngo, T., Nguyen, D. T., Prodhom, B., Reiher, W. E.,
Roux, B., Schlenkrich, M., Smith, J. C., Stote, R., Straub, J., Watanabe, M., Wiorkiewicz-
Kuczera, J., Yin, D. & Karplus, M. (1998). All-atom empirical potential for molecular
modeling and dynamics studies of proteins. J Phys Chem B 102: 3586-3616.
Ndagi, U., Mhlongo, N. N. & Soliman, M. E. (2017). The impact of Thr91 mutation on c-Src
resistance to UM-164: molecular dynamics study revealed a new opportunity for drug
design. Mol Biosyst 13: 1157-1171.
Niu, Y. Z., Shi, D. F., Li, L. L., Guo, J. Y., Liu, H. X. & Yao, X. J. (2017). Revealing inhibition
difference between PFI-2 enantiomers against SETD7 by molecular dynamics
simulations, binding free energy calculations and unbinding pathway analysis. Scientific
Reports 7.
Piovesan, D., Minervini, G. & Tosatto, S. C. (2016). The RING 2.0 web server for high quality
residue interaction networks. Nucleic Acids Res 44: W367-374.
Qi, Y., Zhao, R., Cao, H., Sui, X., Krantz, S. B. & Zhao, Z. J. (2002). Purification and
characterization of protein tyrosine phosphatase PTP‐MEG2. Journal of Cellular
Biochemistry 86: 79-89.
Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., Amin, N.,
Schwikowski, B. & Ideker, T. (2003). Cytoscape: A software environment for integrated
models of biomolecular interaction networks. Genome Research 13: 2498-2504.
Shoichet, B. K. (2004). Virtual screening of chemical libraries. Nature 432: 862-865.
Sittel, F., Jain, A. & Stock, G. (2014). Principal component analysis of molecular dynamics: On
the use of Cartesian vs. internal coordinates. Journal of Chemical Physics 141.
Wang, D., Cheng, Z., Zhao, M., Jiao, C., Meng, Q., Pan, H., Xie, Y., Li, L., Zhu, Y. & Wang, W.
(2019). PTPN9 induces cell apoptosis by mitigating the activation of Stat3 and acts as a
tumor suppressor in colorectal cancer. Cancer management and research 11: 1309.
Wu, J.-w., Yin, L., Liu, Y.-q., Zhang, H., Xie, Y.-f., Wang, R.-l. & Zhao, G.-l. (2019). Synthesis,
biological evaluation and 3D-QSAR studies of 1, 2, 4-triazole-5-substituted carboxylic
27

acid bioisosteres as uric acid transporter 1 (URAT1) inhibitors for the treatment of
hyperuricemia associated with gout. Bioorganic & Medicinal Chemistry Letters 29: 383-
388.
Wu, J.-W., Zhang, H., Duan, Y.-Q., Dong, W.-L., Cheng, X.-C., Wang, S.-Q. & Wang, R.-L.
(2014). Design novel inhibitors for treating cancer by targeting Cdc25B catalytic domain
with de novo design. Combinatorial Chemistry & High Throughput Screening 17: 837-
847.
Zhang, L. S., Wang, S. Q., Xu, W. R., Wang, R. L. & Wang, J. F. (2012). Scaffold-based pan-
agonist design for the PPARalpha, PPARbeta and PPARgamma receptors. PLoS One 7:
e48453.
Zhang, S., Chen, L., Luo, Y., Gunawan, A., Lawrence, D. S. & Zhang, Z.-Y. (2009). Acquisition
of a potent and selective TC-PTP inhibitor via a stepwise fluorophore-tagged
combinatorial synthesis and screening strategy. Journal of the American Chemical Society
131: 13072-13079.
Zhang, S., Liu, S., Tao, R., Wei, D., Chen, L., Shen, W., Yu, Z. H., Wang, L., Jones, D. R.,
Dong, X. C. & Zhang, Z. Y. (2012). A highly selective and potent PTP-MEG2 inhibitor
with therapeutic potential for type 2 diabetes. J Am Chem Soc 134: 18116-18124.
Zhang, X., Wu, J., Liu, Y., Xie, Y., Liu, C., Wang, J. & Zhao, G. (2017). Facile synthetic
approaches to 1-thiocyclopropanecarboxylates. Phosphorus, Sulfur, and Silicon and the
Related Elements 192: 799-811.
Zhang, Z.-Y. (2001). Protein tyrosine phosphatases: prospects for therapeutics. Current Opinion
in Chemical Biology 5: 416-423.
28