Computational Insights into the Different ... - ACS Publications

Subscriber access provided by George Washington University Libraries

Article

Computational insights into the different resistance mechanism of imidacloprid versus dinotefuran in Bemisia tabaci Xiaoqing Meng, Chengchun Zhu, Yue Feng, Weihua Li, Xusheng Shao, Zhiping Xu, Jiagao Cheng, and Zhong Li J. Agric. Food Chem., Just Accepted Manuscript • DOI: 10.1021/acs.jafc.5b05181 • Publication Date (Web): 28 Jan 2016 Downloaded from http://pubs.acs.org on January 28, 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 34

Journal of Agricultural and Food Chemistry

1

Computational

2

mechanism of imidacloprid versus dinotefuran in Bemisia

3

tabaci

4

Xiaoqing Meng,† Chengchun Zhu,† Yue Feng,† Weihua Li,‡ Xusheng Shao,† Zhiping

5

Xu,† Jiagao Cheng,*,†,‡ and Zhong Li†, §

6

†

7

New Drug Design, School of Pharmacy, East China University of Science and

8

Technology, Shanghai 200237, China

9

§

10

insights

into

the

different

resistance

Shanghai Key Laboratory of Chemical Biology and ‡Shanghai Key Laboratory of

Shanghai Collaborative Innovation Center for Biomanufacturing Technology, 130

Meilong Road, Shanghai 200237, China

11 12

*Corresponding Author:

13

(J. Cheng) E-mail: [email protected]

14

Tel: +86-21-64251348

15

Fax: +86-21-64252603

16

ACS Paragon Plus Environment


17

ABSTRACT: Insecticide resistance is a critical problem for pest control and

18

management. For Bemisia tabaci, striking high metabolic resistance (generally

19

conferred by CYP6CM1) was observed for imidacloprid (IMI) and most of other

20

neonicotinoid members. However, dinotefuran (DIN) displayed very low resistance

21

factors, which indicated distinct metabolic properties. Here, molecular modeling

22

methods were applied to explore the different resistance features of IMI versus DIN

23

within the Q type of CYP6CM1. It was found that Arg225 played crucial roles in the

24

binding of IMI-CYP6CM1vQ with a cation-π interaction and two stable H-bonds,

25

however, such interactions were all absent in DIN-CYP6CM1vQ system. The stable

26

binding of IMI with CYP6CM1vQ would facilitate the following metabolic reaction,

27

while the weak binding of DIN might disable its potential metabolism, which should

28

be an important factor for their distinct resistance levels. The findings might facilitate

29

in future design of the anti-resistance neonicotinoid molecules.

30

KEYWORDS: neonicotinoids, insecticides resistance, cytochrome P450, molecular

31

dynamics simulation, B. tabaci

32


Page 2 of 34

Page 3 of 34


33

INTRODUCTION

34

Neonicotinoid insecticides that target on insect nicotinic acetylcholine receptors

35

(nAChRs)1 are by far the most successful insecticides, which sharing more than 25%

36

sales of the global insecticide market,2,3 due to their favourable safety frofile, wide

37

pest spectrum and multiple application methods.4 The first neonicotinoid

38

commercialized was imidacloprid in 1991, followed by seven other members, viz.,

39

nitenpyram, acetamiprid, thiacloprid, thiamethoxam, clothianidin, dinotefuran, and

40

sulfoxaflor.4,5 Neonicotinoids have been successfully used to control not only

41

hemipteran pest species such as aphids, plant- and leafhoppers, bugs and whiteflies,

42

etc, but also coleopteran and some lepidopteran pest species.4 Currently,

43

neonicotinoids have been used in more than 120 countries and areas,2 however, the

44

superiority of neonicotinoids was deeply challenged by the development of

45

resistance.3 Numerous cases of neonicotinoids resistance have been observed in many

46

pest species such as Bemisia tabaci (B. tabaci), Myzus persicae (M. persicae), Aphis

47

gossypii (A. gossypii), Nilaparvata lugens (N. lugens), etc.3 The resistance to

48

neonicotinoids, sometimes, even results in field failures at recommended insecticidal

49

rates, and the increased dosage may even worsens their side-effects, such as high

50

residue, toxicity on bees, etc.6

51

Generally, insecticide resistance is related with two distinct types of mechanism,

52

the increased detoxification effects,7 or the point mutation of target site.8 The

53

neonicotinoid resistance is deeply correlated with the increased detoxification roles

54

commonly conferred by cytochrome P450 monooxygenases.7 For example, the



55

over-expression of CYP6G1 conferred imidacloprid resistance in Drosophila.9 In

56

housefly, a comparative study of P450 gene expression indicated that CYP6G4 is a

57

major insecticide resistance gene involved in neonicotinoid resistance.10 In brown

58

planthopper, N. lugens, over-expression of CYP6AY1 and CYP6ER1 were found

59

contributing to the imidacloprid resistance across Asia.11,12 On the other hand, the first

60

evidence of target-site resistance was the Y151S mutation found in N. lugens,13 yet

61

such mutation has not been observed in field-caught insect populations. Another

62

R81T mutation detected in M. persicae, could draw direct effects on the resistance to

63

neonicotinoids,14,15 though the resistance in M. persicae was also correlated with

64

CYP6CY3 over-expression.16

65

Of all reported neonicotinoid resistance, most cases concerned whitefly B. tabaci,

66

one of the most destructive and invasive sucking crop pests worldwide, with two

67

major widespread biotypes, B and Q.3 Striking neonicotinoid resistance has been

68

widely reported in both biotypes of B. tabaci from several geographic regions

69

particularly against imidacloprid. For example, strains collected in Israel showed up to

70

1,000-fold resistance to thiamethoxam,17 and Q biotype populations collected in Crete

71

displayed 38- to 1,958-fold resistance to imidacloprid.18 Field populations of both B-

72

and Q-biotypes of the B. tabaci collected from southeastern China exhibited 28- to

73

1,900-fold resistance to imidacloprid and 29- to 1,200-fold resistance to

74

thiamethoxam.19 Molecular biology studies identified that the over-expression

75

cytochrome P450 monooxygenase, CYP6CM1, was strongly correlated with the high


Page 4 of 34

Page 5 of 34


76

levels of imidacloprid resistance in B. tabaci, as well as their cross-resistance

77

potential to other neonicotinoid insecticides.20

78

For the two B. tabaci biotypes, the Q biotype has become the dominant species in

79

China and other areas.21,22 In vitro expression of Q biotype version of CYP6CM1

80

protein (CYP6CM1vQ) could rapidly detoxify imidacloprid and most of other

81

neonicotinoid molecules.18 For example, CYP6CM1vQ could catalyze the

82

hydroxylation of imidacloprid to its less toxic 5-hydroxy form.23 Of all

83

commercialized neonicotinoids, however, dinotefuran showed very low resistance

84

level in B. tabaci as compared with imidacloprid and the other neonicotinoid

85

insecticides.20 Dinotefuran has a non-aromatic ring (tetrahydrofuran) structure, which

86

is different from other neonicotinoids (with pyridine or thiazole ring).4 The lack of

87

appreciable cross-resistance to dinotefuran is presumably a reflection of the substrate

88

specificity of the CYP6CM1vQ enzyme, and dinotefuran might have a distinct

89

binding mode or metabolic mechanism within cytochrome P450 monooxygenase.20

90

Considering the development of neonicotinoids resistance, there is a need to

91

explore the distinct metabolic features within different neonicotinoids, for better

92

understanding their intrinsic resistance mechanism. Computational simulations are

93

very important methods and have been successfully used in the action mechanism

94

studies of many pesticides.9,23-27 In the present study, molecular modeling, including

95

the homology modeling, docking, molecular dynamics (MD) simulations, and

96

molecular mechanics-Poisson-Boltzmann surface area (MM-PBSA) calculation, were

97

integrated to explore the distinct binding modes of CYP6CM1vQ with imidacloprid



98

(IMI) and dinotefuran (DIN). It was found that, in IMI-CYP6CM1vQ system, two

99

stable N-H···N H-bonds and a strong cation-π interaction were observed between

100

Arg225 and IMI, however in DIN-CYP6CM1vQ system, such interactions were all

101

absent. The binding of CYP6CM1vQ with IMI was stable, while with DIN was very

102

weak, which was in agreement with the distinct resistance factor (RF) values between

103

IMI and DIN. The structure and mechanistic insights obtained from the present study

104

might facilitate future design of the anti-resistance insecticides.

105

MATERIALS AND METHODS

106

Homology Modeling. The amino acid sequence of B. tabaci CYP6CM1vQ

107

(UniProt ID: B3FQ59) was retrieved from UniProtKB (http://www.uniprot.org/). The

108

human P450 enzyme CYP3A4, which accounting for oxidative metabolism of a wide

109

variety of xenobiotics,28 displays the highest similarity (32% sequence identity) with

110

CYP6CM1vQ among all available structures. Thus, the crystal structures of CYP3A4

111

(PDB entries of 3UA129 at 2.15 Å resolution and 1TQN30 at 2.05 Å resolution) are

112

selected as the templates for building the 3D model of CYP6CM1vQ. Considering the

113

templates structures do not contain residues in N-terminal membrane-binding domain,

114

the first 32 residues at the N-termini of the CYP6CM1vQ were not included in model

115

construction. The sequence alignment between CYP3A4 and CYP6CM1vQ was

116

carried out using the Align Multiple Sequences encoded in the Discovery Studio 3.5

117

software package (Accelrys Inc, 2013). The Build Homology Models Module

118

encoded in the Discovery Studio 3.5 software package (Accelrys Inc., 2013) was

119

applied to generate the CYP6CM1vQ structure. The PROCHECK31 and Profile-3D32


Page 6 of 34

Page 7 of 34


120

approaches were utilized for geometric evaluation. After above validation, a model

121

was finally chosen for further refinement by 6.0 ns MD simulations performed using

122

the AMBER12 package.33 The detailed protocol for the MD simulations was

123

described in subsequent section. The optimized model was subjected to quality

124

assessment with respect to its geometry and energy, and then was used for the

125

subsequent molecular docking.

126

Molecular Docking. The initial 3D structures of IMI and DIN were established by

127

Sybyl 7.0 (Tripos Inc) to assign the standard Tripos atom and bond types, then two

128

ligands were minimized and converted into MOL2 format. Docking study between

129

ligands and protein was performed by GOLD version 5.134 to obtain the starting

130

geometries of IMI-CYP6CM1vQ and DIN-CYP6CM1vQ complexes for simulation.

131

The CYP6CM1vQ protein structure optimized by 6.0 ns MD simulations was used as

132

the starting conformation for docking study. The standard Tripos atom and bond types

133

were assigned for the protein residues. The Fe atom of the heme group was set as the

134

center of the binding site and the binding pocket was defined enclosing the residues

135

within 10.0 Å around the center. The ChemScore scoring function parameterized for

136

heme-containing proteins was used to rank and select the docking poses. Fifty outputs

137

were generated for each docking run. All the output poses were clustered based on the

138

root mean squared deviation (RMSD) values and the poses with lower ChemScore

139

were analyzed in detail in each cluster. Finally, the reasonable binding orientation was

140

selected by considering the lower ChemScore and previously published results that

141

IMI could be hydroxylated by CYP6CM1vQ to its less toxic 5-hydroxy form23 and



142

DIN

would

undergo

143

monooxygenase.7

potential

N-demethylation

by

Page 8 of 34

cytochrome

P450

144

Molecular Dynamics Simulation. The initial models of CYP6CM1vQ complexed

145

with IMI and DIN for MD simulations were obtained from the docking results. MD

146

simulations were performed using the AMBER12 package. Geometrical optimization

147

and electrostatic potential calculation of the ligands were performed by the

148

B3LYP/6-31G(d,p) method using the Gaussian 09 program (Gaussian, Inc.,

149

Wallingford CT, 2009). The RESP35 (Restrained Electrostatic Potential) fitting

150

procedure was utilized to derive the partial atomic charges of IMI and DIN

151

compatible with the standard AMBER force field.

152

An all-atom model of CYP6CM1vQ was generated using the leap module in

153

AMBER12. The AMBER99SB all atom force field was used for the protein, and the

154

general AMBER force field was used as the parameters for ligands. For the force

155

constant parameters involving Fe, we adopted the values that were kindly provided by

156

previous work from Dr. Harris.36 The resulted models were then solvated with water

157

molecules in a truncated hexahedral periodic box, of which the TIP3P37 water model

158

was used. The distance between the box walls and the protein was set to 10.0 Å. All

159

systems for MD simulations were neutralized by adding the corresponding number of

160

counterions. Finally, the CYP6CM1vQ, IMI-CYP6CM1vQ and DIN-CYP6CM1vQ

161

systems have 53888, 56315 and 56316 atoms, respectively.


Page 9 of 34


162

Energy minimization was conducted in three steps. First, movement was allowed

163

only for the solvent and ion molecules with a harmonic constraints applied to the

164

complex. Second, the mainchain atoms of CYP6CM1vQ protein were fixed with the

165

same strength as the above and other atoms were allowed to move. Thirdly, all atoms

166

were minimized without any restraint. In each procedure, 2500 steps with the steepest

167

descent method followed by 2500 steps with the conjugated gradient method were

168

carried out. After the minimization, each system was gradually heated from 0 to 300

169

K over 20 ps under the NVT ensemble condition and equilibrated at 300 K for 20 ps.

170

Finally, 12.0 ns unrestrained MD simulations were conducted at 1 atm and 300 K

171

under the NPT ensemble condition. In the energy minimization and MD simulations,

172

particle mesh Ewald (PME)38 was employed to treat the long-range electrostatic

173

interactions and the SHAKE algorithm39 was applied to constrain the covalent bonds

174

to hydrogen atoms.

175

Binding Free Energy Calculations. Based on the equilibrated dynamic trajectory,

176

the binding free energy of each complex system was calculated using the (MM-PBSA)

177

calculation method,40-42 according to the following equation:

178

∆Gbinding=∆GMM+∆Gsolv-T∆S

(1)

179

where ∆Gbinding is the binding free energy, ∆GMM is the molecular mechanical energy,

180

∆Gsolv is the solvation energy, and T∆S is the entropy contribution. The molecular

181

mechanical energy is calculated by the following equation:

182

∆GMM=∆Gint+∆Gelec+∆Gvdw


(2)


Page 10 of 34

183

where ∆Gint, ∆Gelec, and ∆Gvdw represent internal, electrostatic, and van der Waals

184

energy in the gas phase, respectively. The solvation energy is divided into two

185

components:

186

∆Gsolv=∆Gele,sol+∆Gnonpol,sol

(3)

187

where ∆Gele,sol is the electrostatic contribution to solvation energy, and ∆Gnonpol,sol is

188

the nonpolar solvation term. Here, the polar contribution was calculated by solving

189

the Poisson-Boltzmann equation, whereas the latter is determined using,

190

∆Gnonpol,sol=γ(SASA)+b

(4)

191

where γ represents surface tension and b is constant, whereas SASA is the

192

solvent-accessible surface area (Å2).

193

In this study, 100 snapshots from the last 4.0 ns of production stage were extracted

194

for binding free energy calculations. The polar contribution term of solvation energy

195

was calculated using the PBSA program in AMBER12. The interior dielectric

196

constant was set to 1.0, and the outer dielectric constant was set to 80.0. The

197

solvent-accessible surface area was determined using the LCPO method.43 The

198

coefficient γ and b were set to 0.0072 kcal/(mol·Å2) and 0 respectively, as in the work

199

of Still and co-workers.44 Normal mode analysis45 was conducted to estimate the

200

entropic changes using the nmode program in AMBER12. Because the current

201

CYP6CM1vQ systems were relatively large (about 8000 atoms excluding water and

202

ions) and very memory demanding in calculation, thus only residues within 8.0 Å of

203

the substrate were included for the normal mode calculations. This treatment has been


Page 11 of 34


204

used in many previous studies.36,46,47 The truncated systems were minimized for up to

205

100,000 cycles to give an energy gradient of 0.0001 kcal·mol-1·Å-1 using a

206

distance-dependent dielectric constant of ε = 4r. Finally, 100 snapshots of each system

207

were selected for the entropy calculation.

208

To obtain the detailed interactions between the protein residues and the ligands, the

209

binding free energy was decomposed onto each individual residue. The MM-GBSA

210

program with the ICOSA method48 was used for this purpose. Gas-phase energies,

211

desolvation free energies and molecular mechanism were considered in energy

212

decomposition. The parameter settings were similar to the binding free energy

213

analysis.

214

RESULTS AND DISCUSSION

215

The Homology Model of B. Tabaci CYP6CM1vQ. A Blastp search revealed that

216

human CYP3A4 structures could serve as the potential template for building the

217

CYP6CM1vQ model, which shares the sequence identity about 32% with

218

CYP6CM1vQ. Considering the crystal resolution values and residue completeness,

219

the crystal structures of CYP3A4 (PDB codes: 3UA1 and 1TQN) were eventually

220

selected as the template structures to construct the CYP6CM1vQ model. The final

221

sequence alignment result between CYP6CM1vQ and CYP3A4 (3UA1 and 1TQN)

222

was depicted in Figure S1a, which was used for generating the initial 3D model of

223

CYP6CM1vQ. In PROCHECK evaluation, more residues in the most favored regions,

224

and less residues in disallowed regions implied a good stereochemical quality of the



225

model. For the homology model of CYP6CM1vQ (Figure 1a), 88.3% of the residues

226

were presented in the most favored regions, 10.1% in the additional regions, 1.4% in

227

the generously allowed regions, and only 0.2% in the disallowed regions (Figure 1b).

228

The Verify score by Proﬁle-3D for the CYP6CM1vQ structure was 192.59, which

229

was close to the top score 221.54. Moreover, most residues were reasonable with

230

positive score values, and only few residues that far away from the binding site

231

regions showed small negative profile-3D values (Figure S1b). The above results

232

indicated that the homology model of CYP6CM1vQ was reliable.

233

The modeled CYP6CM1vQ protein was then subjected to exhaustive 6.0 ns MD

234

simulations to examine the stability of the homology model, which was monitored by

235

measuring the RMSD of the protein backbone atoms with respect to the starting

236

structure. Figure 1c showed the RMSD value variation with respect to the simulation

237

time, which converged to about 2.9 Å after 2.0 ns. The results indicated that the built

238

3D models were stable during the MD simulation. The average structure obtained

239

from the last 2.0 ns of equilibration state was used for the further docking analysis.

240

Docking Results. Docking studies were performed to explore the potential

241

different binding modes of CYP6CM1vQ with IMI and DIN. For each ligand, fifty

242

docking outputs were analyzed in details. The final docking poses were obtained by

243

considering the ChemScore values and analyzing the binding modes, which was

244

depicted in Figure S2. The 5-methylene of IMI points to the heme Fe atom (Figure

245

S2a). The binding is consolidated by O-HO hydrogen bonds with residues Arg225

246

and Ser388 similar to the findings by Karunker I.,23 as well as a π-π interaction


Page 12 of 34

Page 13 of 34


247

between the pyridine ring and benzene ring of Phe130 (Figure S2a). It was in

248

agreement with the reports that IMI could be hydroxylated by CYP6CM1vQ to its less

249

toxic 5-hydroxy form.23

250

Comparing with IMI, DIN displayed similar orientation but diﬀerent binding mode,

251

in which the methyl group of DIN faces to the Fe atom, and the furanyl moiety of

252

DIN is located between the Arg225 and Ser321, with an O-HO hydrogen-bond with

253

Ser321 (Figure S2b). It is consistent with the previous reports that DIN would

254

undergo potential N-demethylation by cytochrome P450 monooxygenase.7 From the

255

docking pose and binding mode analyses, IMI might be of stronger interactions with

256

CYP6CM1vQ than DIN, which was also consistent with their ChemScore values (IMI

257

-30.24, and DIN -25.87).

258

Molecular Dynamics Simulation. MD simulations were performed on each

259

complex, not only to refine the ligands binding modes because the docking does not

260

take into account the flexibility of proteins, but to explore the dynamic behavior of the

261

enzyme and the ligands during the long-time MD simulations.

262

Root Mean Square Deviation. A stable MD trajectory is crucial for further

263

analyses. The dynamic flexibility of the MD systems was assessed by measuring the

264

RMSD values of Cα atoms throughout the whole process of simulation. The RMSD

265

values varied with respect to the simulation time were depicted in Figure 2. The

266

RMSD values of two systems displayed a large fluctuation during the first 6.0 ns and

267

reached stability after 6.0 ns. The protein atoms do not undergo significant structural



268

changes with the RMSD values of both systems converged to about 2.8 Å, a relative

269

small deviation from the minimized structure. Meanwhile, the RMSD values of ligand

270

atoms were also calculated and depicted in Figure 2. The data indicated that two

271

systems reached equilibrium states after a small fluctuation in the initial period of the

272

simulation. The magnitude of fluctuations for ligand, together with the backbone of

273

protein, leaded to a conclusion that the simulation produced stable trajectories and

274

provided a reliable basis for further analysis.

275

Position and Orientation Changes of IMI and DIN during Simulation. To

276

determine the mobility of each ligand in the MD simulation, the distance between the

277

heme Fe atom and the carbon atom C16 (Figure S2a, the 5-hydroxyl site of IMI) was

278

monitored during the MD simulation, as well as the distance with the carbon atom

279

C11 of DIN (Figure S2b, potential N-demethylation site). The average distances

280

maintained at 4.2Å and 4.0 Å (Figure S3), respectively for IMI-CYP6CM1vQ and

281

DIN-CYP6CM1vQ systems, which might facilitate their potential hydroxylation23 and

282

N-demethylation reactions catalyzed by cytochrome P450 monooxygenases.7

283

To clarify the conformational variations of the two ligands in the binding pockets

284

during MD simulations, the structures of two complexes were extracted and analyzed

285

from the snapshots at 0, 4, 8 12 ns, respectively. Figure S4 displayed the superposition

286

of IMI and DIN in the snapshots at 0, 4, 8, 12 ns in MD trajectory from the

287

IMI-CYP6CM1vQ and DIN-CYP6CM1vQ systems, respectively. During the entire

288

simulation, IMI displayed small conformation changes (Figure 3a and Figure S4a),

289

indicated that the binding of IMI is very stable, which might facilitate the


Page 14 of 34

Page 15 of 34


290

hydroxylation reaction catalyzed by CYP6CM1vQ.23 However, DIN showed

291

significant conformation variation during the whole simulation (Figure 3b and Figure

292

S4b), accompanied with large conformation changes of guanidine motif, the rotation

293

of the nitro group, and the rotation of tetrahydrofuran ring. Although DIN displayed

294

stable position of CH3 group and short interaction distance (about 4.0 Å) with the

295

heme Fe as compared with IMI (about 4.2 Å), The large conformation variation

296

indicated that DIN is difficult to find a stable binding mode with CYP6CM1vQ

297

during the whole simulation. It might draw detrimental effects on the potential

298

metabolic activity of CYP6CM1vQ with DIN, as compared with IMI.

299

Binding

Free

Energy

Analysis.

The

binding

free

energy

values

of

300

IMI-CYP6CM1vQ and DIN-CYP6CM1vQ system were calculated and analyzed

301

using the MM-PBSA method, which has high accuracy and good computational

302

efficiency in calculating the binding affinities of ligands with their targets. Table 1

303

lists the calculated energies, including the total binding energies and the individual

304

energy components. The CYP6CM1vQ displayed potent interaction with IMI (-13.12

305

kcal/mol), but weak interaction with DIN (-2.97 kcal/mol). It was consistent with their

306

distinct RF values. A recent study reported that the RF value of IMI in the Q biotype

307

strains of B. tabaci (up to 244-fold resistance) was higher than that of DIN (6-fold

308

resistance).20 A detailed binding energy decomposition analysis uncovered that the

309

sum of the electrostatic interaction energies and the van der Waals of

310

IMI-CYP6CM1vQ system (-58.75 kcal/mol and -40.28 kcal/mol, respectively) were

311

more favorable for the ligand binding than that of DIN-CYP6CM1vQ (-24.26



312

kcal/mol and -35.16 kcal/mol, respectively). The distinct electrostatic interaction

313

values contributed most to the significant binding energy difference between IMI and

314

DIN. It might be one of the accounts for the distinct resistance levels between IMI and

315

DIN in B. tabaci.

316

Key Residues Involved in Ligands Binding were Identified by Energy

317

Decomposition. In order to gain a detailed picture of protein-ligand interactions and

318

the contribution of each residue, the binding energies in two complexes were

319

decomposed on per residue located within 5 Å around the ligands by using the

320

MM-GBSA method encoded in the AMBER12 program. The residue-based energy

321

decomposition results characterized the key residues in ligand binding and identified

322

the contributions of different non-bonding interaction forces (the van der Waals and

323

the electrostatic interactions), which was helpful in the understanding of binding

324

mechanism of ligands. The residues with most favorable contributions (more than -1.0

325

kcal/mol) to the binding free energy were displayed in Figure 4a and Table S1. It was

326

found that Arg114, His128, Phe130, Arg225, Ser321, Ala322, Glu325, Ser388, Ile390,

327

and heme played key roles for the binding of IMI with CYP6CM1vQ, whereas

328

Phe130, Arg225, Phe226, Ser321, Ala322, Pro326, Ser388, and heme contributed to

329

the interaction between DIN and CYP6CM1vQ.

330

The interactions of residues Arg114, His128, Arg225, Glu325, Ser388 and Ile390

331

with IMI were evidently more favorable than with DIN, while residues Phe226,

332

Ala322 and Pro326 showed better interactions with DIN than IMI. By detailed energy

333

decomposition analysis, the interaction energy of Arg225 with IMI was up to -13.2


Page 16 of 34

Page 17 of 34


334

kcal/mol, whereas with DIN was only -5.21 kcal/mol. The large interaction energy

335

difference of Arg225 was predominantly originated from the electrostatic interaction,

336

which are -11.35 kcal/mol and -3.56 kcal/mol in IMI-CYP6CM1vQ and

337

DIN-CYP6CM1vQ systems, respectively. As a result, it also indicated that IMI and

338

DIN might be of distinct binding modes with Arg225 when complexed with

339

CYP6CM1vQ.

340

The Binding Modes of CYP6CM1vQ with IMI and DIN. The average structures

341

of the two ligand-enzyme complexes during the last 4.0 ns of the equilibrium phase

342

were extracted and examined to better elucidate the potential difference within the

343

binding modes of IMI and DIN. As depicted in Figure 4(b-c), residue Arg225 plays

344

key roles for the binding of IMI with CYP6CM1vQ. Two N-H···N H-bonds were

345

observed between the positive charged guanidium N-H groups of Arg225 and IMI.

346

Moreover, a total of 2000 snapshots during the equilibrium phase of MD simulations

347

were extracted and analyzed to see whether they were stable in the MD simulations.

348

The monitored H-bond distances were showed in figure S5 and the calculated H-bond

349

occupancy rates were listed in table 2. The two above mentioned N-H···N H-bonds

350

are very stable as revealed by the short H···N distances around 2.0 Å and 2.2 Å

351

(Figure S5a-b), respectively, as well as the high occupancy rate up to 99.80% and

352

97.15% (Table 2), respectively. Apart from the N-H···N H-bonds, a cation-π

353

interaction could be also observed between the guanidium group of Arg225 and the

354

pyridine ring of IMI during the equilibrium state of MD simulations.



355

On the other hand, a stable N-H···O H-bond was observed between the imidazole

356

N-H bond of IMI and the hydroxyl group of Ser388, with H···O distance around 2.0 Å

357

(Figure S5c) and a high H-bond occupancy rate at 98.55% (Table 2). Meanwhile, a

358

weak N-H···N H-bond was detected between the imidazole N-H of His128 and the

359

pyridine N atom of IMI, with the H···N distance about 2.5 Å (Figure S5d) and an

360

occupancy rate approximately at 50.45% (Table 2). Comparing with the docking

361

results, the previous revealed π-π stacking interaction with Phe130 was disappeared

362

during the MD simulation. One of the reasons might be the orientation changes

363

induced by the strong cation-π interaction with Arg225. The H-bonds with Arg225,

364

Ser388, His128, and the cation-π interaction with Arg225 guided and consolidated the

365

binding of IMI with CYP6CM1vQ, which could well facilitate the metabolic reaction

366

of IMI.

367

For DIN-CYP6CM1vQ system (Figure 4d-e), only a short N-H···O H-bond with

368

H···O distance at 1.9 Å was found between the CH3N-H of DIN and the carbonyl O

369

atom of Ala322 (Figure S6a), which was stable during the equilibrium state of MD

370

simulation with H-bond occupancy rates at 99.70% (Table 2). Moreover, a weak

371

N-H···O H-bond was detected between the O atom of the furanyl moiety of DIN and

372

main chain N-H of Phe226 (Figure 4d-e), with H···O distance about 2.5 Å (Figure S6b)

373

and an occupancy rate approximately at 62.05% (Table 2). No other evident H-bonds

374

could be observed, and also DIN could not form a cation-π interaction with Arg225

375

due to the lack of aromatic ring structure. As compared with IMI, DIN displayed

376

weak binding with CYP6CM1vQ, which is in accordance with the above noted


Page 18 of 34

Page 19 of 34


377

binding energy difference and the distinct RF values observed in China between IMI

378

and DIN (244-fold and 6-fold, respectively).20 Thus, the weak binding of DIN with

379

CYP6CM1vQ might lead to a weak metabolic activity of DIN, which on the other

380

hand conferred its low resistance level in Q biotype B. tabaci.

381

All data above revealed that CYP6CM1vQ displayed a quite different binding

382

potency with IMI (-13.12 kcal/mol) versus DIN (-2.97 kcal/mol), which was primarily

383

resulted from the distinct electrostatic contributions. Further energy decomposition

384

analysis revealed that Arg225 conferred the largest binding energy difference between

385

IMI-CYP6CM1vQ and DIN-CYP6CM1vQ systems. Detailed binding mode analysis

386

revealed that Arg225 played a crucial role in the binding of IMI-CYP6CM1vQ, of

387

which a cation-π interaction and two stable N-H···N H-bonds were observed between

388

Arg225 and IMI. However, no cation-π interaction and H-bonds were detected in the

389

binding of DIN with Arg225 in CYP6CM1vQ. It indicated that IMI forms stable and

390

potent binding with CYP6CM1vQ, which might well facilitate the following

391

metabolic reaction, in contrast, the binding between DIN and CYP6CM1vQ was very

392

weak, which might disable the potential metabolism of DIN. The different binding

393

potency and binding modes of IMI versus DIN within CYP6CM1vQ might be an

394

important factor contributing to their distinct resistance levels. It could also be

395

hypothesized that Arg225 contributed most to the distinct resistant features between

396

IMI and DIN in Q biotype of B. tabaci. The findings from the present study might be

397

helpful to future design of the anti-resistance neonicotinoid insecticides.



398

ASSOCIATED CONTENT

399

Supporting Information

400

Figures S1-S6 and Table S1. This information is available free of charge via the

401

Internet at http://pubs.acs.org.

402

Funding

403

This work was financial supported by the National Natural Science Foundation of

404

China (21172070, 21572059), National High Technology Research Development

405

Program of China (2011AA10A207) and the Fundamental Research Funds for the

406

Central Universities.

407

Notes

408

The authors declare no competing financial interest.

409

ABBREVIATIONS USED

410

B. tabaci, Bemisia tabaci; nAChRs, nicotinic acetylcholine receptors; IMI,

411

imidacloprid; DIN, dinotefuran; MD, molecular dynamics; RF, resistance factor;

412

MM-PBSA, molecular mechanics-Poisson-Boltzmann surface area; NCBI, National

413

Center for Biotechnology Information; RESP, Restrained Electrostatic Potential; PME,

414

particle-mesh Ewald; RMSD, root mean square deviation.

415


Page 20 of 34

Page 21 of 34


416

REFERENCES

417

(1) Jeschke, P.; Nauen, R.; Beck, M. E., Nicotinic acetylcholine receptor agonists: a milestone

418

for modern crop protection. Angew. Chem. Int. Ed. Engl. 2013, 52, 9464-85.

419

(2) Nauen, R.; Jeschke, P.; Copping, L., In focus: Neonicotinoid insecticides editorial. Pest

420

Manag. Sci. 2008, 64, 1081-1081.

421

(3) Bass, C.; Denholm, I.; Williamson, M. S.; Nauen, R., The global status of insect resistance to

422

neonicotinoid insecticides. Pestic. Biochem. Physiol. 2015, 121, 78-87.

423

(4) Jeschke, P.; Nauen, R.; Schindler, M.; Elbert, A., Overview of the status and global strategy

424

for neonicotinoids. J. Agric. Food Chem. 2010, 59, 2897-2908.

425

(5) Cutler, P.; Slater, R.; Edmunds, A. J.; Maienfisch, P.; Hall, R. G.; Earley, F. G.; Pitterna, T.;

426

Pal, S.; Paul, V. L.; Goodchild, J., Investigating the mode of action of sulfoxaflor: a fourth‐

427

generation neonicotinoid. Pest Manag. Sci. 2013, 69, 607-619.

428

(6) Fairbrother, A.; Purdy, J.; Anderson, T.; Fell, R., Risks of neonicotinoid insecticides to

429

honeybees. Environ. Toxicol. Chem. 2014, 33, 719-731.

430

(7) Casida, J. E., Neonicotinoid metabolism: compounds, substituents, pathways, enzymes,

431

organisms, and relevance. J. Agric. Food Chem. 2010, 59, 2923-2931.

432

(8) Crossthwaite, A. J.; Rendine, S.; Stenta, M.; Slater, R., Target-site resistance to

433

neonicotinoids. J. Chem. Biol. 2014, 7, 125-128.

434

(9) Jones, R. T.; Bakker, S. E.; Stone, D.; Shuttleworth, S. N.; Boundy, S.; McCart, C.; Daborn,

435

P. J.; van den Elsen, J. M., Homology modelling of Drosophila cytochrome P450 enzymes

436

associated with insecticide resistance. Pest Manag. Sci. 2010, 66, 1106-1115.



437

(10) Højland, D. H.; Vagn Jensen, K. M.; Kristensen, M., A comparative study of P450 gene

438

expression in field and laboratory Musca domestica L. strains. Pest Manag. Sci. 2014, 70,

439

1237-1242.

440

(11) Ding, Z.; Wen, Y.; Yang, B.; Zhang, Y.; Liu, S.; Liu, Z.; Han, Z., Biochemical mechanisms

441

of imidacloprid resistance in Nilaparvata lugens: over-expression of cytochrome P450 CYP6AY1.

442

Insect Biochem. Mol. Biol. 2013, 43, 1021-1027.

443

(12) Bass, C.; Carvalho, R.; Oliphant, L.; Puinean, A.; Field, L.; Nauen, R.; Williamson, M.;

444

Moores, G.; Gorman, K., Overexpression of a cytochrome P450 monooxygenase, CYP6ER1, is

445

associated with resistance to imidacloprid in the brown planthopper, Nilaparvata lugens. Insect

446

Mol. Biol. 2011, 20, 763-773.

447

(13) Liu, Z.; Williamson, M. S.; Lansdell, S. J.; Han, Z.; Denholm, I.; Millar, N. S., A nicotinic

448

acetylcholine receptor mutation (Y151S) causes reduced agonist potency to a range of

449

neonicotinoid insecticides. J. Neurochem. 2006, 99, 1273-1281.

450

(14) Hirata, K.; Kiyota, R.; Matsuura, A.; Toda, S.; Yamamoto, A.; Iwasa, T., Association

451

between the R81T mutation in the nicotinic acetylcholine receptor β1 subunit of Aphis gossypii

452

and the differential resistance to acetamiprid and imidacloprid. J. Pestic. Sci 2015, 40, 25-31.

453

(15) Bass, C.; Puinean, A. M.; Andrews, M.; Cutler, P.; Daniels, M.; Elias, J.; Paul, V. L.;

454

Crossthwaite, A. J.; Denholm, I.; Field, L. M., Mutation of a nicotinic acetylcholine receptor β

455

subunit is associated with resistance to neonicotinoid insecticides in the aphid Myzus persicae.

456

BMC Neurosci. 2011, 12, 51.

457

(16) Voudouris, C. C.; Kati, A. N.; Sadikoglou, E.; Williamson, M.; Skouras, P. J.; Dimotsiou, O.;

458

Georgiou, S.; Fenton, B.; Skavdis, G.; Margaritopoulos, J. T., Insecticide resistance status of


Page 22 of 34

Page 23 of 34


459

Myzus persicae in Greece: long‐term surveys and new diagnostics for resistance mechanisms.

460

Pest Manag. Sci. 2015.

461

(17) Kontsedalov, S.; Abu-Moch, F.; Lebedev, G.; Czosnek, H.; Horowitz, A. R.; Ghanim, M.,

462

Bemisia tabaci biotype dynamics and resistance to insecticides in Israel during the years

463

2008–2010. J. Integr. Agr. 2012, 11, 312-320.

464

(18) Roditakis, E.; Grispou, M.; Morou, E.; Kristoffersen, J. B.; Roditakis, N.; Nauen, R.; Vontas,

465

J.; Tsagkarakou, A., Current status of insecticide resistance in Q biotype Bemisia tabaci

466

populations from Crete. Pest Manag. Sci. 2009, 65, 313-322.

467

(19) Wang, Z.; Yan, H.; Yang, Y.; Wu, Y., Biotype and insecticide resistance status of the

468

whitefly Bemisia tabaci from China. Pest Manag. Sci. 2010, 66, 1360-1366.

469

(20) Qiong, R.; Xu, Y.-h.; Chen, L.; Zhang, H.-y.; Jones, C. M.; Devine, G. J.; Gorman, K.;

470

Denholm, I., Characterisation of neonicotinoid and pymetrozine resistance in strains of Bemisia

471

tabaci (Hemiptera: Aleyrodidae) from China. J. Integr. Agr. 2012, 11, 321-326.

472

(21) Simón, B.; Cenis, J.; De La Rúa, P., Distribution patterns of the Q and B biotypes of Bemisia

473

tabaci in the Mediterranean Basin based on microsatellite variation. Entomol. Exp. Appl. 2007,

474

124, 327-336.

475

(22) Pan, H.; Chu, D.; Ge, D.; Wang, S.; Wu, Q.; Xie, W.; Jiao, X.; Liu, B.; Yang, X.; Yang, N.,

476

Further spread of and domination by Bemisia tabaci (Hemiptera: Aleyrodidae) biotype Q on field

477

crops in China. J. Econ. Entomol. 2011, 104, 978-985.

478

(23) Karunker, I.; Morou, E.; Nikou, D.; Nauen, R.; Sertchook, R.; Stevenson, B. J.; Paine, M. J.;

479

Morin, S.; Vontas, J., Structural model and functional characterization of the Bemisia tabaci



480

CYP6CM1vQ, a cytochrome P450 associated with high levels of imidacloprid resistance. Insect

481

Biochem. Mol. Biol. 2009, 39, 697-706.

482

(24) Zheng, N.; Cheng, J.; Zhang, W.; Li, W.; Shao, X.; Xu, Z.; Xu, X.; Li, Z., Binding Difference

483

of Fipronil with GABAARs in Fruitfly and Zebrafish: Insights from Homology Modeling,

484

Docking, and Molecular Dynamics Simulation Studies. J. Agric. Food Chem. 2014, 62,

485

10646-10653.

486

(25) Hao, G.-F.; Wang, F.; Li, H.; Zhu, X.-L.; Yang, W.-C.; Huang, L.-S.; Wu, J.-W.; Berry, E.

487

A.; Yang, G.-F., Computational discovery of picomolar Q o site inhibitors of cytochrome bc 1

488

complex. J. Am. Chem. Soc. 2012, 134, 11168-11176.

489

(26) Lertkiatmongkol, P.; Jenwitheesuk, E.; Rongnoparut, P., Homology modeling of mosquito

490

cytochrome P450 enzymes involved in pyrethroid metabolism: insights into differences in

491

substrate selectivity. BMC Res. Notes 2011, 4, 321.

492

(27) Zhang, J.; Zhuang, S.; Tong, C.; Liu, W., Probing the molecular interaction of triazole

493

fungicides with human serum albumin by multispectroscopic techniques and molecular modeling.

494

J. Agric. Food Chem. 2013, 61, 7203-7211.

495

(28) Danielson, P., The cytochrome P450 superfamily: biochemistry, evolution and drug

496

metabolism in humans. Curr. Drug Metab. 2002, 3, 561-597.

497

(29) Sevrioukova, I. F.; Poulos, T. L., Structural and mechanistic insights into the interaction of

498

cytochrome P4503A4 with bromoergocryptine, a type I ligand. J. Biol. Chem. 2012, 287,

499

3510-3517.


Page 24 of 34

Page 25 of 34


500

(30) Yano, J. K.; Wester, M. R.; Schoch, G. A.; Griffin, K. J.; Stout, C. D.; Johnson, E. F., The

501

structure of human microsomal cytochrome P450 3A4 determined by X-ray crystallography to

502

2.05-Å resolution. J. Biol. Chem. 2004, 279, 38091-38094.

503

(31) Laskowski, R. A.; MacArthur, M. W.; Moss, D. S.; Thornton, J. M., PROCHECK: a program

504

to check the stereochemical quality of protein structures. J. Appl. Crystallogr. 1993, 26, 283-291.

505

(32) Liithy, R.; Bowie, J.; Eisenberg, D., Assessment of protein models with three-dimensional

506

profiles. Nature. 1992, 356, 83-85.

507

(33) Case, D. A.; Berryman, J. T.; Betz, R. M.; Cerutti, D. S.; Cheatham, T. E.; III; Darden, T. A.;

508

Duke, R. E.; Giese, T. J.; Gohlke, H.; Goetz, A. W.; Homeyer, N.; Izadi, S.; Janowski, P.; Kaus,

509

J.; Kovalenko, A.; Lee, T. S.; LeGrand, S.; Li, P.; Luchko, T.; Luo, R.; Madej, B.; Merz, K. M.;

510

Monard, G.; Needham, P.; Nguyen, H.; Nguyen, H. T.; Omelyan, I.; Onufriev, A.; Roe, D. R.;

511

Roitberg, A.; Salomon-Ferrer, R.; Simmerling, C. L.; Smith, W.; Swails, J.; Walker, R. C.; Wang,

512

J.; Wolf, R. M.; Wu, X.; Kollman, D. M. Y. a. P. A., AMBER12. University of California, San

513

Francisco 2012.

514

(34) Verdonk, M. L.; Cole, J. C.; Hartshorn, M. J.; Murray, C. W.; Taylor, R. D., Improved

515

protein-ligand docking using GOLD. Proteins. 2003, 52, 609-23.

516

(35) Gilson, M. K.; Sharp, K. A.; Honig, B. H., Calculating the electrostatic potential of molecules

517

in solution: method and error assessment. J. Comput. Chem. 1988, 9, 327-335.

518

(36) Harris, D. L.; Park, J. Y.; Gruenke, L.; Waskell, L., Theoretical study of the ligand–CYP2B4

519

complexes: effect of structure on binding free energies and heme spin state. Proteins. 2004, 55,

520

895-914.



521

(37) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L., Comparison

522

of simple potential functions for simulating liquid water. J. Integr. Agr. 1983, 79, 926-935.

523

(38) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G., A smooth

524

particle mesh Ewald method. J. Chem. Phys. 1995, 103, 8577-8593.

525

(39) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J., Numerical integration of the cartesian

526

equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput.

527

Phys. 1977, 23, 327-341.

528

(40) Massova, I.; Kollman, P. A., Combined molecular mechanical and continuum solvent

529

approach (MM-PBSA/GBSA) to predict ligand binding. Perspect. Drug Discov. Des. 2000, 18,

530

113-135.

531

(41) Kollman, P. A.; Massova, I.; Reyes, C.; Kuhn, B.; Huo, S.; Chong, L.; Lee, M.; Lee, T.;

532

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

533

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

534

(42) Wang, J.; Morin, P.; Wang, W.; Kollman, P. A., Use of MM-PBSA in reproducing the

535

binding free energies to HIV-1 RT of TIBO derivatives and predicting the binding mode to HIV-1

536

RT of efavirenz by docking and MM-PBSA. J. Am. Chem. Soc. 2001, 123, 5221-5230.

537

(43) Weiser, J.; Shenkin, P. S.; Still, W. C., Approximate atomic surfaces from linear

538

combinations of pairwise overlaps (LCPO). J. Comput. Chem. 1999, 20, 217-230.

539

(44) Still, W. C.; Tempczyk, A.; Hawley, R. C.; Hendrickson, T., Semianalytical treatment of

540

solvation for molecular mechanics and dynamics. J. Am. Chem. Soc. 1990, 112, 6127-6129.

541

(45) Brooks, B. R.; Janežič, D.; Karplus, M., Harmonic analysis of large systems. I. Methodology.

542

J. Comput. Chem. 1995, 16, 1522-1542.


Page 26 of 34

Page 27 of 34


543

(46) Weis, A.; Katebzadeh, K.; Söderhjelm, P.; Nilsson, I.; Ryde, U., Ligand affinities predicted

544

with the MM/PBSA method: dependence on the simulation method and the force field. J. Med.

545

Chem. 2006, 49, 6596-6606.

546

(47) Li, W.; Ode, H.; Hoshino, T.; Liu, H.; Tang, Y.; Jiang, H., Reduced catalytic activity of P450

547

2A6 mutants with coumarin: a computational investigation. J. Chem. Theory Comput. 2009, 5,

548

1411-1420.

549

(48) Gohlke, H.; Kiel, C.; Case, D. A., Insights into protein–protein binding by binding free

550

energy calculation and free energy decomposition for the Ras–Raf and Ras–RalGDS complexes.

551

J. Mol. Biol. 2003, 330, 891-913.

552



553

Figure Captions

554

Figure 1. (a) Homology model of B. tabaci CYP6CM1vQ. (b) The Ramachandran plots of

555

CYP6CM1vQ structure. (c) RMSD of the backbone of modeled CYP6CM1vQ protein.

556

Figure 2. RMSD plots for the backbone of IMI-CYP6CM1vQ and DIN-CYP6CM1vQ during MD

557

simulation.

558

Figure 3. Binding conformations IMI (a) and DIN (b) at 0, 4, 8, 12 ns of MD simulation.

559

Figure 4. (a) Ligands-residue interaction spectrum of IMI-CYP6CM1vQ and DIN-CYP6CM1vQ

560

complexes (only residues located within 5 Å of ligand were calculated); (b-c) The

561

binding mode of IMI with CYP6CM1vQ; (d-e) The binding mode of DIN with

562

CYP6CM1vQ. The 3D Ligand interaction images were created by PyMol (DeLano

563

Scientific), and the 2D ligand interaction diagrams were generated with MOE.

564


Page 28 of 34

Page 29 of 34


565

Table 1. The Binding Free Energies (kcal/mol) of IMI, DIN with CYP6CM1vQ

CYP6CM1vQ

∆ Gele

∆ Gvdw

∆ Gnonp,sol

∆ Gele,sol

-T∆ S

∆ Gbinding

IMI

-58.75

-40.28

-4.53

75.06

15.38

-13.12

DIN

-24.26

-35.16

-4.16

41.74

18.88

-2.97

566

567

Table 2. H-bonds with Occupancy Rates > 50% within the Two Systems During the Equilibrium

568

Phase of MD Simulations. Only H-bond

ligands

H-bond donor

H-bond acceptor

Occupancy rate (%)

IMI

IMI: H12

Ser388:OG

98.55

His128: HE2

IMI: N25

50.45

Arg225: HH12

IMI: N7

99.80

Arg225: HH22

IMI: N7

97.15

DIN: H20

Ala322: O

99.70

Phe226: H

DIN: O2

62.05

DIN




Page 30 of 34

Page 31 of 34


IMI DIN IMI_CYP6CM1vQ DIN_CYP6CM1vQ

4

RMSD (Angstrom)

3

2

1

0 0

4

Time (ns)

8


12



Page 32 of 34

Page 33 of 34




Table of Contents Graphic


Page 34 of 34

Computational Insights into the Different ... - ACS Publications

Recommend Documents