Drug Design for mPGES-1 from Traditional Chinese Medicine Database: A Screening, Docking, QSAR, Molecular Dynamics, and Pharmacophore Mapping Study

and sharing with colleagues.

Other uses, including reproduction and distribution, or selling or

licensing copies, or posting to personal, institutional or third party

websites are prohibited.

In most cases authors are permitted to post their version of the

article (e.g. in Word or Tex form) to their personal website or

institutional repository. Authors requiring further information

regarding Elsevier’s archiving and manuscript policies are

encouraged to visit:

Drug design for mPGES-1 from traditional Chinese medicine database:

A screening, docking, QSAR, molecular dynamics, and pharmacophore

mapping study

Tung-Ti Chang

, Mao-Feng Sun

, Yung-Hao Wong

, Shun-Chieh Yang

, Kuan-Chung Chen

,

Hsin-Yi Chen

, Fuu-Jen Tsai

, Calvin Yu-Chian Chen

*

a

Laboratory of Computational and Systems Biology, School of Chinese Medicine, China Medical University, Taichung 40402, Taiwan

b_{Department of Chinese Pediatrics, China Medical University Hospital, Taiwan} c_{Department of Acupuncture, China Medical University Hospital, Taiwan} d

Department of Bioinformatics, Asia University, Taichung 41354, Taiwan

e

Department of Medical Genetics, Pediatrics and Medical Research, China Medical University Hospital and College of Chinese Medicine, China Medical University, Taichung 40402, Taiwan

f

Computational and Systems Biology, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

1. Introduction

Prostaglandins, products of prostanoids synthesis, are

auto-crines that are produced in various parts of human body.

Prostaglandin E

(PGE

), the most abundant subtype of

prosta-glandins, is the product from the action of prostaglandin E

synthases on prostaglandin H

(

Hara et al., 2010

). To date, there are

three known classes of prostaglandin E synthases, namely cytosolic

PGE synthase (cPGES) and two membrane-bound synthases,

mPGES-1 and mPGES-2 (

Hara et al., 2010

). The mPGES-1 protein,

unlike the other PGE synthases, is not constitutively expressed and

can be induced in response to inflammatory stimuli.

In the past, non-steroidal anti-inflammatory drugs (NSAIDs)

have been designed to target COX-1 and COX-1, the upstream

enzyme for producing prostaglandin H

. However, long-term

suppression of prostanoid biosynthesis by using COX-1 and COX-2

inhibitors can have severe side effects, including gastrointestinal

injury and renal irritation (

Koeberle and Werz, 2009

). Current

evidences suggest that suppression of mPGES-1 activity can be

considered as an alternative anti-inflammatory approach. Since

mPGES-1 is functionally coupled to COX-2 and is also being

responsible for excessive PGE

, its inhibitions are, therefore,

related to inflammation, pain, fever, atherosclerosis and

tumori-genesis.

Journal of the Taiwan Institute of Chemical Engineers 42 (2011) 580–591

A R T I C L E I N F O

Article history: Received 23 July 2010

Received in revised form 15 November 2010 Accepted 26 November 2010

Available online 26 January 2011

Keywords:

Microsomal prostaglandin E2synthase 1

(mPGES-1)

Traditional Chinese medicine (TCM) Docking

Molecular dynamics (MD)

Qualitative structure–activity relationship (QSAR)

A B S T R A C T

To search for new anti-inflammatory that can replace the current COX-1 and COX-2 inhibitors, virtual screening by molecular docking of traditional Chinese medicine (TCM) molecules into microsomal prostaglandin E2synthase (mPGES-1) glutathione binding site was performed. To compare the top

ranking derivatives with other mPGES-1 inhibitors, we constructed QSAR models using comparative molecular force field analysis (CoMFA) and comparative molecular similarity indices analysis (CoMSIA). The CoMFA model had a non-cross-validated coefficient (r2_{) and a cross-validated coefficient (q}2_{) of}

0.960 and 0.597. The r2_{and q}2_{for CoMSIA (S + H + D) was 0.931 and 0.719, respectively. The top three}

TCM derivatives all can map into the respective steric, hydrophobic and hydrogen bond donor force fields. The top ranking TCM molecules were taken for de novo design; the top three de novo products were further analyzed using molecular dynamics simulation and qualitative structure–activity relationship (QSAR) model. Derivative, 2-O-caffeoyl tartaric acid-Evo_2, glucogallin-Evo_1 and 4-O-feruloylquinic acid-Evo_7, all had conserved hydrogen bond networks to key residues Arg38 and Arg70 during the 20 ns molecular dynamics simulation. In addition, all derivative–protein complexes had total energy lower the control–protein complex. Combining the results from molecular dynamics simulation and CoMFA/ CoMSIA, we suggest 2-O-caffeoyl tartaric Evo_2, glucogallin-Evo_1 and 4-O-feruloylquinic acid-Evo_7 as potent mPGES-1 inhibitors.

ß2010 Taiwan Institute of Chemical Engineers. Published by Elsevier B.V. All rights reserved.

* Corresponding author at: Laboratory of Computational and Systems Biology, School of Chinese Medicine, China Medical University, Taichung 40402, Taiwan. Tel.: +886 4 22053366x3326/+1 617 353 7123.

E-mail addresses:[email protected],[email protected](C.-C. Chen).

Contents lists available at

ScienceDirect

Journal of the Taiwan Institute of Chemical Engineers

j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / j t i c e

1876-1070/$ – see front matter ß 2010 Taiwan Institute of Chemical Engineers. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.jtice.2010.11.009

To develop novel anti-inflammatory agents, TCM

Database@-Taiwan (

http://tcm.cmu.edu.tw

), the current world largest small

molecule database on traditional Chinese medicine, was employed

in docking. Natural compounds isolated from TCM are gaining

attention in the last few, and studies have been done already on

compounds, such as emodin or gallic acid and many more, to

investigate their potential effects to cellular process (

Chen, 2009d;

Chen et al., 2010b; Lin et al., 2010; Liu et al., 2010; Lo et al., 2010

). In

addition to screening of TCM database, we used both

structure-based and ligand-structure-based approaches to further validate and support

the bioactivity of the lead candidates. Both strategies have been

implemented in drug researches in the past, and we have

successfully used both techniques in designing viral,

anti-inflammatory, and anti-tumor agents (

Chang et al., 2010; Chen,

2008, 2009a,b,c,e,f, 2010a,b; Chen and Chen, 2007, 2010; Chen

et al., 2009a,b, 2010a; Huang et al., 2010a,b,c

).

2. Materials and methods

2.1. Docking

A total of 20,000 traditional Chinese medicine compounds

were

downloaded

from

TCM

Database@Taiwan

(

http://

tcm.cmu.edu.tw

) and docked into the glutathione binding

site of mPGES-1. All the ligands were pre-treated with force

field of CHARMm, and all the missing hydrogen were added. The

protein model used for docking was downloaded from Protein

Data Bank (PDB: 3DWW (

Jegerschold et al., 2008

)). The nature

substrate for mPGES-1, glutathione (g-

-glutamyl-

-cysteinyl-glycine), which co-crystallized with mPGES-1 by electron

crystallography, was used as the control molecule. The binding

location of glutathione found in the protein crystal was set as

the docking site.

Table 1

TCM docking results. Only the top 10 candidates and controls are shown.

Compound DS -PMF04 Jain Ludi1 Ludi2 Ludi 3

2-O-caffeoyl tartaric acid 215.079 144.63 5.42 870 698 789

Chicoric acid 206.092 177.22 5.00 915 694 826

Mumefural 201.985 136.75 8.50 1028 822 799

2-O-feruloyl tartaric acid 198.739 145.13 5.16 833 661 693

Rosmarinic acid 148.434 145.59 5.40 971 743 921 Quinic acid 143.961 106.56 1.49 557 487 475 Genipinic acid 142.772 110.31 2.53 530 467 651 Digallic acid 142.547 147.27 1.67 953 708 848 5-O-feruloylquinic acid 142.462 149.26 6.90 1011 784 865 4-O-feruloylquinic acid 140.488 142.14 7.45 946 746 878 Glutathione 66.787 127.28 6.39 564 447 412

DS: Dock Score; PMF: Potential of Mean Force.

Table 2

Docking results for de novo products. The top 5 candidates are shown.

Compound DS -PMF04 Jain Ludi1 Ludi2 Ludi3

2-O-caffeoyl tartaric acid-Evo_2 222.198 136.13 5.65 701 569 656

Glucogallin-Evo_1 1 169.762 158.84 6.12 1062 795 974

4-O-feruloylquinic acid-Evo_7 167.056 152.93 8.01 999 783 919 1-Caffeoylquinic acid-Evo_3 165.916 150.93 5.90 774 591 674 4-O-feruloylquinic acid-Evo_5 159.929 160.19 7.34 1071 822 950

Glutathione 66.787 127.28 6.39 564 447 412

DS: Dock Score; PMF: Potential of Mean Force.

Table 3

Structures of glutathione, the top three derivatives and the parental compounds of the top three derivatives. Glutathione 2-O-caffeoyl tartaric acid-Evo_2

[TD$INLINE] [TD$INLINE]

Glucogallin-Evo_1 4-O-feruloylquinic acid-Evo_7

[TD$INLINE] [TD$INLINE]

Table 4

The structure of the training set and the test set.

[TD$INLINE] Index R1 R2 R3 R4 R5 1 60_{-Cl Cl NC H NC} 2 H H Cl H F 3 60_{-Cl H Cl H F} 4 50_{-Cl H Cl H F} 5 60_{-Br H Cl H F} 6a 40_{-Br H Cl H F} 7 50_{-Br H Cl H F} 8 70_{-Br H Cl H F} 9 60_{-Br Br Cl H F} 10a ₆0_{-Br iPr Cl H F} 11 60_{-Br OMe Cl H F} 12 60_{-Br COMe Cl H F} 13 60_{-Br HO(CF} 3)2C Cl H F 14a 60_{-Br HO(CH} 3)2C Cl H F 15 60_{-Me Me Cl H F} 16 60_{-COMe COMe Cl H F} 17 60_-C(CH 3)2OH HO(CF3)2C Cl H F 18 60_{-Cl HO(CH} 3)2C NC H NC 19 60_{-Cl MeSO} 2 NC H NC 20 60_{-Cl NCCH} 2CH2CH2O NC H NC 21a 60_{-Cl [TD$INLINE] NC H NC} 22 60_{-Cl p-MeSO} 2C6H4 NC H NC 23 60_{-Cl MeOCH} 2CBBC NC H NC 24 60_{-Cl 3-PyridylCBBC NC H NC} 25 60_{-Cl 4-PyridylCBBC NC H NC} 26 60_{-Cl [TD$INLINE] NC H NC} 27a NC H NC 28 60_{-Et HO(CH} 3)2CCBBC NC H NC 29 [TD$INLINE] HO(CH3)2CCBBC NC H NC 30 60_{-Cl HO(CH} 3)2CCH2CH2 NC H NC 31a 60_{-Cl HO(CH} 3)2CCH2 NC H NC 32 [TD$INLINE] HO(CH3)2CCH2 NC H NC 33 60_-OCH2CH 2CH2CF3 HO(CH3)2CCH2 NC H NC 34 60_-OCH 2CH2CH(CH3)2 HO(CH3)2CCH2 NC H NC 35 60_-OCH 2CH(CH3)2 HO(CH3)2CCH2 NC H NC 36 [TD$INLINE] HO(CH3)2CCH2 NC H NC 37 [TD$INLINE] HO(CH3)2CCH2 NC F NC

a _{Test set compounds.}

T.-T. Chang et al. / Journal of the Taiwan Institute of Chemical Engineers 42 (2011) 580–591 582

The LigandFit program (

Venkatachalam et al., 2003

) within

Discovery Studio 2.5 was used to dock the TCM ingredients.

LigandFit is a receptor-rigid docking algorithm that utilizes Monte

Carlo simulation for generating ligand poses and a shape-matching

filter for comparing ligand poses with binding site shape.

Candidate ligand poses are then positioned into the binding site

and followed by rigid body energy minimization.

Dock Score (

Venkatachalam et al., 2003

), Potential of Mean

Force (PMF04) (

Muegge and Martin, 1999

), Jain (

Jain, 1996

) and

Ludi (

Bohm, 1998

) were calculated for predicting receptor–ligand

binding affinities. Dock Score was used as the primary scoring

function for ranking the poses of each ligand and for selecting the

poses. This scoring function is expressed as the sum of ligand

internal energy and receptor–ligand interaction energy. The other

scores, PMF04, Jain and Ludi, were calculated mainly for

references. PMF04 calculates pairwise interactions of all

inter-atomic pairs of a receptor–ligand system whereas Jain is based on

lipophilic interactions, polar attractive and repulsive interactions,

ligand entropy and solvation state of the protein and ligand. Ludi,

on the other hand, is the sum of contribution from ideal hydrogen

bond, ionic interaction, lipophilic interaction, and loss of internal

degrees of freedom, translational entropy and rotational entropy of

the ligand.

2.2. De novo design and Lipinski’s Rule of Five

The top ranking TCM compounds from docking were input into

De Novo Evolution protocol of Discovery Studio 2.5 for generating

derivatives. This protocol based on LUDI (

Bohm, 1992

), which

determines interaction sites suitable for hydrogen bond or

hydrophobic interaction, fits Ludi fragments best complementing

interaction sites and connects fitted fragments to a single molecule.

[()TD$FIG]

Fig. 1. Docking poses of (a) glutathione, (b) 2-O-caffeoyl tartaric acid-Evo_2, (c) glucogallin-Evo_1, and (d) 4-O-feruloylquinic acid-Evo_7 in mPGES-1 glutathione binding site.

Table 5

Statistical data for CoMFA and CoMSIA.

CoMFA CoMSIA Cross-validation Non-cross-validation

ONC q2 cv r2 SEE F ONC 6 S 6 0.446 0.792 0.306 15.238 q2 cv 0.597 H 6 0.733 0.956 0.141 86.480 SEE 0.135 D 6 0.167 0.432 0.506 3.043 r2 0.960 A 6 0.420 0.786 0.311 14.711 F 95.361 S + H 5 0.724 0.948 0.151 90.623 S + D 5 0.291 0.662 0.383 9.775 S + A 5 0.459 0.807 0.289 20.896 H + D 4 0.722 0.905 0.199 61.949 H + A 6 0.480 0.924 0.185 48.584 D + A 6 0.442 0.846 0.263 22.053 S + H + D 0.719 0.931 0.176 54.319 S + H + A 5 0.642 0.915 0.192 54.043 S + D + A 5 0.501 0.836 0.267 25.455 H + D + A 5 0.654 0.930 0.175 65.965 S + H + D + A 6 0.630 0.940 0.164 62.960

ONC: optimal number of components; SEE: standard error of estimate; F: F-test value; PLS: partial least squares; S: steric; H: hydrophobic; D: hydrogen bond donor; A: hydrogen bond acceptor.

The evolution mode for generating derivatives was set to full

evolution that builds molecules in an evolutionary fashion. The

maximum number of generations, survivors and population size

were set to 2, 3 and 10, respectively. Ludi scoring function (Ludi 3)

was used to evaluate the fitness of the newly generated derivatives.

The generated derivatives, before re-docked back to protein to

evaluate protein–ligand interactions, were screened with the

Lipinski’s Rule of Five (

Lipinski et al., 2001

). The derivatives that

exceed the number of hydrogen bond donor and acceptors, the

upper molecular weight limit, and the octanol–water coefficient

were excluded from further analysis. An ADMET analysis was also

performed to rule out derivatives that poses potential toxicity to

human body. This is done in Discovery Studio 2.5.

2.3. CoMFA and CoMSIA

The structures of 37 phenanthrene imidazole based compounds

used in this study were compiled from

Giroux et al. (2009)

. The

inhibitory potencies of these phenanthrene imidazoles were

assessed in vitro using recombinant human mPGES-1 enzyme,

and the reported IC

values were means of at least two

experiments (

Mancini et al., 2001

). The compounds were further

divided into a training set of 31 compounds and a test set of 6

compounds. The selection for training set and test set was totally

random. Alignment of the training set molecules was performed

using atom-fit module of SYBYL 8.0.

To build QSAR model using comparative force field analysis

(CoMFA), the steric and electrostatic field descriptors were

calculated using Lennard–Jones potential and Coulombic potential,

respectively. For comparative similarity indices analysis (CoMSIA),

the steric, electrostatic, hydrophobic and hydrogen bond donor

and acceptors were calculated with Gaussian function, in contrast

to distance-dependent dielectric method used in CoMFA. Partial

least squares (PLS) regression was utilized to analyze the 3D-QSAR

models, and all CoMFA and CoMSIA fields were taken as

independent variables.

2.4. Molecular dynamics simulation

Selected TCM-mPGES-1 complexes were taken for molecular

dynamics simulation. Each complex was solvated in a cubic water

box before energetically minimized using 500 steps of Steepest

Descent and 500 steps of Conjugate Gradient method. The system

was then heated from 50 K to 310 K without constraint for 50 ps. The

equilibration step was conducted for 200 ps without constraint. The

final production step was conducted for 20 ns in NVT ensemble with

snapshots save every 2.5 ps. The time step was set to 1 fs. Particle

mesh Ewald (PME) was used for electrostatic calculation. Hydrogen

bond frequency, energy trajectory, and hydrogen bond distance

were calculated for analyzing the ligand–protein system.

3. Results and discussions

3.1. Docking and de novo design

The binding affinities of TCM ligands to the receptor were

predicted based on scoring functions. We used Dock Score as the

[()TD$FIG]

Fig. 2. Predicted pIC50vs. experiment pIC50for (a) CoMFA and (b) CoMSIA.

Table 6

Predicted activity of the training set and the test set.

Compound pIC50 CoMFA CoMSIA

Predicted Residual Predicted Residual 1 9.00 8.79 0.21 8.75 0.25 2 7.44 7.80 0.36 7.63 0.19 3 8.30 8.49 0.19 8.19 0.11 4 7.96 7.72 0.24 7.68 0.29 5 8.70 8.69 0.01 8.55 0.15 6a _{6.96 7.75} 0.79 7.60 0.64 7 7.60 7.53 0.07 7.76 0.16 8 8.10 7.98 0.12 8.15 0.05 9 8.52 8.57 0.04 8.79 0.27 10a 8.40 8.31 0.09 8.69 0.29 11 8.22 8.23 0.01 8.42 0.20 14a 8.10 8.01 0.09 8.14 0.04 15 9.00 8.76 0.24 8.93 0.07 17 8.15 8.73 0.58 8.35 0.20 18 8.15 8.27 0.12 8.02 0.14 19 6.91 6.92 0.01 6.78 0.13 20 7.35 7.37 0.02 7.45 0.10 21a 8.40 8.50 0.10 8.52 0.12 22 8.05 8.10 0.05 8.16 0.11 24 9.00 8.99 0.01 8.98 0.02 25 9.00 8.80 0.21 8.88 0.12 26 8.70 8.69 0.01 8.61 0.09 27a 8.70 8.71 0.01 8.79 0.09 28 9.00 8.97 0.03 8.96 0.04 30 9.00 9.10 0.10 8.94 0.06 31a 9.00 9.02 0.02 8.95 0.05 32 9.00 9.01 0.01 9.15 0.15 33 9.00 9.06 0.06 9.02 0.02 34 9.00 9.09 0.09 9.01 0.01 35 9.00 8.96 0.04 9.04 0.04 36 8.70 8.66 0.04 8.69 0.01 37 8.52 8.60 0.08 8.52 0.00 a

Test set compounds.

T.-T. Chang et al. / Journal of the Taiwan Institute of Chemical Engineers 42 (2011) 580–591 584

primary scoring function for ranking the ligands as this algorithm

has been shown to correlate positively with actual bioactivity. TCM

molecules ranked top in the docking are shown in

Table 1

; the top

10 compounds all have Dock Score higher than the control

glutathione.

The top-ranking TCM derivatives, which all passed the ADMET

screen, are shown in

Table 2

. As shown in

Table 2

, derivative of

2-O-caffeoyl tartaric acid derivative, 4-O-feruloylquinic acid

deriva-tive and glucogallin all have substantial increase in binding

affinities as compared to the parental compounds. The structures

of the top three derivatives and glutathione are shown in

Table 3

.

[()TD$FIG]

Fig. 3. The contour maps of CoMFA and CoMISA with (a) glutathione, (b) 2-O-caffeoyl tartaric acid-Evo_2, (c) glucogallin-Evo_1, and (d) 4-O-feruloylquinic acid-Evo_7 in the mPGES-1 binding site. Green and yellow contours indicate where steric groups would be favored or disfavored. Cyan and purple contours indicate where would be hydrophobic or hydrophilic. (For interpretation of the references to color in the figure caption, the reader is referred to the web version of the article.)

[()TD$FIG]

Fig. 4. The all atoms RMSD (top) and ligand RMSD (bottom) of mPGES-1–ligand complexes.

[()TD$FIG]

Fig. 5. The total energy for mPGES-1 complexes during 20 ns simulation. T.-T. Chang et al. / Journal of the Taiwan Institute of Chemical Engineers 42 (2011) 580–591 585

The docking conformations of these three derivatives are shown in

Fig. 1

. A conserved hydrogen bond network is formed between the

mPGES-1 binding site residues and all the TCM molecules and the

control. The key interacting residues are Arg38, Arg70, Arg73,

Arg110 and Arg126 (

Fig. 1

).

3.2. CoMFA/CoMSIA

The structure and the experimental IC

of the training set and

the test set are shown in

Table 4

. The core atoms used in alignment

are shown in

Table 4

as well. The results for CoMFA and CoMSIA are

[()TD$FIG]

Fig. 6. Time dependent hydrogen bond distances between glutathione and mPGES-1 residues.

Table 7

Hydrogen bond statistics summary for glutathione.

Ligand Amino acid Max. distance Min. distance Average distance H-bond occupancy

O16 A:ARG110:HH11 3.81 2.01 3.04 5.33% O16 A:ARG110:HH12 3.50 1.83 2.71 16.06% O16 A:ARG110:HH22 4.20 2.27 3.09 0.48% N7 A:ARG126:HH11 5.25 1.94 3.66 0.33% N7 A:ARG126:HH12 4.86 2.24 3.13 0.38% O12 A:ARG126:HH12 5.01 2.17 3.60 1.38% O5 A:ARG126:HH12 3.15 1.70 2.33 76.91% O12 A:ARG126:HH21 5.61 2.36 3.24 0.05% O12 A:ARG126:HH22 4.53 1.84 2.78 19.60% O5 A:ARG126:HH22 2.96 1.58 2.19 92.25% O19 A:ARG70:HH11 4.14 1.60 2.21 83.70% O20 A:ARG70:HH11 4.28 1.59 3.15 11.40% O19 A:ARG73:HE 6.37 1.97 2.83 25.40% O20 A:ASN74:HD22 5.15 2.39 3.60 0.04% H37 A:ASN74:OD1 5.44 1.73 2.88 18.35% H26 A:GLU77:OE2 4.19 2.18 3.44 0.99% H27 A:GLU77:OE1 4.51 1.96 2.70 23.90% H32 A:GLU77:OE1 7.98 1.90 3.57 38.46% H32 A:GLU77:OE2 7.09 1.75 3.15 52.56% H36 A:GLN134:OE1 5.93 2.44 4.09 0.01% H27 A:HIS113:NE2 4.90 1.86 3.04 12.25% H28 A:HIS113:NE2 4.41 2.00 3.29 3.61% H36 A:TYR130:O 5.55 1.70 4.40 7.84% O12 B:ARG38:HH12 4.78 1.67 2.55 56.71% S11 B:ARG38:HH12 6.21 1.78 4.61 11.05% O12 B:ARG38:HH22 4.61 2.33 3.47 0.23% S11 B:CYS68:HG 8.71 2.46 6.47 0.01% H32 B:CYS68:SG 8.47 1.95 6.35 0.81% H26 B:HIS72:NE2 4.92 2.29 3.06 0.39% H32 B:HIS72:NE2 5.90 1.88 3.57 4.80% O16 B:TYR28:HH 5.40 2.31 4.47 0.05% H36 B:TYR28:OH 5.50 1.73 2.48 61.39%

T.-T. Chang et al. / Journal of the Taiwan Institute of Chemical Engineers 42 (2011) 580–591 586

shown in

Table 5

. The CoMFA model has a non-cross-validated

coefficient of 0.960 and a cross-validated coefficient of 0.597

(

Fig. 2

). As for the CoMSIA model, we selected the model containing

steric, hydrophobic and hydrogen bond donor descriptor. For

CoMSIA, the non-validated coefficient is 0.931 and the

cross-validation coefficient is 0.719 (

Fig. 2

). The activities predicted by

CoMFA and CoMSIA, for both the training set and the test set, are

shown in

Table 6

. Most predicted pIC

only deviate less than

one-fifth of log order from the experimental pIC

. These values suggest

that the CoMFA and CoMSIA models are reliable.

[()TD$FIG]

Fig. 7. Time dependent hydrogen bond distances between 2-O-caffeoyl tartaric acid-Evo_2 with mPGES-1 residues.

[()TD$FIG]

Fig. 8. Time dependent hydrogen bond distances between glucogallin-Evo_1 and mPGES-1 residues.

[()TD$FIG]

Fig. 9. Time dependent hydrogen bond distances between 4-O-feruloylquinic acid-Evo_7 and mPGES-1 residues.

Table 8

Hydrogen bond statistics summary for 2-O-caffeoyl tartaric acid-Evo_2.

Ligand Amino acid Max. distance Min. distance Average distance H-bond occupancy

H27 A:HIS72:O 6.338 1.617 2.72 24.04% H7 A:TYR117:OH 3.994 1.83 3.48 0.15% O1 A:ARG126:HE 3.964 1.818 2.99 11.39% O10 A:ARG70:HH11 4.068 1.868 3.09 1.29% O10 A:ARG70:HH12 3.836 1.808 2.59 39.41% O11 A:ARG70:HH12 4.334 1.649 2.29 82.53% O6 A:ARG70:HH12 3.907 1.78 3.9071 11.41% O10 A:ARG73:HH11 2.558 1.543 1.97 99.96% O11 A:ARG73:HH11 3.172 1.677 2.07 94.36% O11 A:ARG73:HH12 3.978 2.138 2.70 11.70% O6 A:TYR117:HH 4.136 1.765 2.73 12.43% O12 B:ARG38:HH11 4.262 1.815 3.37 9.23% O14 B:ARG38:HH11 4.321 1.692 3.43 17.70% O1 B:ARG38:HH12 3.697 1.76 2.46 58.15% O11 B:ARG38:HH12 5.923 1.821 4.70 4.74% O12 B:ARG38:HH12 3.37 1.636 2.46 55.35% O3 B:ARG38:HH12 2.879 1.588 1.92 99.71% O3 B:ARG38:HH22 4.091 1.633 2.40 59.25% O26 B:ARG73:HH12 4.551 2.215 3.40 5.61% O11 B:CYS68:HG 8.382 2.119 5.75 2.65% O14 B:CYS68:HG 6.206 1.651 2.52 72.66% O1 B:LYS42:HZ1 5.493 1.633 3.69 6.36% O3 B:LYS42:HZ1 3.87 1.654 3.15 6.91% O1 B:LYS42:HZ2 4.392 1.591 2.68 57.26% O3 B:LYS42:HZ2 3.737 1.6 2.39 56.01% O1 B:LYS42:HZ3 3.915 2.139 3.18 0.20% O3 B:LYS42:HZ3 3.81 1.552 2.13 83.90% O1 A:ARG126:HH21 4.056 2.392 3.10 0.15% O14 B:ARG38:HH12 4.263 2.441 3.49 0.03% O26 B:ARG73:HH11 4.93 1.71 2.37 77.23% O29 B:ARG73:HH22 4.31 2.453 3.44 0.05% O11 A:ARG70:HH22 5.298 2.244 3.70 0.34% H27 A:MET76:SD 5.739 2.136 2.92 13.66% O11 A:ARG70:HH11 5.511 2.095 3.70 0.13% O12 B:ARG38:HH22 4.715 2.391 3.38 0.08%

T.-T. Chang et al. / Journal of the Taiwan Institute of Chemical Engineers 42 (2011) 580–591 588

To compare our TCM derivatives with experimentally tested

mPGES-1 inhibitors, we mapped the TCM derivatives into the

CoMFA/CoMSIA contours. We only show the contour map from

CoMSIA as it has covered the force field in CoMFA. As illustrated in

Fig. 3

, in addition to glutathione, all TCM derivatives map well into

the CoMSIA contours. The hydrophilic contours are indicated by

purple color, and the hydrophobic contours are shown in cyan. The

green and yellow contours represent where adding steric group

would be favored or disfavored. The top three TCM derivatives all

contain hydrophilic groups and can form hydrogen bond

interac-tions to Arg70 and Arg38. These hydrogen bond interacinterac-tions are all

mapped onto the purple contours, which are in well agreement

with the QSAR model. The contour maps can also be superimposed

onto to corresponding residues of the mPGES-1 binding site. This

further supports the docking result of the top TCM derivatives.

3.3. Molecular dynamics simulation

The top TCM derivatives and the control were selected for

molecular dynamics simulation. The whole complex RMSD (top)

Table 9

Hydrogen bond statistics summary for Glucogallin-Evo_1.

Ligand Amino acid Max. distance Min. distance Average distance H-bond occupancy

O39 A:ARG110:HH12 3.85 1.97 3.14 0.08% O39 A:ARG110:HH21 3.89 2.32 3.13 0.01% O39 A:ARG110:HH22 3.61 1.79 2.81 3.86% O24 A:ARG70:HH11 4.31 1.79 2.46 66.54% O20 A:ARG70:HH12 3.72 1.71 2.09 96.45% O22 A:ARG70:HH12 3.73 1.78 3.04 2.48% O24 A:ARG70:HH12 4.01 2.20 2.73 7.93% O20 A:ARG70:HH22 4.34 1.74 2.25 91.29% O22 A:ARG70:HH22 5.74 2.50 4.27 0.01% O24 A:ARG73:HH11 3.99 2.29 2.79 1.13% O41 A:ASN74:HD22 4.96 2.18 4.13 0.23% O36 A:GLN134:HE22 5.32 2.01 4.05 2.98% O20 A:TYR117:HH 5.48 2.24 3.40 1.59% O20 B:ARG38:HH12 2.74 1.66 2.03 99.94% O21 B:ARG38:HH12 3.52 1.96 2.59 34.83% O22 B:ARG38:HH12 4.12 2.17 2.65 9.48% O20 B:ARG38:HH22 3.94 1.61 2.20 91.53% O21 B:ARG38:HH22 3.56 1.60 2.21 91.88% H37 B:ALA31:O 4.17 2.05 3.49 1.00% H40 A:GLU77:OE1 5.01 2.19 4.32 0.25% H40 B:TYR28:OH 5.34 1.90 4.42 3.00% H42 A:ASN74:OD1 4.53 1.63 2.02 93.35% H42 A:GLU77:OE1 2.95 1.74 2.19 97.49% H42 A:GLU77:OE2 4.00 2.41 3.42 0.08% H8 A:GLU77:OE2 4.58 1.84 3.69 1.51% O21 A:ARG70:HH22 4.37 2.18 3.08 1.10% O21 B:CYS68:HG 5.74 2.00 4.06 1.94% H37 B:TYR28:OH 5.64 2.48 4.53 0.01% O7 B:HIS72:HD1 6.49 1.70 2.78 68.90% H37 A:TYR130:O 5.93 2.28 3.45 0.69% Table 10

Hydrogen bond statistics summary for 4-O-feruloylquinic acid-Evo_7.

Ligand Amino acid Max. distance Min. distance Average distance H-bond occupancy

H11 B:TYR28:OH 3.68 1.64 2.67 24.65% H40 A:GLU77:OE2 3.94 2.44 3.00 0.10% H44 B:HIS72:NE2 6.27 2.07 5.57 0.28% H46 A:HIS113:NE2 3.05 1.81 2.33 83.70% O10 A:ARG110:HH12 3.16 1.53 1.86 97.19% O10 A:ARG110:HH22 2.90 1.78 2.32 87.91% O43 A:ARG126:HH11 6.92 2.18 3.15 7.41% O45 A:ARG126:HH11 5.66 1.84 4.48 17.06% O21 A:ARG126:HH12 4.21 1.84 2.61 59.19% O43 A:ARG126:HH12 6.07 1.59 2.46 80.46% O45 A:ARG126:HH12 4.67 1.73 3.54 18.75% O45 A:ARG126:HH21 6.84 1.95 3.39 10.60% O21 A:ARG126:HH22 4.43 1.65 2.49 58.54% O45 A:ARG126:HH22 6.24 1.72 2.74 68.58% O30 A:ARG70:HH11 3.47 2.09 2.69 10.54% O31 A:ARG70:HH11 2.44 1.61 1.86 100.00% O30 A:ARG70:HH12 3.49 2.18 2.90 1.55% O31 A:ARG70:HH12 3.90 2.15 3.20 0.05% O31 A:TYR117:HH 4.45 2.26 3.60 0.04% O5 A:GLN134:HE21 4.20 2.25 3.06 0.04% O21 B:ARG38:HH11 4.98 2.25 3.80 0.44% O30 B:ARG38:HH12 3.78 1.61 2.09 92.28% O30 B:ARG38:HH22 2.54 1.62 1.92 99.99% O32 B:ARG38:HH22 3.98 1.91 3.34 0.25% O10 B:TYR28:HH 4.20 1.79 2.75 10.73% O5 B:TYR28:HH 5.47 2.00 3.48 1.00%

and the ligand RMSD (bottom) are shown in

Fig. 4

. All the ligand–

protein complexes have reached equilibration during the

simula-tion period. However, the glutathione–mPGES-1 complex has an

increase in whole molecule and ligand RMSD value at 18 ns. In

physiological state, binding of glutathione to mPGES-1 induces a

conformational change in binding site residues, allowing for

effective conversion of PGH

to PGE

. We, therefore, attribute this

increase ligand RMSD to changes in protein conformation. The

energy trajectories of for all ligand–mPGES-1 complexes are shown

in

Fig. 5

. All derivative–protein complexes have total energy lower

than the control.

For the control, the most significant fluctuation in hydrogen

bond distances appears in the beginning of the simulation run and

at 18 ns (

Fig. 6

). The shift in interaction distance may be because

the system is yet to adapt equilibration. At 18 ns, more hydrogen

bonds are formed between the control and Arg38; on the contrary,

interactions to Glu77 decrease (

Table 7

and

Supplementary video

1

). As illustrated earlier, this is due to changes in both protein and

ligand conformation.

For all the TCM derivatives, 2-O-caffeoyl tartaric acid-Evo_2,

glucogallin-Evo_1 and 4-O-feruloylquinic acid-Evo_7, all have very

stable interactions to binding site residues, with no obvious change

in binding distances (

Figs. 7–9

). The carboxyl and carbonyl group of

2-O-caffeoyl tartaric acid-Evo_2 have H-bond interaction with

residue Cys68 or Arg38 before the 4 ns of simulation (

Supplemen-tary video 2

). Afterwards, at the following stable stage of

simulation, the interaction was formed between 2-O-caffeoyl

tartaric acid-Evo_2 and Arg38 shifts to between 2-O-caffeoyl

tartaric acid-Evo_2 and Arg73 (

Table 8

and

Supplementary video

2

). The interaction between glucogallin-Evo_1 hydroxyl group and

Arg38 and Arg70 persist throughout the simulation time (

Fig. 8

and

Supplementary video 3

). In addition, after equilibration hydrogen

bond interactions are also found between glucogallin-Evo_1 and

Asn74 and His7 (

Table 9

), which further strengthens the

interaction between the derivative and binding site residues

(

Supplementary video 3

). Similar to other TCM derivatives,

4-O-feruloylquinic acid-Evo_7 also established stable H-bonds with

Arg38 and Arg70 as observed from hydrogen bond statistics (

Table

10

), hydrogen bond distance (

Fig. 9

) and

Supplementary video 4

.

Overall, the TCM derivatives all have stable and continuous

hydrogen bond interaction to binding site residues; these

interactions keep the derivative firmly bonded to mPGES-1

glutathione binding site.

4. Conclusion

To identify new anti-inflammatory molecules, docking of TCM

molecules into mPGES-1 glutathione bind site was performed.

Subsequent de novo design generated 2-O-caffeoyl tartaric acid_2,

glucogallin_1 and 4-O-feruloylquinic acid_7 that serve great

potential as next generation of anti-inflammatory compounds.

These derivatives were selected for further analyses using QSAR

and molecular dynamics simulation. Both the CoMFA and CoMSIA

models investigated had high non-cross-validation coefficient of

r

of 0.960 and 0.719 and high cross-validation coefficient q

of

0.597 and 0.719. The top three derivatives were able to map onto

the steric, hydrophobic and hydrogen bond donor fields of

CoMSIA, meeting the requirements of being effective mPGES-1

inhibitor. Analyses of molecular dynamics simulation showed

that the TCM derivatives have continuous and stable hydrogen

bond interactions to key binding site residues, Arg38 and Arg70. In

addition, extensive hydrogen bond networks were formed

between the TCM derivatives and the mPGES-1 binding site

residues. This event suggests that the top three TCM derivatives

can bind strongly to mPGES-1 residues, making them very potent

candidates.

Acknowledgements

The research was supported by grants from the National Science

Council of Taiwan (NSC 99-2221-E-039-013), China Medical

University (CMU98-TCM, CMU99-TCM, CMU99-S-02) and Asia

University (CMU98-ASIA-09). This study is also supported in part

by Taiwan Department of Health Clinical Trial and Research Center

of Excellence (DOH99-TD-B-111-004) and Taiwan Department of

Health Cancer Research Center of Excellence

(DOH99-TD-C-111-005). We are grateful to the National Center of High-performance

Computing for computer time and facilities.

Appendix A. Supplementary data

Supplementary data associated with this article can be found, in

the online version, at

doi:10.1016/j.jtice.2010.11.009

.

References

Bohm, H. J., ‘‘The Computer Program Ludi: A New Method for the De Novo Design of Enzyme Inhibitors,’’ J. Comput. Aided Mol. Des., 6, 61 (1992).

Bohm, H. J., ‘‘Prediction of Binding Constants of Protein Ligands: A Fast Method for the Prioritization of Hits Obtained from De Novo Design or 3d Database Search Programs,’’ J. Comput. Aided Mol. Des., 12, 309 (1998).

Chang, T. T., H. J. Huang, K. J. Lee, H. W. Yu, H. Y. Chen, F. J. Tsai, M. F. Sun, and C. Y. C. Chen, ‘‘Key Features for Designing Phosphodiesterase-5 Inhibitors,’’ J. Biomol. Struct. Dyn., 28, 309 (2010).

Chen, C. Y. C., ‘‘Discovery of Novel Inhibitors for C-Met by Virtual Screening and Pharmacophore Analysis,’’ J. Chin. Inst. Chem. Engrs., 39, 617 (2008).

Chen, C. Y. C., ‘‘Chemoinformatics and Pharmacoinformatics Approach for Exploring the Gaba-a Agonist from Chinese Herb Suanzaoren,’’ J. Taiwan Inst. Chem. Engrs., 40, 36 (2009a).

Chen, C. Y. C., ‘‘Computational Screening and Design of Traditional Chinese Medicine (Tcm) to Block Phosphodiesterase-5,’’ J. Mol. Graph. Model., 28, 261 (2009b). Chen, C. Y. C., ‘‘De Novo Design of Novel Selective Cox-2 Inhibitors: From Virtual

Screening to Pharmacophore Analysis,’’ J. Taiwan Inst. Chem. Engrs., 40, 55 (2009c). Chen, C. Y. C., ‘‘Magnolol Encapsulated by Different Acyl Chain Length of Liposomes on Inhibiting Proliferation of Smooth Muscle Cells,’’ J. Taiwan Inst. Chem. Engrs., 40, 380 (2009d).

Chen, C. Y. C., ‘‘Pharmacoinformatics Approach for Mpges-1 in Anti-Inflammation by 3d-Qsar Pharmacophore Mapping,’’ J. Taiwan Inst. Chem. Engrs., 40, 155 (2009e). Chen, C. Y. C., ‘‘Weighted Equation and Rules—A Novel Concept for Evaluating Protein–

Ligand Interaction,’’ J. Biomol. Struct. Dyn., 27, 271 (2009f).

Chen, C. Y. C., ‘‘Bioinformatics, Chemoinformatics, and Pharmainformatics Analysis of Her2/Hsp90 Dual-Targeted Inhibitors,’’ J. Taiwan Inst. Chem. Engrs., 41, 143 (2010a). Chen, C. Y. C., ‘‘Virtual Screening and Drug Design for Pde-5 Receptor from Traditional

Chinese Medicine Database,’’ J. Biomol. Struct. Dyn., 27, 627 (2010b).

Chen, C. Y., Y. H. Chang, D. T. Bau, H. J. Huang, F. J. Tsai, C. H. Tsai, and C. Y. C. Chen, ‘‘Discovery of Potent Inhibitors for Phosphodiesterase 5 by Virtual Screening and Pharmacophore Analysis,’’ Acta Pharmacol. Sin., 30, 1186 (2009a).

Chen, C. Y., Y. H. Chang, D. T. Bau, H. J. Huang, F. J. Tsai, C. H. Tsai, and C. Y. C. Chen, ‘‘Ligand-based Dual Target Drug Design for H1N1: Swine Flu-a Preliminary First Study,’’ J. Biomol. Struct. Dyn., 27, 171 (2009b).

Chen, Y. C. and K. T. Chen, ‘‘Novel Selective Inhibitors of Hydroxyxanthone Derivatives for Human Cyclooxygenase-2,’’ Acta Pharmacol. Sin., 28, 2027 (2007). Chen, C. Y. and C. Y. C. Chen, ‘‘Insights into Designing the Dual-Targeted Her2/Hsp90

Inhibitors,’’ J. Mol. Graph. Model., 29, 21 (2010).

Chen, C. Y., H. J. Huang, F. J. Tsai, and C. Y. C. Chen, ‘‘Drug Design for Influenza a Virus Subtype H1N1,’’ J. Taiwan Inst. Chem. Engrs., 41, 8 (2010a).

Chen, C. Y., P. L. Kuo, Y. H. Chen, J. C. Huang, M. L. Ho, R. J. Lin, J. S. Chang, and H. M. Wang, ‘‘Tyrosinase Inhibition, Free Radical Scavenging, Antimicroorganism and Antican-cer Proliferation Activities of Sapindus Mukorossi Extracts,’’ J. Taiwan Inst. Chem. Engrs., 41, 129 (2010b).

Giroux, A., L. Boulet, C. Brideau, A. Chau, D. Claveau, B. Cote, D. Ethier, R. Frenette, M. Gagnon, J. Guay, S. Guiral, J. Mancini, E. Martins, F. Masse, N. Methot, D. Riendeau, J. Rubin, D. Xu, H. Yu, Y. Ducharme, and R. W. Friesen, ‘‘Discovery of Disubstituted Phenanthrene Imidazoles as Potent, Selective and Orally Active Mpges-1 Inhibitors,’’ Bioorg. Med. Chem. Lett., 19, 5837 (2009).

Hara, S., D. Kamei, Y. Sasaki, A. Tanemoto, Y. Nakatani, and M. Murakami, ‘‘Prosta-glandin E Synthases: Understanding Their Pathophysiological Roles through Mouse Genetic Models,’’ Biochimie, 92, 651 (2010).

Huang, H. J., C. Y. Chen, H. Y. Chen, F. J. Tsai, and C. Y. C. Chen, ‘‘Computational Screening and Qsar Analysis for Design of Amp-Activated Protein Kinase Agonist,’’ J. Taiwan Inst. Chem. Engrs., 41, 352 (2010a).

Huang, H. J., K. J. Lee, H. W. Yu, C. Y. Chen, C. H. Hsu, H. Y. Chen, F. J. Tsai, and C. Y. C. Chen, ‘‘Structure-based and Ligand-based Drug Design for Her 2 Receptor,’’ J. Biomol. Struct. Dyn., 28, 23 (2010b).

Huang, H. J., K. J. Lee, H. W. Yu, H. Y. Chen, F. J. Tsai, and C. Y. C. Chen, ‘‘A Novel Strategy for Designing the Selective Ppar Agonist by The Sum of Activity Model,’’ J. Biomol. Struct. Dyn., 28, 187 (2010c).

T.-T. Chang et al. / Journal of the Taiwan Institute of Chemical Engineers 42 (2011) 580–591 590

Jain, A. N., ‘‘Scoring Noncovalent Protein–Ligand Interactions: A Continuous Differen-tiable Function Tuned to Compute Binding Affinities,’’ J. Comput. Aided Mol. Des., 10, 427 (1996).

Jegerschold, C., S. C. Pawelzik, P. Purhonen, P. Bhakat, K. R. Gheorghe, N. Gyobu, K. Mitsuoka, R. Morgenstern, P. J. Jakobsson, and H. Hebert, ‘‘Structural Basis for Induced Formation of the Inflammatory Mediator Prostaglandin E2,’’ Proc. Natl. Acad. Sci. U.S.A., 105, 11110 (2008).

Koeberle, A. and O. Werz, ‘‘Inhibitors of the Microsomal Prostaglandin E-2 Synthase-1 as Alternative to Non Steroidal Anti-Inflammatory Drugs (Nsaids)—A Critical Review,’’ Curr. Med. Chem., 16, 4274 (2009).

Lin, M. L., Y. C. Lu, J. G. Chung, S. G. Wang, H. T. Lin, S. E. Kang, C. H. Tang, J. L. Ko, and S. S. Chen, ‘‘Down-Regulation of Mmp-2 through the P38 Mapk-Nf-Kappa B-Depen-dent Pathway by Aloe-Emodin Leads to Inhibition of Nasopharyngeal Carcinoma Cell Invasion,’’ Mol. Carcinog., 49, 783 (2010).

Lipinski, C. A., F. Lombardo, B. W. Dominy, and P. J. Feeney, ‘‘Experimental and Computational Approaches to Estimate Solubility and Permeability in Drug Dis-covery and Development Settings,’’ Adv. Drug Deliv. Rev., 46, 3 (2001).

Liu, J. F., W. H. Yang, Y. C. Fong, S. C. Kuo, C. S. Chang, and C. H. Tang, ‘‘Bfpp, a Phloroglucinol Derivative. Induces Cell Apoptosis in Human Chondrosarcoma Cells through Endoplasmic Reticulum Stress,’’ Biochem. Pharmacol., 79, 1410 (2010).

Lo, C., T. Y. Lai, J. H. Yang, J. S. Yang, Y. S. Ma, S. W. Weng, Y. Y. Chen, J. G. Lin, and J. G. Chung, ‘‘Gallic Acid Induces Apoptosis in A375.S2 Human Melanoma Cells through Caspase-Dependent and -Independent Pathways,’’ Int. J. Oncol., 37, 377 (2010).

Mancini, J. A., K. Blood, J. Guay, R. Gordon, D. Claveau, C. C. Chan, and D. Riendeau, ‘‘Cloning, Expression, and up-Regulation of Inducible Rat Prostaglandin E Synthase during Lipopolysaccharide-induced Pyresis and Adjuvant-induced Arthritis,’’ J. Biol. Chem., 276, 4469 (2001).

Muegge, I. and Y. C. Martin, ‘‘A General and Fast Scoring Function for Protein–Ligand Interactions: A Simplified Potential Approach,’’ J. Med. Chem., 42, 791 (1999). Venkatachalam, C. M., X. Jiang, T. Oldfield, and M. Waldman, ‘‘Ligandfit: A Novel

Method for the Shape-Directed Rapid Docking of Ligands to Protein Active Sites,’’ J. Mol. Graph. Model., 21, 289 (2003).