Computational Study Exploring the Interaction Mechanism of

Subscriber access provided by - Access paid by the | UCSB Libraries

Article

A Computational study exploring the interaction mechanism of benzimidazole derivatives as potent cattle bovine viral diarrhea virus inhibitors Jinghui Wang, Yinfeng Yang, Yan Li, and Yonghua Wang J. Agric. Food Chem., Just Accepted Manuscript • DOI: 10.1021/acs.jafc.6b01067 • Publication Date (Web): 29 Jun 2016 Downloaded from http://pubs.acs.org on June 30, 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 Agricultural and Food Chemistry 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 42

Journal of Agricultural and Food Chemistry

1

A Computational study exploring the interaction mechanism of

2

benzimidazole derivatives as potent cattle bovine viral diarrhea virus

3

inhibitors

4

Jinghui Wang 1, 2, †, Yinfeng Yang 1, 2, †, Yan Li 1, 2, Yonghua Wang 1,*

5

1. Key Laboratory of Xinjiang Endemic Phytomedicine Resources, Pharmacy School, Shihezi

6

University, Ministry of Education, Shihezi 832002, China; [email protected];

7

[email protected]; [email protected]

8

2. Key Laboratory of Industrial Ecology and Environmental Engineering (MOE), Faculty of

9

Chemical, Environmental and Biological Science and Technology, Dalian University of

10

Technology, Dalian, Liaoning, 116024, P R China.

11

† Co-first authors.

12

* Corresponding author: Tel.: +86-029-87092262; E-mail: [email protected]

13

Abstract: Bovine viral diarrhea virus (BVDV) infections are prevailing in cattle populations on

14

a worldwide scale. The BVDV RNA-dependent RNA polymerase (RdRp), as a promising target

15

for new anti-BVDV drug development, has attracted increasing attention. To explore the

16

interaction mechanism of 65 benzimidazole scaffold-based derivatives as BVDV inhibitors,

17

presently, a computational study was performed based on a combination of 3D-QSAR, molecular

18

docking and molecular dynamics (MD) simulations. The resultant optimum CoMFA and CoMSIA

19

models present proper reliabilities and strong predictive abilities (with Q2 = 0. 64, R2ncv = 0.93,

20

R2pred = 0.80 and Q2 = 0. 65, R2ncv = 0.98, R2pred = 0.86, respectively). In addition, there was good

21

concordance between these models, molecular docking and MD results. Moreover, the MM-PBSA

ACS Paragon Plus Environment


22

energy analysis reveals that the major driving force for ligand binding is the polar salvation

23

contribution term. Hopefully, these models and the obtained findings could offer better

24

understanding of the interaction mechanism of BVDV inhibitors as well as benefit the new

25

discovery of more potent BVDV inhibitors

26

Keywords: BVDV RdRp; 3D-QSAR; docking; molecular dynamics; MM-PBSA.

27

1 INTRODUCTION

28

Livestock (cattle, pigs and sheep) infections by members of Pestivirus genus in

29

the family flaviviridae generate heavy losses in the global economy and there is an

30

strong need to develop effective and inexpensive anti-pestiviruses agents.1 Bovine

31

viral diarrhea virus (BVDV), a prototype virus of the Pestivirus genus and a major

32

pathogen of cattle, is also responsible for the economic loss of herds all over the

33

world.2 The prevention and control of BVDV infections therefore cause of global

34

economical concern.3 In addition, BVDV, as a small enveloped, positive-stranded

35

virus, comprise a single RNA genome, which is translated into a single polyprotein.4

36

Since the RNA replication of BVDV RNA-dependent RNA polymerase (RdRp) can

37

be initiated by a de novo mechanism without a primer and the RdRp of BVDV plays a

38

critical role in viral replication,5,6 BVDV RdRp therefore serves as one of the most

39

promising drug targets for the development of new anti-BVDV agent.

40

During the past eight years, many screening efforts have been devoted to target

41

BVDV RdRp, leading to the identification and optimization of inhibitors with high

42

potency. Up to now, a series of inhibitors with different structures are collected from

43

the literature4,7-10 to inhibit the RdRp of BVDV, displaying various degrees of potency


Page 2 of 42

Page 3 of 42


44

and selectivity. However, there is no specific antiviral drug that inhibits BVDV

45

infection. Therefore, it is very important to develop antiviral lead compounds for

46

treating expensive animals in breeding programs,4 which specifically suppress the

47

virus replication. Additionally, the development of new therapeutic alternatives for

48

anti-BVDV is still difficult for researchers.

49

To facilitate the discovery of new drugs, Computer-aided drug design (CADD)

50

methods including especially the three-dimensional quantitative structure activity

51

relationship (3D-QSAR) method as the most significant applications11,12 of chemo

52

metrics, provides valuable information for future design of novel inhibitors, which

53

interact with a specific target.13-15 Currently, employing comparative molecular field

54

analysis (CoMFA)16 and comparative molecular similarity indices analysis

55

(CoMSIA)17 approaches, a total of 65 benzimidazole scaffold-based derivatives as

56

BVDV inhibitors18-20 were collected to generate 3D-QSAR models. In addition,

57

docking and molecular dynamics (MD) analysis were carried out to provide deep

58

insight into the QSAR21-22 and to explore the possible binding modes of these

59

inhibitors at the bonding site of BVDV proteins.23 The information derived from our

60

work may be helpful for the future rational design of novel BVDV inhibitors with

61

improved activity.

62

2 METHODS

63

2.1 Database Preparation and Splitting the Dataset

64

Presently, a series of 65 benzimidazole derivatives (Table S1 in Supporting

65

Information) was investigated as BVDV inhibitors.18-20 In order to use these in vitro



66

EC50 values as dependent variables in the following QSAR investigations, all the

67

inhibitory activities of these chemicals were converted into the corresponding pEC50

68

(- log EC50) values. Since the division of the training and test sets for the whole date

69

set is a crucial step for the construction of QSAR models, random selection, which is

70

comprised of rebuilding the QSAR models and assessing the resultant statistics, is a

71

popular utilized method to evaluate the robustness of models. 24,25 In the current work,

72

65 compounds were divided into training and test sets in a proportion of 4:1 for

73

generating and validating the model, respectively. To evaluate the predictive power of

74

the models, the training and test sets molecules were chosen randomly but their pEC50

75

values are uniformly distributed in terms of the value range of the whole set.

76

2.2 Molecular modeling and Quantum Chemical Calculations

77

Using the SYBYL 6.9 software package26, molecular modeling works with the

78

database alignment were carried out. To compute the partial atomic charges,

79

Gasteiger-Hückel charges were used.

80

2.3 Alignment of Compounds

81

For generating the optimal 3D-QSAR model, we use three different alignment

82

rules, i.e., Alignment-I, -II and -III. The first alignment strategy is an atom-based

83

approach, where the highest active molecule 39 (Fig. 1A) was selected as the template

84

for all compounds. Then, applying the “align database” program in SYBYL6.9, the

85

rest molecules were aligned to the common scaffold of molecule 39 to generate the

86

effective and reliable QSAR models (Fig. 1B).

87

Fig. 1


Page 4 of 42

Page 5 of 42


88

Alignment-II (Fig. 1C) is a receptor-based approach, which was generated based

89

on the derived conformation of template molecule derived from docking simulation.

90

The third alignment rule (Fig. 1D), is still a receptor-based approach. In this method,

91

all conformations of the compounds were firstly derived from dock analysis, and then

92

subjected to the process of Alignment-I, i.e., the conformation of ligand 39 was also

93

chosen as the template molecule to fit all others aligning the receptor-based

94

conformations.

95

2.4 CoMFA and CoMSIA Descriptors

96

The descriptor fields of CoMFA were obtained by employing Tripos force field.

97

For calculating the steric and electrostatic field energies, a steric probe, i.e., a sp3

98

carbon as with a van der Waals radius of 0.152 nm and an electrostatic probe, i.e., a

99

single positive +1 charge were used, respectively. For CoMSIA analysis, the CoMSIA

100

model was derived by five descriptor fields related to steric (S), electrostatic (E),

101

hydrophobic (H), H-bond donor (D), and H-bond acceptor (A). 27

102

2.5 Generation and Validation of QSAR Models

103

Using Partial least squares (PLS) analysis, the CoMFA and CoMSIA fields are

104

linearly correlated with the activity values. The leave-one out (LOO) method was

105

adopted to carry out the initial cross-validation. In order to evaluate the predictive

106

ability of obtained models, we performed the cross-validation analysis. Employing the

107

cross-validated Q2, the optimum number of components (ONC) with the lowest

108

standard error of prediction (SEP) was generated. Using the ONC, the

109

non-cross-validation analysis was performed, generating the final PLS regression



110

models for 3D-QSAR. Through several statistical parameters like the non- cross-

111

validated correlation coefficient (R2ncv), F value and the standard error of estimate

112

(SEE), the obtained non-cross-validation results were assessed. Besides, the

113

cross-validated correlation coefficient (R2pred) was evaluating the predictive capability

114

of the derived CoMFA and CoMSIA models.

115

2.6 Docking

116

For validating the compound-target associations and generating the most

117

probable binding conformation, the GOLD (version 5.1) with a widely used genetic

118

algorithm was performed. In order to predict the binding positions of ligand, we

119

optimize the GOLD Score ﬁtness function by considering van der Waals, H-bonding

120

and intra-molecular strain.28 In our studies, the X-ray crystal structure of the BVDV

121

RdRp (PDB code 1S49) was obtained from Protein Data Bank.29,30

122

2.7 Molecular Dynamics (MD)

123

Using AMBER 10.031, the MD simulations of all the molecules were carried out.

124

Under periodic boundary conditions, the simulation process employed a NPT

125

ensemble at 300 K.32 The Berendsen thermostat approach was adopted to set the

126

temperature to a constant. The values of 4.5×10-5 bar-1 for the isothermal

127

compressibility were given and the Parrinello-Rahman scheme was used to maintain

128

the pressure at 1 bar.33 To ensure the stability of the whole system, all simulations

129

were run lasted for 5 ns with a 2 fs time step.

130

2.8 MM–PBSA Calculation

131

Through all-atom molecular mechanics and Poisson-Boltzmann solvation area


Page 6 of 42

Page 7 of 42


132

(MM-PBSA) method,34 a commonly used method in determining energetics on each

133

complex

134

The following equations were used to calculate the binding free energy:

system, the binding free energy of compounds to protein was calculated.

∆Gbing = ∆EMM + ∆Gsol

135

∆GMM = ∆Einternal + ∆Eelectrostatic + ∆Evdw

136

∆Gsol = ∆GPB/GB + ∆GSA

137

(1) (2) (3)

138

Where ∆GMM is the interaction energy computed with the molecular mechanics

139

method and ∆Gsol is the solvation energy. ∆EMM denotes the molecular mechanical

140

free energy, which consists of van der Waals interaction energy (∆Evdw), the

141

electrostatic interaction energy (∆Eelectrostatic), and the internal energy of bonds (∆Eval).

142

∆Gsol encodes the solvation free energy, which is the sum of electrostatic solvation.

143

Since the negligible contribution of the entropy to the binding free energy to the same

144

protein, which is described in many literatures, 35,36presently, the entropy contribution

145

was not computed.

146

3 RESULTS

147

3.1 3D-QSAR Models

148

To assess a QSAR model’s reliability for activity prediction of unknown

149

compounds, various statistical parameters containing the Q2, R2ncv, SEE, ONC, the

150

R2pred and F-statistic values were analyzed.37 General speaking, a cross-validated Q2 >

151

0.5 is considered as the evidence of a 3D-QSAR model with satisfactory

152

predictive capability.38 Additionally, the foundation of a reliable QSAR model



153

involves the high R2ncv and F values as well as the low SEE values.39 Currently, to

154

derive the reliable models, three alignment methods, i.e., Alignment-I, -II, and -III,

155

were adopted. The statistical results of the obtained QSAR models are summarized in

156

Table 1. Clearly, among all the CoMFA and CoMSIA models, Alignment-II and -III

157

led to statistically unacceptable models with Q2 < 0.5, while Alignment-I generated

158

the models with more significant statistical values, i.e., larger Q2, R2ncv, F value, R2pred

159

and lower SEE values. Therefore, we mainly focused on Alignment-I-based models,

160

which were further used for the following discussion.

161

3.1.1 CoMFA Results

162

By analyzing the generated 3D-QSAR models, we found that through a

163

combinational use of S and E fields, the best CoMFA model was generated. The

164

optimum number of components was 4 with Q2 leave-one-out value 0.64 while the

165

PLS analysis gave a good R2ncv value of 0.93, F = 140.33, and a low SEE of 0.17.

166

These results demonstrate the stability and robustness of this model.

167

Table 1

168

Currently, to validate the constructed 3D-QSAR model, 15 BVDV inhibitors in

169

the test set were used. The predicted R2pred value for the test set from CoMFA was

170

0.80, indicating that this model has strong predictive capability. The radar correlation

171

plot of the experimental values versus the predicted activities employing the training

172

and test sets for the optimal CoMFA model was depicted in Fig. S1 (Supporting

173

Information). Clearly, a well accordance between the predicted pEC50 values and the

174

actual pEC50 values was found in a tolerable error range. These statistical parameters


Page 8 of 42

Page 9 of 42


175

suggest that during drug design process, the generated CoMFA model contour maps

176

should be taken into consideration.

177

3.1.2 CoMSIA Results

178

For CoMSIA models, using five field descriptors, i.e., S, E, H, D and A fields, all

179

31 possible combinations of descriptors were computed to build various CoMSIA

180

models. Here, the combination of S, E, H and A fields leads to the optimal CoMSIA

181

result with larger Q2, R2ncv, R2pred and lower SEE values, which was used for the

182

following discussion. Table 1 shows the statistical result of this CoMSIA model.

183

As shown in Table 1, the combination of S, E, H and A gives a good LOO

184

cross-validated Q2 of 0.65 using 10 optimum components, indicating the strong ability

185

of predicting the EC50 values of this model. In addition, a high R2ncv of 0.98 and a SEE

186

value (0.10), as well as an F-value (186.11) are obtained, presenting a good agreement

187

of the actual and the predicted pEC50 values for the non-cross-validated model. The

188

contributions fields of S, E, H and A were shown to be 14.80%, 33.10%, 29.60% and

189

22.60%, respectively.

190

Fig. 2

191

Similar to the CoMFA model, 15 compounds in the same test set were used to test

192

the reliabilities of the generated model. As shown in Table 1, an excellent value of

193

0.86 R2pred was obtained for the CoMSIA model, demonstrating its capability for

194

accurately predicting the affinity of ~86% compounds. Fig. 2 shows the scattered plot

195

of the predicted activity data versus the actual data with small variations for the

196

training and test sets. The good correlations prove the robustness and good predictive



197

ability of the optimal CoMSIA model.

198

3.2 Contour Map Analysis

199

For viewing the field effects on the target features, the 3D contour maps of the

200

optimum CoMFA and CoMSIA models were constructed and analyzed, offering good

201

insight into the new drug design.40 To aid in visualization, the most active inhibitor

202

(compound 39) along with the CoMFA and CoMSIA contour maps is depicted in Figs.

203

3 and 4, respectively.

204

Fig. 3

205

In the case of CoMFA model, Fig. 3A depicts its steric contour maps. Green and

206

yellow contours represent regions where bulky groups favor and disfavor the activity

207

of inhibitors, respectively. As observed in Fig. 3A, a large positive green isopleth

208

around the -CH3CO and a little green contour near the –NO2 show sterically favored

209

positions, meaning that chemicals with bulky groups at these positions should present

210

higher bioactivity than those without or with less bulky substituents in that location.

211

For example, a comparison of compounds 27 and 33 shows that group changed from

212

-CH3CO to -CH3CH2CO in this position might improve the activity. Moreover,

213

compound 29 ((CH2)4N-CH2CO) with relatively lower activity as compared to

214

compound 45 (CH3N(CH2CH2)2N-CH2CO) might partly due to the steric interaction.

215

Additionally, compound 23 possessing -5-NO2 group has lower potency than

216

compound 24 (-COHCH3) possibly owing to the green isopleth. Also, a big yellow

217

map near the benzimidazole, embedding the positions 9-12 and the NH-bridge at the

218

position-15, indicate that big substituents in these sites is unfavorable for the activity.


Page 10 of 42

Page 11 of 42


219

The contour maps of electrostatic descriptor are presented in Fig. 3B, where the

220

blue and red regions suggest the preference of electropositive groups and

221

electronegative substitutions, respectively. As a matter of fact, most of the excellent

222

derivatives (22, 40-47) possess an electron-donating -NH2, -NHCH3, -NCH3CH3 and

223

-OCH3 groups around the -CH3CO. As clearly observed from Fig. 3B, a red contour

224

with medium size at position-13 shows that the molecular with electron-withdraw

225

substituents at this site have higher activity. Actually, those with high

226

electron-withdrawing N or O atoms at the position-13 (1-33) are mostly inactive

227

compounds, which are partly due to the presence of red contours around here.

Fig. 4

228

229

By comparing with CoMFA model, the optimal CoMSIA model generated from

230

our work not only computes the descriptors of steric and electrostatic fields, but also

231

utilizes the hydrophobic and HB acceptor fields to correlate with the activity of

232

inhibitors. Figs. 4A and B depict the steric and electrostatic contours of CoMSIA

233

model with favoring and disfavoring activity, respectively, which are nearly

234

accordance to the corresponding CoMFA contours (Figs. 3A and B). Thus, they are

235

not discussed here and we mainly discuss the CoMSIA hydrophobic contours and the

236

H-bond acceptor contours below.

237

The hydrophobic contours of CoMSIA model are presented in Fig. 4C, where the

238

yellow and gray areas suggest the favorable and unfavorable hydrophobic feature,

239

respectively. As clearly seen from Fig. 4C, a yellow contour with optimal sized length

240

around position-1 plays a significant role for increasing the inhibitory activity of the



241

molecular, owing to the occupancy of hydrophobic groups. Moreover, trailed by the

242

yellow region, one big gray contour observed in Fig. 4C indicates that the

243

benzimidazole ring with the optimal chain length is required attached to the position-1,

244

demonstrating the favor for hydrophilic groups at this location. Additionally, the gray

245

maps with the optimal sized length presented around posotions-10 and -16 suggest

246

that the occupancy of hydrophobic substitutions will improve the activity. In fact,

247

after comparing compounds 32-52 (pEC50 > 5), a pattern for the hydrophobic groups

248

existed at these positions are found. On the contrary, the lower biological activities of

249

compounds 1-3 (pEC50 < 5) are possibly due to the favored region of the hydrophilic

250

groups.

251

In terms of the H-bond acceptor contour maps of CoMSIA model, the purple and

252

red regions displayed in Fig. 4D represent the increased and decreased activity of

253

chemicals with H-bond acceptor substituents, respectively. Clearly, a bulky favorable

254

purple contour presented near the -COCH3 of the position-16 moiety of compound 39

255

suggests the need for acceptor substitutions at this location to maintain the maximum

256

potency of compounds. By comparing the compounds 41 and 26, a conclusion can be

257

drawn that an addition of a H-bond acceptor group -OCH3 at this position would

258

generate a higher activity than -NH2, which is exemplified by the improved activity

259

from 4.79 to 6.00. In addition, two unfavorable red contours are found near the linker

260

between postions-11 and -15, implying the increased inhibitory activity of molecules

261

if they possess a H-bond donor substitute at these sites. This is in agreement with the

262

observations between molecules 41 and 25, in which compound 41with a -NH2 group


Page 12 of 42

Page 13 of 42


263

displays an improved activity than compound 25 with the NO2 in this position

264

exhibiting significantly decreased potency.

265

3.3 Flexible Docking Analysis

266

The obtained compounds were docked into BVDV RdRp using the GOLD

267

docking algorithm. Presently, all 65 chemicals were docked into the active site of

268

protein 10 times for searching the possible binding conformation. The most potent

269

compound 39 with the best pose was chosen for exploring the interaction mechanism

270

of BVDV inhibitors and used for further analysis. Fig. 5 shows the docking results of

271

chemical 39 into the active site of the receptor with the key residues located at the

272

binding pocket. In addition, Fig. S2 depicts the binding gesture of compound 39

273

inside the active site cavity. Clearly, compound 39 is suitably situated at the GTP

274

binding site, showing a good agreement with the X-ray crystal data of BVDV RdRp,

275

supporting the validation of our docking study. The topmost pose of compound 39 is

276

stabilized by hydrogen bonding and electrostatic interactions, as well as hydrophobic

277

interactions with the residues at the binding pocket.

278

As shown in Fig. 5, fitting compound 39 into the difference electron density map,

279

its orientation is supported by a basic environment containing the phosphate groups

280

(His499, Arg517, Ly525, Arg529) and a polar environment including the ribose OH

281

groups (Ser548) as well as the hydrophobic environment around the guanine base

282

(Phe549, Met552, Tyr553, Ile522), which is also almost similar to the GTP binding

283

site in BVDV RdRp.29,30 From docking studies, we observed that compound 39 and

284

GTP lie in the same pocket of BVDV RdRp, presenting the contributions of key



285

active-site residues to the ligand-receptor interaction. The binding interactions showed

286

a good binding affinity for interaction with BVDV RdRp active binding site.

Fig. 5

287

288

In addition to the above interactions with active-site residues, hydrogen

289

bonding analyses are also carried out in our docking study. It is observed from Fig. 5

290

that the core of the ligand is anchored in the active site through four

291

hydrogen-bonding interactions with the BVDV RdRp. In detail, the oxygen and

292

nitrogen atoms in the positions-15, 16 participate in the hydrogen bonding with

293

Arg529 (-O···HN, 3.14 Å) and Tyr553 (-O···HN, 2.66 Å). These hydrogen bonds are

294

perfectly filled in the center of the BVDV RdRp pocket, and play essential roles in

295

improving the inhibitors’ activity. Additionally, the nitrogen atom in the scaffold of

296

benzimidazole ring bound in the hydrogen bonding with Arg517 at a distance of 2.59

297

Å, showing a good concordance with the results of Rossmann et al. 30 In addition, the

298

-NO2 group at the terminal of compound 39 solidly form two hydrogen bonds with the

299

Arg267 (-O···HN, 2.70 and 3.03 Å).

300

Moreover, in order to further explore the binding position of the BVDV

301

derivatives, the binding modes of the lowest active compound 64 and medium affinity

302

inhibitor 29 were also investigated (Fig. S3). As demonstrated in Fig. S3A, molecule

303

64 fitted into the binding cavity which consists of His499, Ser498, Arg285, Asp268,

304

Arg267, Arg529, Lys525, Ile522, Ser548, Phe549, Met552 and Tyr553. Additionally,

305

there are two potential H-bonds formed by compound 64 with Arg529. Through a

306

comparison of the binding site of chemical 39, we find that compound 64 appears at


Page 14 of 42

Page 15 of 42


307

the almost same binding site with compound 39 and when a big group is introduced,

308

the binding mode changes clearly affect the binding affinity. Due to the strong steric

309

hindrance of this substituent, the H-bond is broken, and the nitrogen atom of the

310

substituent in turn forms two H-bonds with residue Arg529. This may be the reason

311

that compound 64 is the lowest active molecule.

312

Fig. S3B shows the binding conformation of compound 29. Obviously,

313

hydrogen-bonding and hydrophobic interactions with particular active-site residues

314

are crucial for the stabilization of the docked pose of the medium active molecule 29

315

into BVDV RdRp. In total, there are three H-bonds were found between the ligand

316

and residues Arg529 and Ser498. Compared with the binding mode of compound 39,

317

we observed that both inhibitors adopted the same binding pocket and the substituent

318

at the position-17 affects the binding energies. This is not surprising because both

319

compounds contain the same scaffold and similar substituents. However, when this

320

position of compound 39 was substituted by -OCH3, the binding energies increases

321

remarkably, owing to the hydrophobic interactions with residues Phe549, Met552,

322

Tyr553 and Ile522. Therefore, molecular optimization near this region would further

323

improve bioactivity.

324

Taken together, all the findings yield valuable insights into the molecular

325

mechanisms of the BVDV inhibitors, which will facilitate medicinal chemists in

326

identifying the scaffolds surpassing the congeneric series of chemical space used for a

327

general QSAR methodology.41

328

3.4 MD Analysis



329

To assess the reasonability of the binding mode predicted by GOLD, the initial

330

ligand geometries generated from docking model were used to perform MD protocols.

331

Herein, a 5 ns MD simulation was carried out to investigate the molecular interactions

332

involved in ligand 39 and the receptor, with main emphasis on exploring the

333

conformational changes occurred in the obtained complex.

334

Additionally, to check the dynamic stability of the simulated system, the

335

root-mean-square deviation (RMSD) of the MD trajectories for the initial complex

336

was calculated. Fig. 6A depicts the RMSDs of the system simulation for the BVDV

337

RdRp complex. Clearly, the RMSD value for the protein-ligand complex trajectory

338

reaches 1.6 to 3 Å from initial to 1.1 ns and undergoes a conformational change

339

during the MD simulation. After 2 ns, the RMSD value attains a plateau about 4Å,

340

suggesting that the model thereafter is in equilibrating state from 1.1 to 2 ns. Besides,

341

it is observed that after 2 ns of simulation, a clear metastable conformation for the

342

molecular system was found.

343

Fig. 6

344

Furthermore, a superimposition of the initial structure from docking and the

345

average structure of the complex for the last 1 ns from MD simulations was shown in

346

Fig. 6B. Obviously, the average structure and the initial docked structure of the

347

ensemble are superimposed well. In addition, after a certain time, a more consistency

348

in conformation of the docked complex of compound 39 was observed during MD

349

process, which indicates that the ligand is stable during the course of simulation.

350

Since hydrogen-bonding interactions deeply affect the specificity, adsorption and


Page 16 of 42

Page 17 of 42


351

metabolization of drugs, they play important roles in the development of drug

352

design.42 For exploring the positional conformational variations in the ligand of the

353

binding pocket, hydrogen-bonding interactions with the main native contacts in both

354

MD and docking simulations were compared (Fig. 7).43 As seen from Fig.7, a total of

355

three H-bonds appeared in both simulation systems for the protein-ligand complex. In

356

docking simulation, a backbone H-bond is observed with Tyr553, whereas in MD

357

simulation, a side chain H-bond is noted with the same residue. In addition, the NH

358

groups of Arg267 form two stable H-bonds with the ligand in docking study, showing

359

high occupancy, which are consistent with the docking results (Table S2). Moreover,

360

in MD simulation, compound 39 forms three stable hydrogen bonds with Asp268.

361

Similar to the results obtained from docking, the majority of crucial residues in the

362

binding site including Phe549, Tyr553, Ile522, Lys525, Arg517, Val521, Asp268,

363

Val269, Arg267, His499, Arg529, Glu575, Leu677, Met552 and Ser548 are observed

364

during MD simulation, showing the reliability of our docking model. Actually, there

365

are some small changes existed in the MD simulation, but the conformation of the

366

ligand in the binding pocket is still stable, indicating the rationality and reasonability

367

of our docking study. All these findings derived from our dry lab work would be

368

helpful in designing novel BVDV inhibitors. The key structural features with the

369

requirement for improving the activity of chemicals are shown in Fig. S4.

370

371 372

Fig. 7 3.5 Binding Free Energy Component Analysis With an attempt to evaluating the relationship between the biological activities of



373

inhibitors with their binding affinities, we calculated the binding free energy of

374

compound 39 with the RdRp of BVDV employing the MM-PBSA approach. The

375

MM-PBSA is composed of four terms, i.e., the van der Waals interaction energy

376

(∆Evdw), the electrostatic energy (∆Eelectrostatic), the non-polar solvation free energy

377

(∆GSA), and the polar solvation free energy (∆GPB/GB). Table 2 shows the binding free

378

energies results with the detailed contribution of each energy component.

Table 2

379

380

As shown in Table 2, the total binding free energies of the binding system in

381

drug-receptor complex simulations is -6.91 kcal·mol-1, revealing that compound 39

382

prefers binding to BVDV RdRp. Moreover, a further observation of the individual

383

free energy component suggests the important contributions of the van der Waals

384

(-24.76) and electrostatic (-42.55) interactions in the ligand-binding system.

385

Additionally, the polar interaction terms (∆Eelectrostatic + ∆GPB/GB) made a slightly

386

favorable contribution to the ligand–protein complex. With respect to the nonpolar

387

interaction contributions (∆Evdw + ∆GSA), they offer the dominant driving forces for

388

the binding, indicating that the association of the receptor-ligand complex is

389

predominantly driven by the highly favorable nonpolar interaction contribution rather

390

than by the polar ones.

391

4 DISCUSSION

392

4.1 Comparing Contours with MD Results

393

For testing the reliability of the built models, the contour maps were used for a

394

comparison with the MD results. As depicted in Figs. 3A and 4A, the large positive


Page 18 of 42

Page 19 of 42


395

yellow contour maps of CoMFA and CoMSIA models suggest that the bulky/steric

396

groups around the -CH3CO are not favored, which coincides well with the results

397

derived from MD, where the large substitutions added here will lead to a decreased

398

activity. Actually, MD studies have demonstrated that there are many residues

399

including Arg529, His499 and Tyr553 at this region, leading to a small volume of the

400

binding pocket. Therefore, too large groups introduced here may result in a steric

401

clash, which is in accordance with the steric yellow contours of the previous CoMFA

402

and CoMSIA models. In electrostatic contour maps, it is worth mentioning that

403

several electropositive residues including His499, Arg529 and Arg517 are around the

404

position-15, indicating that this region is favorable for electronegative groups, which

405

is in good accordance with the electrostatic blue contour by CoMFA and CoMSIA

406

analysis (shown in Figs. 3B and 4B). Similarly, the binding pocket also has a small

407

quantity of electronegative charge at position-9 due to the amino acid residue Asp268,

408

thus forming a favorable interaction with the negatively charged substitution of the

409

ligand, which can be deduced from the blue contours seen in above CoMFA and

410

CoMSIA analysis.

411

In addition, Fig. 7 shows that position-1 points toward the hydrophobic binding

412

pocket formed by Val521, Val269 residues. As a matter of fact, the groups at this

413

position were almost hydrophobic, which may explain the presence of yellow contour

414

maps in above CoMSIA study. Moreover, the H-bond acceptor contours of CoMSIA

415

correlate well with the H-bond interactions involved in affecting the binding and

416

activity of ligand. In detail, the oxygen atom of -NO2 at the terminal of the compound



417

39 solidly forms five bonds (Fig. 7) with the Arg267 (-O···HN, 3.38 and 2.93 Å) and

418

Asp268 (-O···HN, 3.38 and 2.93 Å), which is strongly supported by the purple

419

contour at this position as observed from the previous CoMSIA model (Fig. 4D).

420

Overall, during MD simulation, all the steric, electrostatic, hydrophobic and

421

H-bond acceptor interactions between the ligand and BVDV RdRp are in good

422

agreement with the contour maps described previously, suggesting the reasonability of

423

the QSAR models generated presently. All these results may offer useful information

424

for the discovery of new potent BVDV inhibitors in the future.

425

4.2 Comparison of Binding Mode Analysis

426

To further explore the binding mode of the BVDV inhibitors and determine

427

whether a similar mechanism of BVDV RdRp existed, docking procedure was carried

428

out with the developed BVDV inhibitors. Firstly, we directly compare the binding site

429

of compound 39 with the crystal structure of the BVDV RdRp bound to GTP. Then,

430

the findings of previous works have been compared with the current results. All

431

obtained results could provide useful insights for developing novel BVDV inhibitors

432

with improved potency.

433

4.2.1 Comparison with the crystal structure of the BVDV RdRp bound to GTP

434

By comparing the docking results of compound 39 with the crystal structure of

435

the BVDV RdRp bound to GTP, we found that both ligands are deeply located inside

436

the BVDV RdRp pocket and their binding modes are almost similar, proving the

437

reasonability of our docking process. As seen from Fig. S5, like chemical 39 in

438

docking study, the binding site for the X-ray crystal data of compound GTP was also


Page 20 of 42

Page 21 of 42


439

enclosed in the difference electron density map, which consists of the phosphate

440

groups (His499, Arg517, Ly525, and Arg529), a polar environment including the

441

ribose OH groups (Ser498) and the hydrophobic environment around the guanine base

442

(Tyr581). Additionally, for understanding the flexibilities of BVDV RdRp and

443

exploring the protein-ligand interactions, the H-bonds for ligands 39 and GTP were

444

compared as well. In Fig. S5, a H-bond formed with Arg517 was observed in X-ray

445

crystallography of the molecule with the BVDV protein, which coincides well the

446

results in the docking study. Due to the difference of the structure of the two

447

compounds and the strong steric hindrance of -NHCOCH3 group near the

448

positions-15, 16, 17 in our docking result, compound 39 has only a slight

449

conformational drift compared with the X-ray crystal data of compound GTP.

450

Moreover, four new H-bonds were formed in the docking study (Fig. 5). In detail, the

451

oxygen and nitrogen atoms of the ligand participate in the hydrogen bonding with

452

Arg529 and Tyr553 and the -NO2 group of compound 39 solidly form two hydrogen

453

bonds with the Arg267. Although the difference of the H-bonds between the

454

compound 39 and GTP of the two binding modes is small, the orientation of the

455

ligand still can be observed. Therefore, the almost similar orientation of the ligands in

456

the binding cavity for both the current work and the X-ray crystal data once again

457

demonstrates the reliability of our docking model.

458

4.2.2 Comparison with five known BVDV Inhibitors

459

Until now, there are several docking reports on BVDV RdRp receptor, making a

460

great contribution for the exploitation and optimization of small-molecule inhibitors



461

with high potency. As far as we know, five kinds of BVDV inhibitors including BPIP,

462

LZ37, TSC, AG110 and γ-carboline based on the molecule modeling were reported in

463

recent literatures and their putative binding sites were used for a comparison with our

464

docking results. Fig. S6 summarizes the structures of these molecules and the major

465

structural features impacting the activity drawn by them are summarized in Table S3,

466

where the similarities and differences of these studies are easily observed.

467

Obviously, four binding modes including Phe224, Glu291, Asn264 and Ala392

468

positions are found in Table S3. In fact, these binding modes share common

469

characteristics with the prominent hydrogen-bonding interactions and most of the

470

ligands occupy approximately the same area in the hydrophobic domains. In terms of

471

the differences in these binding modes, they are mainly the residues which participate

472

in the hydrogen-bonding and hydrophobic interactions.

473

In the case of the first binding mode, molecular docking studies have revealed

474

that Phe224 is situated in a small cavity near the top of the finger domain region of

475

the BVDV RdRp, forming possible interactions with the inhibitor BPIP.44,49 In the

476

work of Paeshuyse et al., 44 the hydrophobic region of BPIP interacting with residues

477

Ala221, Ala222, Phe224, and an H-bond between Phe224 and the imidazole ring of

478

BPIP are observed. Additionally, the γ-carboline formed a H-bond with Ala392 and an

479

aromatic-ring stacking interaction with Phe224.45 Docking of LZ37 close to Phe224

480

revealed the important roles of the hydrophobic and hydrogen-bonding interactions in

481

the binding site.46 With regard to Glu291 position, docking of AG110 into the receptor

482

showed the hydrophobic contacts of AG110 with Leu225, Leu244, Ile254, Lys255


Page 22 of 42

Page 23 of 42


483

and Tyr257.47 Moreover, a H-bond was formed by the side chain of Glu291with the

484

N7 of AG110.47 As a result, AG110 was stabilized and located between the side chains

485

and the hydrophobic interactions.

486

For Asn264 position, the docking of TSC close to Asn264 presented several

487

H-bonds interactions between the BVDV RdRp and TSC.48 In detail, a total of four

488

H-bonds were found between the TSC and the residues Asn264, Glu265 and

489

Lys266,48 indicating that the H-bonds are essential in determining the potency of the

490

TSC. Moreover, by docking TSC into BVDV RdRp, a similar interaction with the

491

Asn264 position was observed in the binding mode. The thiocarbonyl group of TSC is

492

located at the binding pocket which contains the protein residues Phe165 and Ala392

493

and two strong H-bonds are observed between residues Thr152, Thr160 and the N4 of

494

the thiosemicarbazone group.48 All these interactions are therefore required for the

495

high-affinity binding.

496

In summary, by a comparison of the four binding modes with our docking study, a

497

similar interaction mechanism was found among them, indicating that they

498

complement each other and validate our generated models. All the findings could give

499

a new strategy for the future design of novel inhibitors with high biological activities.

500

4.3 Design of New BVDV Inhibitors

501

Evaluating a QSAR model is not only by how well it fits the data, but more

502

importantly by how well it predicts data. The contour maps, docking studies and MD

503

results offer region-specific substituents on the scaffolds of compounds, which are

504

crucial for improving the activity of BVDV inhibitors. Based on these information,



505

we tried to design a set of new compounds in silico, with steric, electrostatic,

506

hydrophobic and HB acceptor substitutions in their favored regions. From molecular

507

modeling we know that appropriate substituents mainly in region 1-3 (Fig. S4) are

508

favored for the activity. Therefore, we focused the modifications on regions 1-3 in the

509

reference molecule 39 to design potent BVDV inhibitors. Interestingly, eight

510

molecules (Nos. ND01-08, Fig. S7) showed potent inhibitory activity against BVDV

511

RdRp. As shown in Fig. S7, it is obvious that all eight designed molecules show better

512

pEC50 than that of the reference molecule 39.

513

Since our current study mainly concentrates on the in sillico investigation of the

514

SAR of benzimidazole derivatives as BVDV inhibitors and the exploration of their

515

possible inhibition mechanism, the further experimental evaluation for validating the

516

inhibitory activity of these novel-designed chemicals is needed.

517

Supporting Information

518

Supporting Information Available: [Tables S1-S3: (S1) overview of the

519

Structures of BVDV inhibitors and their pEC50 values; (S2) Hydrogen bonds analysis

520

from docking and MD simulation; (S3) molecule docking results of BVDV RdRp

521

reported in other papers. Figs. S1-S7: (S1) Graph of actual and predicted biological

522

activity of CoMFA model; (S2~S5) Binding information of BVDV RdRp from

523

docking and MD simulation; (S6) Five BVDV inhibitors binding in the BVDV RdRp

524

docking model from the recent literature; (S7) Structures of novel-designed BVDV

525

compounds with their predicted inhibitory activities.] This material is available free of

526

charge via the Internet at http://pubs.acs.org.


Page 24 of 42

Page 25 of 42


527

ACKNOWLEDGEMENTS

528

Financial support from the National Natural Science Foundation of China (Nos.

529

11201049, 10801025 and 30973590) and the National High Technology Research and

530

Development Program (“863”) of China (No. 2009AA02Z205).

531

CONFLICT OF INTEREST

532

The authors confirm no conflicts of interest.

533

REFERENCE

534

1. Tonelli, M.; Vettoretti, G.; Tasso, B.; Novelli, F.; Boido, V.; Sparatore, F.; Busonera, B.;

535

Ouhtit, A.; Farci, P.; Blois, S. Acridine derivatives as anti-BVDV agents. Antivir. Res. 2011,

536

91, 133-141.

537 538 539 540

2. Gao, Y.; Zhao, X.; Sun, C.; Zang, P.; Yang, H.; Li, R.; Zhang, L. A transgenic ginseng vaccine for bovine viral diarrhea. Virol. J. 2015, 12, 1-7. 3. Xue, W.; Zhang, S.; Minocha, H.C. Characterization of a putative receptor protein for bovine viral diarrhea virus. Vet. Microbiol. 1997, 57, 105-118.

541

4. Chai, H.; Lim, D.; Chai, H.; Jung, E. Molecular modeling of small molecules as BVDV

542

RNA-dependent RNA polymerase allosteric inhibitors. B. Kor. Chem. Soc. 2013, 34, 837-850.

543

5. Baharuddin, A.; Amir Hassan, A.; Chye Sheng, G.; Bakhtiar Nasir, S.; Othman, S.; Yusof, R.;

544

Othman, R.; Abdul Rahman, N. Current approaches in antiviral drug discovery against the

545

Flaviviridae family. Curr. Pharm. Design. 2014, 20, 3428-3444.

546

6. Malet, H.; Massé, N.; Selisko, B.; Romette, J.-L.; Alvarez, K.; Guillemot, J.C.; Tolou, H.;

547

Yap, T.L.; Vasudevan, S.G.; Lescar, J. The flavivirus polymerase as a target for drug discovery.

548

Antivir. Res. 2008, 80, 23-35.



549

7. Puerstinger, G.; Paeshuyse, J.; De Clercq, E.; Neyts, J., Antiviral 2, 5-disubstituted imidazo [4,

550

5-c] pyridines: From anti-pestivirus to anti-hepatitis C virus activity. Bioorg. Med. Chem. Let.

551

2007, 17, 390-393.

552

8. Tonelli, M.; Boido, V.; La Colla, P.; Loddo, R.; Posocco, P.; Paneni,M. S.; Fermeglia, M.;

553

Pricl, S. Pharmacophore modeling, resistant mutant isolation, docking, and MM-PBSA

554

analysis: Combined experimental/computer-assisted approaches to identify new inhibitors of

555

the bovine viral diarrhea virus (BVDV). Bioorg. Med. Chem. 2010, 18, 2304-2316.

556

9. Giliberti, G.; Ibba, C.; Marongiu, E.; Loddo, R.; Tonelli, M.; Boido, V.; Laurini, E.; Posocco,

557

P.; Fermeglia, M.; Pricl, S. Synergistic experimental/computational studies on arylazoenamine

558

derivatives that target the bovine viral diarrhea virus RNA-dependent RNA polymerase.

559

Bioorg. Med. Chem. 2010, 18, 6055-6068.

560 561

10. Sako, K.; Aoyama, H.; Sato, S.; Hashimoto, Y.; Baba, M. γ-Carboline derivatives with anti-bovine viral diarrhea virus (BVDV) activity. Bioorg. Med. Chem. 2008, 16, 3780-3790.

562

11. Wang, J.; Li, Y.; Yang, Y.; Zhang, J.; Du, J.; Zhang, S.; Yang, L. Profiling the interaction

563

mechanism of indole-based derivatives targeting the HIV-1 gp120 receptor. RSC. Adv. 2015,

564

5, 78278-78298.

565

12. Wang, J.; Li, Y.; Yang, Y.; Zhang, S.; Yang, L. Profiling the Structural Determinants of

566

Heteroarylnitrile Scaffold-Based Derivatives as Falcipain-2 Inhibitors by In Silico Methods.

567

Curr. Med. Chem. 2013, 20, 2032-2042.

568 569 570

13. Guariento, S.; Bruno, O.; Fossa, P.; Cichero, E. New insights into PDE4B inhibitor selectivity: CoMFA analyses and molecular docking studies. Mol. Divers. 2016, 20, 77-92. 14. Cichero, E.; Fossa, P. Docking-based 3D-QSAR analyses of pyrazole derivatives as HIV-1


Page 26 of 42

Page 27 of 42


571

non-nucleoside reverse transcriptase inhibitors. J. Mol. Model. 2012, 18, 1573-1582.

572

15. Cichero, E.; Buffa, L.; Fossa, P. 3, 4, 5-Trisubstituted-1, 2, 4-4H-triazoles as WT and Y188L

573

mutant HIV-1 non-nucleoside reverse transcriptase inhibitors: docking-based CoMFA and

574

CoMSIA analyses. J. Mol. Model. 2011, 17, 1537-1550.

575

16. Zou, X.-J.; Lai, L.-H.; Jin, G.-Y.; Zhang, Z.-X., Synthesis, fungicidal activity, and 3D-QSAR

576

of pyridazinone-substituted 1, 3, 4-oxadiazoles and 1, 3, 4-thiadiazoles. J. Agric. Food Chem.

577

2002, 50, 3757–3760.

578

17. Huang, M.-Z.; Luo, F.-X.; Mo, H.-B.; Ren, Y.-G.; Wang, X.-G.; Ou, X.-M.; Lei, M.-X.; Liu,

579

A.-P.; Huang, L.; Xu, M.-C., Synthesis and herbicidal activity of isoindoline-1, 3-dione

580

substituted benzoxazinone derivatives containing a carboxylic ester group. J. Agric. Food

581

Chem. 2009, 57, 9585-9592.

582

18. Tonelli, M.; Paglietti, G.; Boido, V.; Sparatore, F.; Marongiu, F.; Marongiu, E.; La Colla, P.;

583

Loddo, R. Antiviral Activity of Benzimidazole Derivatives. I. Antiviral Activity of 1‐

584

Substituted‐2[(Benzotriazol‐1/2‐yl) methyl] benzimidazoles. Chem. Biodivers. 2008, 5,

585

2386-2401.

586

19. Tonelli, M.; Simone, M.; Tasso, B.; Novelli, F.; Boido, V.; Sparatore, F.; Paglietti, G.; Pricl, S.;

587

Giliberti, G.; Blois, S. Antiviral activity of benzimidazole derivatives. II. Antiviral activity of

588

2-phenylbenzimidazole derivatives. Bioorg. Med. Chem. 2010, 18, 2937-2953.

589

20. Tonelli, M.; Novelli, F.; Tasso, B.; Vazzana, I.; Sparatore, A.; Boido, V.; Sparatore, F.; La

590

Colla, P.; Sanna, G.; Giliberti, G. Antiviral activity of benzimidazole derivatives. III. Novel

591

anti-CVB-5, anti-RSV and anti-Sb-1 agents. Bioorg. Med. Chem. 2014, 22, 4893-4909.

592

21. Cichero, E.; D'Ursi, P.; Moscatelli, M.; Bruno, O.; Orro, A.; Rotolo, C.; Milanesi, L.; Fossa, P.



593

Homology modeling, docking studies and molecular dynamic simulations using graphical

594

processing unit architecture to probe the type ‐ 11 phosphodiesterase catalytic site: a

595

computational approach for the rational design of selective inhibitors. Chem. Biol. Drug Des.

596

2013, 82, 718-731.

597

22. Cichero, E.; Ligresti, A.; Allarà, M.; di Marzo, V.; Lazzati, Z.; D’Ursi, P.; Marabotti, A.;

598

Milanesi, L.; Spallarossa, A.; Ranise, A. Homology modeling in tandem with 3D-QSAR

599

analyses: a computational approach to depict the agonist binding site of the human CB2

600

receptor. Eur. J. Med. Chem., 2011, 46, 4489-4505.

601

23. Ghasemi, J.B.; Nazarshodeh, E.; Abedi, H. Molecular docking, 2D and 3D-QSAR studies of

602

new indole-based derivatives as HCV-NS5B polymerase inhibitors. J. Iran. Chem. Soc. 2015,

603

12, 1789-1799.

604

24. Ambure, P.; Roy, K. Exploring structural requirements of leads for improving activity and

605

selectivity against CDK5/p25 in Alzheimer's disease: an in silico approach. RSC. Adv. 2014, 4,

606

6702-6709.

607 608

25. Vesanto, J.; Alhoniemi, E. Clustering of the self-organizing map. IEEE Trans. Neural Network, 2000, 11, 586-600.

609

26. Good, A.C.; Hermsmeier, M.A. Measuring CAMD technique performance. 2. How “druglike”

610

are drugs? Implications of random test set selection exemplified using druglikeness

611

classification models. J. Chem. Inf. Model. 2007, 47, 110-114.

612

27. Kothandan, G.; Gadhe, C.G.; Madhavan, T.; Cho, S.J. Binding site analysis of CCR2 through

613

in silico methodologies: Docking, CoMFA, and CoMSIA. Chem. Biol. Drug. Des. 2011, 78,

614

161-174.


Page 28 of 42

Page 29 of 42


615 616 617 618

28. GOLD

User

Guide

&

Tutorials,

CCDC

Software

Ltd,

2005,

http://www.ccdc.cam.ac.uk/support/documentation/gold/4_0/gold/toc.html, 27/01/2012. 29. Westbrook, J.; Feng, Z.; Chen, L.; Yang, H.; Berman, H. M., The protein data bank and structural genomics. Nucleic Acids Res. 2003, 31, 489-491.

619

30. Choi, K.H.; Groarke, J.M.; Young, D.C.; Kuhn, R.J.; Smith, J.L.; Pevear, D.C.; Rossmann,

620

M.G., The structure of the RNA-dependent RNA polymerase from bovine viral diarrhea virus

621

establishes the role of GTP in de novo initiation. Proc. Natl. Acad. Sci. USA. 2004, 101,

622

4425-4430.

623

31. Do Kwon, Y.; Finzi, A.; Wu, X.; Dogo-Isonagie, C.; Lee, L.K.; Moore, L.R.; Schmidt, S.D.;

624

Stuckey, J.; Yang, Y.; Zhou, T. Unliganded HIV-1 gp120 core structures assume the

625

CD4-bound conformation with regulation by quaternary interactions and variable loops. Proc.

626

Natl. Acad. Sci. USA. 2012, 109, 5663-5668.

627 628

32. Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys. 1981, 52, 7182-7190.

629

33. Lin, J.-H.; Perryman, A.L.; Schames, J.R.; McCammon, J.A. Computational drug design

630

accommodating receptor flexibility: the relaxed complex scheme. J. Am. Chem. Soc. 2002,

631

124, 5632-5633.

632

34. Kollman, P.A.; Massova, I.; Reyes, C.; Kuhn, B.; Huo, S.; Chong, L.; Lee, M.; Lee, T.; Duan,

633

Y.; Wang, W. Calculating structures and free energies of complex molecules: combining

634

molecular mechanics and continuum models. Accounts. Chem. Res. 2000, 33, 889-897.

635

35. Zhou, Z.; Madrid, M.; Evanseck, J.D.; Madura, J.D. Effect of a bound non-nucleoside RT

636

inhibitor on the dynamics of wild-type and mutant HIV-1 reverse transcriptase. J. Am. Chem.



637

Page 30 of 42

Soc. 2005, 127, 17253-17260.

638

36. Zhao, Y.; Li, W.; Zeng, J.; Liu, G; Tang, Y. Insights into the interactions between HIV-1

639

integrase and human LEDGF/p75 by molecular dynamics simulation and free energy

640

calculation. Proteins 2008, 72, 635-645.

641

37. Malla, P.; Kumar, R.; Kumar, M. Validation of Formylchromane Derivatives as Protein

642

Tyrosine Phosphatase 1B Inhibitors by Pharmacophore Modeling, Atom‐Based 3D‐QSAR

643

and Docking Studies. Chem. Biol. Drug. Des. 2013, 82, 71-80.

644 645 646

38. Wang, J.; Li, F.; Li, Y.; Yang, Y.; Zhang, S.; Yang, L. Structural features of falcipain-3 inhibitors: an in silico study. Mol. BioSyst. 2013, 9, 2296-2310. 39. AbdulHameed, M. D. M.; Hamza, A.; Liu, J.; Zhan, C.-G. Combined 3D-QSAR modeling and

647

molecular

648

3-phosphoinositide-dependent protein kinase-1. J. Chem. Inf. Model. 2008, 48, 1760-1772.

649

40. Lu, X.; Zhao, L.; Xue, T.; Zhang, H. Design of novel focal adhesion kinase inhibitors using

650 651 652

docking

study

on

indolinone

derivatives

as

inhibitors

of

3D-QSAR and molecular docking. Med. Chem. Res. 2014, 23, 1976-1997. 41. Nair, P. C.; Sobhia, M. E., Fingerprint directed scaffold hopping for identification of CCR2 antagonists. J. Chem. Inf. Model. 2008, 48, 1891-1902.

653

42. Gupta, M.K.; Misra, K. Atom-based 3D-QSAR, molecular docking and molecular dynamics

654

simulation assessment of inhibitors for thyroid hormone receptor α and β. J. Mol. Model.

655

2014, 20, 1-19.

656

43. Wang, J.; Li, Y.; Yang, Y.; Du, J.; Zhang, S.; Yang, L. In silico research to assist the

657

investigation of carboxamide derivatives as potent TRPV1 antagonists. Mol. BioSyst. 2015,

658

11, 2885-2899.


Page 31 of 42


659

44. Paeshuyse, J.; Leyssen, P.; Mabery, E.; Boddeker, N.; Vrancken, R.; Froeyen, M.; Ansari, I.

660

H.; Dutartre, H.; Rozenski, J.; Gil, L. H. A novel, highly selective inhibitor of pestivirus

661

replication that targets the viral RNA-dependent RNA polymerase. J. Virol. 2006, 80,

662

149-160.

663

45. Salim, M. T.; Goto, Y.; Hamasaki, T.; Okamoto, M.; Aoyama, H.; Hashimoto, Y.; Musiu, S.;

664

Paeshuyse, J.; Neyts, J.; Froeyen, M., Highly potent and selective inhibition of bovine viral

665

diarrhea virus replication by γ-carboline derivatives. Antivir. Res. 2010, 88, 263-268.

666

46. Paeshuyse, J.; Letellier, C.; Froeyen, M.; Dutartre, H.; Vrancken, R.; Canard, B.; De Clercq,

667

E.; Gueiffier, A.; Teulade, J.-C.; Herdewijn, P. A pyrazolotriazolopyrimidinamine inhibitor of

668

bovine viral diarrhea virus replication that targets the viral RNA-dependent RNA polymerase.

669

Antivir. Res. 2009, 82, 141-147.

670

47. Paeshuyse, J.; Chezal, J.-M.; Froeyen, M.; Leyssen, P.; Dutartre, H.; Vrancken, R.; Canard,

671

B.; Letellier, C.; Li, T.; Mittendorfer, H. The imidazopyrrolopyridine analogue AG110 is a

672

novel, highly selective inhibitor of pestiviruses that targets the viral RNA-dependent RNA

673

polymerase at a hot spot for inhibition of viral replication. J. Virol. 2007, 81, 11046-11053.

674

48. Castro, E. F.; Fabian, L. E.; Caputto, M. E.; Gagey, D.; Finkielsztein, L. M.; Moltrasio, G. Y.;

675

Moglioni, A. G.; Campos, R. H.; Cavallaro, L. V. Inhibition of bovine viral diarrhea virus

676

RNA synthesis by thiosemicarbazone derived from 5, 6-dimethoxy-1-indanone. J. Virol. 2011,

677

85, 5436-5445.

678

49. Vrancken, R.; Paeshuyse, J.; Haegeman, A.; Puerstinger, G.; Froeyen, M.; Herdewijn, P.;

679

Kerkhofs, P.; Neyts, J.; Koenen, F., Imidazo [4, 5-c] pyridines inhibit the in vitro replication of

680

the classical swine fever virus and target the viral polymerase. Antivir. Res. 2008, 77, 114-119.



681 682 683 684 685 686 687 688 689 690 691


Page 32 of 42

Page 33 of 42


692 693

Fig. 1. (A) Compound 39 was used as a template for alignments. The common substructure is

694

shown in blue. Alignment-I, -II and -III of all the compounds are shown in panels (B), (C) and (D),

695

respectively. Molecules are colored in white for common C, blue for N, red for O, yellow for S,

696

cyan for H atoms, respectively.

697 698 699 700 701 702 703 704 705



706 707

Fig. 2. The ligand-based correlation plots of the predicted versus the actual pEC50 values using the

708

training (filled black dot) and the test (filled red triangles) set compounds based on the CoMSIA

709

model from random splitting.

710 711 712 713 714 715 716 717 718


Page 34 of 42

Page 35 of 42


719 720

Fig. 3 CoMFA StDev*Coeff contour maps put onto the most active molecule 39. (A) Steric

721

contour map, green and yellow contours illustrate regions where steric bulk has favorable and

722

unfavorable effects on the activity, respectively. (B) Electrostatic contour map in combination, in

723

which blue contours indicate regions where electropositive groups increase the activity, while red

724

contours indicate regions where electronegative groups increase the activity. The default values of

725

80% contribution for favored and 20% for disfavored regions were set for the visualization of the

726

CoMFA contour maps.

727

728

729

730

731

732

733

734



735 736

Fig. 4. CoMSIA StDev*Coeff contour maps put onto the most active molecule 39. (A) Steric

737

contour map, green and yellow contours illustrate regions where steric bulk has favorable and

738

unfavorable effects on the activity, respectively. (B) Electrostatic contour map in combination, in

739

which blue contours indicate regions where electropositive groups increase the activity, while red

740

contours indicate regions where electronegative groups increase the activity. (C) Hydrophobic

741

contour map. Yellow shows regions in which the introduction of hydrophobic groups is favored;

742

white indicates regions where the introduction of hydrophobic groups is disfavored; (D) HB

743

acceptor contour map. Purple shows regions where the introduction of H-accepting groups is

744

favored; red indicates regions where the introduction of H-acceptors is disfavored. The default

745

values of 80% contribution for favored and 20% for disfavored regions were set for the

746

visualization of the CoMSIA contour maps.

747 748 749 750


Page 36 of 42

Page 37 of 42


751 752

Fig. 5. Stereo view of docking interactions of compound 39 with active site amino acid residues of

753

BVDV RdRp.

754 755 756 757 758 759 760 761 762 763 764 765 766



767 768

Fig. 6. (A) View of superimposed backbone atoms of the average structure for the MD simulations

769

and the initial structure of the docking for the complex. (B) Plot of the root-mean-square deviation

770

(RMSD) of docked complex/ligand versus the MD simulation time in the MD-simulated structures.

771

Compound 39 is represented as carbon-chain in green for the initial complex and carbon-chain in

772

maroon for the average structure, respectively.

773

774

775

776

777

778 779


Page 38 of 42

Page 39 of 42


780 781

Fig. 7. Plot of the MD-simulated structures of the binding site with compound 39. H-bonds are

782

shown as dotted red lines; active site amino acid residues are represented as sticks.

783 784 785 786 787 788 789 790 791 792 793 794 795



Page 40 of 42

Table 1. Summary of CoMFA and CoMSIA results

796 PLS statistics

Alignment-I

Alignment-II

Alignment-III

CoMFA

CoMSIA

CoMFA

CoMSIA

CoMFA

CoMSIA

0.64

0.65

0.20

0.24

0.40

0.37

ONC

4

10

1

1

1

1

SEP

0.38

0.40

0.54

0.53

0.47

0.48

R2ncv

0.93

0.98

0.54

0.55

0.53

0.53

SEE

0.17

0.10

0.41

0.41

0.42

0.42

140.33

186.11

57.33

58.84

53.23

54.46

0.80

0.85

-0.09

0.21

0.581

0.59

Q2

F R2pred

Contribution (%) S

54.50

14.80

36.90

17.90

55.60

22.50

E

45.50

33.10

63.10

23.40

44.40

31.60

H

-

29.60

-

29.50

-

22.90

A

-

22.60

-

29.20

-

23.00

797

Q2 = cross-validated correlation coefficient using the LOO methods. R2ncv = non-cross-validated

798

correlation coefficient. SEE = standard error of estimate. F = R2ncv / (1-R2ncv). R2pred = predicted

799

correlation coefficient for the test set of compounds. SEP = standard error of prediction. NC =

800

optimal number of components.

801 802 803 804 805 806 807


Page 41 of 42


808

Table 2. Contributions of various energy components to the binding free energy (kcal mol-1)

809

Binding free energy and its components for the studied inhibitor with BVDV. Component ∆Evdw ∆Eelectrostatic ∆GPB/GB ∆GSA ∆Egas(EMM) ∆Gsol ∆Gbind

Receptor

Ligand

Complex

Deata

Mean

Std

Mean

Std

Mean

Std

Mean

Std

-4783.33 -40717.37 -9271.96 157.77 -45500.71 -9114.19 -54614.89

3.71 10.08 8.21 0.15 10.76 8.17 7.65

-1.08 -12.21 -27.63 1.90 -13.31 -25.72 -39.03

0.13 0.23 0.07 0.00 0.24 0.07 0.23

-4809.17 -40772.15 -9236.90 157.38 -45581.32 -9079.52 -54660.84

3.74 10.17 8.18 0.15 10.80 8.14 7.65

-24.76 -42.55 62.69 -2.30 -67.31 60.40 -6.91

0.20 0.41 0.50 0.02 0.45 0.50 0.37

810

∆Evdw, energies shown as contributions from van der Waals energy; ∆Eelectrostatic, electrostatic

811

energy; ∆GPB is the polar solvation energy with the PB model. ∆Gbind is the sum of ∆Evdw +

812

∆Eelectrostatic + ∆GPB/GB + ∆GSA. Delta is equal to mean (complex) - [mean (ligand) + mean

813

(receptor)]; Std, Standard error of mean values.

814



119x45mm (72 x 72 DPI)


Page 42 of 42

Computational Study Exploring the Interaction Mechanism of

Recommend Documents