SMD-Based Interaction-Energy Fingerprints Can Predict Accurately

Subscriber access provided by UNIVERSITY OF SASKATCHEWAN LIBRARY

Computational Chemistry

SMD-based Interaction-energy Fingerprints Can Predict Accurately the Dissociation Rate Constants of HIV-1 Protease Inhibitors Shuheng Huang, Duo Zhang, Hu Mei, MuliadiYeremia Kevin, Sujun Qu, Xianchao Pan, and Laichun Lu J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.8b00567 • Publication Date (Web): 13 Nov 2018 Downloaded from http://pubs.acs.org on November 14, 2018

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 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 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.

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

Journal of Chemical Information and Modeling

1 2 3

SMD-based Interaction-energy Fingerprints Can Predict Accurately the Dissociation Rate Constants of HIV-1 Protease Inhibitors

4 5 6

Shuheng Huang†,‡,§, Duo Zhang‡,§, Hu Mei*,†,‡, MuliadiYeremia Kevin‡, Sujun Qu‡, Xianchao Pan‡,¶, Laichun Lu*,†,‡

7 8 9

† Key

Laboratory of Biorheological Science and Technology (Ministry of Education),

Chongqing University, Chongqing 400044, China

10

‡ College

11

¶

12

University, Luzhou, Sichuan, 646000, China

of Bioengineering, Chongqing University, Chongqing 400044, China

Department of Medicinal Chemistry, College of Pharmacy, Southwest Medical

13

1

ACS Paragon Plus Environment

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

Page 2 of 25

14

ABSTRACT

15

Recent researches have increasingly suggested that the crucial factors affecting drug

16

potencies are related not only to the thermodynamic properties but also to the kinetic

17

properties. Therefore, in silico prediction of ligand-binding kinetic properties,

18

especially the dissociation rate constant (koff), has aroused more and more attentions.

19

However, there are still a lot of challenges that need to be addressed. In this paper,

20

steered

21

decomposition was employed to predict the dissociation rate constants of 37 HIV-1

22

protease inhibitors (HIV-1 PIs). For the first time, a predictive model of the

23

dissociation rate constant was established by using the interaction-energy fingerprints

24

sampled along the ligand dissociation pathway. Based on the key fingerprints

25

extracted, it can be inferred that the dissociation rates of 37 HIV-1 PIs are basically

26

determined in the first half of the dissociation processes, and that the H-bond

27

interactions with active-site Asp25 and van der Waals interactions with flap-region

28

Ile47 and Ile50 have important influences on the dissociation processes. In general,

29

the strategy established in this paper can provide an efficient way for the prediction of

30

dissociation rate constants as well as the unbinding mechanism researches.

molecular

dynamics

(SMD)

combined

with

residue-based

energy

31 32

Keywords: Steered molecular dynamics, HIV-1 protease, Inhibitors, Dissociation rate

33

constant, Ligand-receptor interaction fingerprints

34 35

1. INTRODUCTION

36

Human immunodeficiency virus type 1 protease (HIV-1 PR) is one of the most crucial

37

enzymes in the life cycle of HIV-1, which produces mature enzymes and structural

38

proteins during virus replication process by cleaving polypeptide precursors.1 HIV-1

39

PR is a C2-symmetric, homodimeric protein composed of two chemically matching

40

monomers of 99 residues. Each monomer contains an α-helix near the C-terminus and

41

two antiparallel β-sheets.2 The C2-symmetric active site is composed of two tripeptide

42

sequences, i.e. Asp25/Asp25*-Thr26/Thr26*-Gly27/Gly27*, and is gated by two

43

extended β hairpin loops (residues 46-56), known as flap regions (Fig.1).3, 4 The main

44

function of the flap regions is to control the entrance of substrates and other small

45

molecules into the active site. The flap regions are generally believed to be partially

46

closed (semi-open form) in unbinding state (Fig 1a). However, when the substrates or

47

other small molecules enter into the active site, the flap regions will be closed to

2


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


48

prevent molecules out of the active site (Fig 1b) until protease hydrolysis prompts the

49

flap regions opening again.5, 6

50

There is a growing interest in the discovery of HIV-1 protease inhibitors (HIV-1

51

PIs) in the last three decades. HIV-1 PIs can competitively occupy the binding site of

52

HIV-1 PR substrates, and prevent the synthesis of functional proteins required by

53

virus assemble.7 US Food and Drug Administration has approved 9 HIV-1 PIs as

54

clinical medicines for AIDS treatments, namely Saquinavir (Saq), Amprenavir (Amp),

55

Atazanavir (Ataz), Darunavir (Dar), Indinavir (Ind), Lopinavir (Lop), Nelfinavir

56

(Nelf), Ritonavir (Rit), and Tipranavir (Tip).8

57

For a long time, the affinity has been regarded as a key indicator of drug potency,

58

which is often measured by thermodynamic properties, e.g., equilibrium dissociation

59

constant (KD). However, Recent researches have increasingly showed that the kinetic

60

properties, i.e. dissociation rate constant or drug-target residence time, are more

61

important for the potencies of HIV-1 PIs.9 Maschera et al.10 investigated the

62

association (kon) and dissociation rate constants (koff) of Saq for wild and mutant

63

HIV-1 PRs. The results showed that the kon of Saq were similar between the wild and

64

mutant HIV-1 PRs, while a significant difference in koff was observed. That is to say,

65

the reduced affinities for the mutant HIV-1 PRs are mainly due to the relatively fast

66

dissociation rate constants. By using biosensor technology, Shuman et al.11 measured

67

the kinetic properties of five clinical inhibitors (Amp, Ind, Nelf, Rit, and Saq) by

68

using drug-resistant variants of HIV-1 PRs. The results showed that single-amino-acid

69

substitution in the drug-resistant variants can increase the dissociation rates of the

70

HIV-1 PIs. In general, the researches revealed that the dissociation rate constant is the

71

most crucial factor determining the potencies of HIV-1 PIs.

72

Up to date, the kinetic properties are mainly determined by fluorescence

73

measurement,12 capillary electrophoresis,13 affinity chromatography,14 high pressure

74

spectroscopy,15 and surface plasmon resonance methods,16,

75

methods are still confronted with great difficulties (i.e., time-consumed, high cost, and

76

large measurement errors), which limits drug R&D in a large degree. With the

77

development of molecular simulation techniques, more and more in silico approaches

78

have been used to predict qualitatively or quantitatively the kinetic properties of small

79

molecules. Schaal et al.18 utilized comparative molecular field analysis (CoMFA) to

80

predict kinetic rate constants of 34 HIV-1 PIs. After molecular alignments and partial

81

least squares (PLS) modeling, the cross-validated determination coefficients (R2) were

3


17

etc.. However, these


82

0.72 and 0.48 for the optimal koff and kon models, respectively. By using grid-based

83

VolSurf method, Qu et al.19 established PLS models successfully for predicting kon,

84

koff, and KD of 37 cyclic and linear HIV-1 PIs. In particular, the predictive power of

85

the koff model was quite satisfied with the R2, cross-validated R2, and predictive R2 of

86

0.70, 0.70, and 0.77, respectively.

87

Many attempts have proven to be successful in the estimation of dissociation rates

88

or drug residence time by using the methods based on molecular dynamics

89

simulations.20 Li et al.21 investigated the dissociation processes of three HIV-1 PIs,

90

i.e., AHA001, XK263 and ABT538, by steered molecular dynamics (SMD) and

91

umbrella sampling techniques. The results showed that the strength of intermolecular

92

H-bonds plays a crucial role in the dissociation processes, and that the predicted

93

kinetic parameters were closely consistent with the experimental results. Buch et al.22

94

constructed a Markov state model (MSM) to predict the kinetic properties of

95

benzamidine, and the predicted kon and koff were in good agreement with the

96

experimental values. By using metadynamics method, Sun et al.23 predicted drug

97

residence time by using the optimized structures derived from holo-state proteins. In

98

comparison with the MD methods based on Poisson processes, this strategy gave a

99

comparable or even better prediction accuracies.

100

Kokh et al.24 applied τ-random acceleration molecular dynamics (τRAMD) for

101

predicting the residence time of 70 inhibitors of human N-terminal domain of heat

102

shock protein 90α. A strong correlation (R=0.81) was observed between the predicted

103

and measured residence time. Mollica et al.25 performed multiple-replica scaled-MD

104

simulations to predict the ligand residence time of Heat Shock Protein 90,

105

Glucose-Regulated Protein, and Adenosine A2A receptor, respectively. The

106

correlation coefficients between the predicted and observed residence time were 0.95,

107

0.85, and 0.95, respectively. This protocol was further applied to predict the residence

108

time of 7 glucokinase activators, and achieved satisfied results.26

109

In spite of these exploratory studies, accurate prediction of the unbinding kinetics

110

remains a challenging task. As a simple and efficient biased molecular simulation

111

method, steered molecular dynamics has been proved to be a powerful tool in the drug

112

design and molecular simulation domains,27 e.g. ligand unbinding mechanism,28

113

membrane translocation29 and prediction of drug binding affinity.30, 31 In this study,

114

the ligand-receptor interaction fingerprints extracted from SMD trajectories were, for

115

the first time, employed to predict the dissociation rate constants of 37 HIV-1 PIs with

4


Page 4 of 25

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


116

excellent prediction results. The framework presented in this paper can provide an

117

efficient way for quantitative prediction of the dissociation rate constants as well as

118

the dissociation mechanism researches.

119 120

2. METHODS

121

2.1. Molecular docking

122

A dataset of 37 HIV-1 PIs with different structural skeletons was derived from

123

Shumanet et al.32 (Table 1). There are five categories of the 37 HIV-1 PIs: non-B268

124

analogues, cyclic sulfamide analogues, P1/P1’ analogues of B268, P2/P2’ and central

125

hydroxy analogues of B268, and 7 clinical medicines. It can be seen that the kon vary

126

by more than 4 orders of magnitude and the koff by 3 orders of magnitude.

127

Based on the crystal structure of HIV-1 PR with a co-crystallized ligand

128

AHA001 (PDB: 1AJX, resolution 2.0 Å), Surflex-dock was firstly used to derive the

129

conformations of HIV-1 PR-inhibitor complexes. Before molecular docking, all the

130

37 HIV-1 PIs were charged by MMFF94 method and then optimized by Tripos force

131

field with conjugate gradient minimizer (Sybyl 8.1). The maximum iteration steps and

132

energy gradient were set to 5000 times and 0.05 kcal/mol·Å, respectively. In the

133

structural preparation of HIV-1 PR, Asp25/Asp25* were firstly protonated according

134

to previous experimental and theoretical researches.33,

135

minimization of HIV-1 PR was performed by using AMBER FF99 force field with

136

conjugate gradient algorithm. Surflex-dock employs an idealized ensemble of CH4,

137

NH, and CO probes, namely protomol, as a guide to generate putative poses of

138

ligands. In this paper, the protomol was generated based on the co-crystalized ligand

139

AHA001 with a threshold of 0.5 and bloat of 0.0 Å.35,

140

conformations were set to 5 and ring flexibility was considered in the docking

141

processes.

142

2.2. Molecular dynamics

34

36

Then, a quick global

The additional starting

143

Molecular dynamics is one of the most widely used methods in computer-aided

144

drug design and computational biology.37 MD can simulate the movement process of

145

molecular system at atomic level, and can in principle characterize the relevant

146

thermodynamics and kinetics properties.38, 39

147

Based on the optimal docking conformation for each HIV-1 PR complex, MD

148

simulation was performed to obtain an equilibrium conformation in solvent

149

environment. Each complex was firstly solvated by TIP3P water box with the

5



Page 6 of 25

150

dimensions of 90Å×90Å×90Å ， and added counter ion Cl- to neutralize the whole

151

system. Before MD simulations, 2500 steps of steepest descent followed by 2500

152

steps of conjugate gradient energy minimizations were performed. Herein, MD

153

simulations with periodic boundary conditions40, 41 were performed by NAMD 2.10

154

with CHARMM22 force field. At first, each system was gradually heated from 0 to

155

310 K within 5000 ps in the NVT ensemble, where the HIV-1 PR complex was

156

harmonically constrained by a force of 50 kcal/mol·Å2. Then, a 10-ns MD simulation

157

was performed in the NPT ensemble (310 K, 1 atm) without any constraint. The

158

integration time step was set to 2 fs, and the cutoff of van der Waals and electrostatic

159

interactions was set to 10 Å. The Particle Mesh Ewald (PME) method was used to

160

calculate the long-range electrostatic interactions,42 and the SHAKE algorithm43 was

161

used to constrain the covalent bonds with H-atoms.

162

2.3. Steered molecular dynamics

163

Steered molecular dynamics, firstly proposed by Schulten et al.44, 45 in 1998, is a

164

well-established method for non-equilibrium MD simulations. In SMD simulations, a

165

hypothetical force (or harmonic potential) was exerted to a given ligand, which can

166

accelerate dissociation process along a predefined reaction coordinate in constant

167

velocity or constant force modes.

168

Before SMD simulation, a series of computational experiments were performed on

169

a reference molecule A016 for exploring the optimal force constant k by using eight

170

gradual values (i.e. 20, 25, 30, 35, 40, 50, 60 and 70 kcal/mol·Å2). Each force rate

171

was simulated more than twice to calculate the unbinding speed. Based on the results

172

of SMD simulations of a reference molecule A016, a 6-ns SMD simulation with a

173

constant velocity of 2.5 Å/ns was performed. At this moderate pulling velocity, the

174

ligands and the nearby residues could be relaxed efficiently while reducing the

175

computation time significantly. Meanwhile, all the HIV-1 PIs can be completely

176

pulled out of the binding site of HIV-1 PR after 6-ns SMD simulations.

177

In this paper, a 6-ns SMD simulation was performed for each HIV-1 PR complex,

178

of which the starting conformation was derived from the lowest-energy conformation

179

in the last 2-ns MD trajectory. Previous studies proved that the HIV-1 PIs tended to

180

laterally slide out from the binding site in the dissociation process.46,

181

lateral direction depicted by the vector from the center of mass (COM) of Arg8 to the

182

COM of Arg8* was defined as pulling direction (Fig.2). Herein, the COM of the

183

ligand was harmonically constrained to move at a constant velocity (2.5 Å/ns) in the

6


47

Thus, the

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


184

specified direction. To avoid translation and rotation movements of the receptor, the

185

Ca atoms of Thr74/Thr74* and Leu90/Leu90* were restrained by a harmonic

186

potential with a constant force of 5 kcal/mol·Å2. The SMD trajectories were sampled

187

at a time interval of 2 ps.

188

2.4. Extraction of interaction-energy fingerprints

189

Fig.3 shows the flowchart of interaction-energy fingerprint extraction and PLS

190

modeling. In this paper, the interaction-energy fingerprints were calculated between

191

the ligand and residues within 5 Å distance away from the surface of dissociation

192

channel. A total of 48 residues were involved in the extraction of interaction-energy

193

fingerprints (Table S1). The fingerprints were then sampled from the 6-ns SMD

194

trajectory at a time interval of 300 ps, which resulted in a total of 2880 (20×48×3)

195

interaction fingerprints, i.e., 960 electrostatic interaction energies, 960 van der Waals

196

interaction energies, and 960 total interaction energies for each HIV-1 PR complex.

197

2.5 Partial least-squares modeling

198

Partial least-squares (PLS) regression48,

49

is a statistical method that combines

199

principal component analysis (PCA) and multiple regression (MLR). This method is

200

especially suited for treating the variables with strongly collinear, noisy, and

201

numerous X variables. In PLS modeling, both X and Y variables are bilinearly

202

decomposed and projected into new principal component spaces (Eq.(1-2)).

203

X=TPT+E

Eq.(1)

204

Y=UQT+F

Eq.(2)

205

Where, T and U are scores of the X and Y matrices, respectively; P is loading

206

matrix of X and Q is weight matrix of Y; E and F are residual matrices of X and Y.

207

The aim of PLS method is to find an inner linear relationship between T and U

208

matrices (Eq.(3)):

209

U=CT+G

Eq.(3)

210

where C is the regression coefficient matrix and G is the residual matrix. For details,

211

please refers to the references.48, 49

212

Before PLS modeling, the 37 HIV-1 PR complexes were randomly divided into 29

213

training samples and 8 test samples. Herein, 2880 interaction fingerprints were used

214

as independent variables and the -log(koff) was used as a target variable. Then, forward

215

selection (FS) was used for variable selection, of which the entry probability was set

216

to 0.05. Based on the variables screened from FS, PLS modeling was performed. The

7



217

determination coefficients (R2), three-fold cross-validated R2 (Q2), and R2 for the test

218

set (R2 pred) were used for model evaluations.

219 220

3. RESULTS AND DISCUSSION

221

3.1. Molecular docking and MD simulations

222

To validate the protocol of Surflex-dock, the co-crystallized AHA001 was

223

re-docked into the binding site of HIV-1 PR (PDB ID: 1AJX). The total score of the

224

optimal docking conformation of AHA001 was 13.57, and the root mean squared

225

deviation (RMSD) of the heavy atoms was 0.77 Å (Fig.4a). As shown in Fig.4b, the

226

optimal docking conformation of AHA001 forms strong H-bond interactions with

227

Asp25/Asp25*, Gly27/Gly27*, and Ile50/Ile50*, which are consistent with the

228

experimental observations. Thus, it can be concluded that Surflex-dock can reproduce

229

the native conformation of AHA001 very well. Besides, no significant correlation was

230

observed between the docking scores (Table S2) and -log(KD) of the 37 HIV-1 PIs.

231

The reasons may be that the empirical scoring function of Surflex-dock cannot

232

estimate accurately the binding free energies of the 37 HIV-1 PIs, or the Surflex-dock

233

is unable to take fully consideration of the ligand/receptor induced fit effects.

234

Fig.5 shows the total energy and RMSD evolutions in the MD simulations of 5

235

representative HIV-1 PR complexes. The results indicate that all the 5 complexes

236

reached equilibrium states after 4 ns simulations. The detailed MD results of 37

237

HIV-1 PR complexes are shown in Fig.S1-S3. For each complex, the lowest-energy

238

conformation in the last 2-ns was used for the following SMD simulations.

239

3.2. Steered molecular dynamics

240

In SMD simulations, small spring force constant k can reduce simulation errors and

241

enhance simulation accuracy.44,

50

242

optimized by using a reference molecule A016. From the plot of distance vs. time

243

(Fig.6), it can be seen that A016 tends to move at a constant speed of 2.5 Å/ns when

244

SMDk increased to 35 kcal/mol·Å2. When applying this optimal SMDk to the rest 36

245

HIV-1 PR complexes, similar results were obtained.

Herein, the force constant k (SMDk) was firstly

246

Fig.7 shows the SMD results of Saq and Ind. As shown in Fig.7a, both of the

247

inhibitors moved approximately at the constant velocity of 2.5 Å/ns by using the

248

SMDk of 35 kcal/mol·Å2. Fig.7b shows the monitored pulling forces exerted on Saq

249

and Ind in the SMD processes. It can be observed that Saq has a force peak of about

250

650 pN at 1500 ps, which is much larger than that of Ind (180 pN at 2500 ps). In

8


Page 8 of 25

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


251

comparison with Ind, Saq tends to overcome relatively higher energy barriers during

252

the dissociation process, which is in consistent with the longer residence time of Saq

253

(Table 1).

254

3.3. Interaction-energy fingerprint extraction and PLS modeling

255

Based on the SMD trajectories, a total of 2880 fingerprint descriptors were

256

extracted for each HIV-1 PR complex, which was comprised of 960 electrostatic

257

energies, 960 van der Waals energies, and 960 total energies. By using all of the 2880

258

fingerprint descriptors, a PLS model was obtained with the R2, Q2, and R2 pred were

259

0.789, 0.76, and 0.357, respectively. It is obvious that the external predictive power is

260

very poor.

261

To eliminate the irrelevant variables and enhance predictive power, PLS combined

262

with forward selection was performed. A total of 11 interaction-energy fingerprints

263

were obtained (Table 2), which were mainly related to the active-site residues (Asp25,

264

Asp30) and flap residues (Ile47, Gly49, Ile50, Phe54). From Table 2, it can be

265

observed that the predictive performances are increased significantly with the

266

increased numbers of variables. In consideration of the predictive power and model

267

explainability, Model 5 with 5 variables was chosen as the optimal PLS model, of

268

which the R2, Q2 and R2 pred were 0.749, 0.742, and 0.834, respectively. It is should

269

be noted that the predictive power of Model 5 is quite satisfactory in comparison with

270

the available models.18, 19

271

The variables included in the optimal model are T1500-Ile47*, E900-Ile50,

272

V3000-Asp25, T3000-Ile47 and V3000-Ile50*. Herein, T, E, and V represent the total,

273

electrostatic, and van der Waals interaction energies, respectively, and the subscript

274

indicates the SMD time (ps). As mentioned above, Asp25, Ile47, Ile47*, Ile50, and

275

Ile50* were determined as the key residues effecting dissociation rates, of which

276

hydrophilic Asp25 is located in the active site and hydrophobic Ile47/Ile47* and

277

Ile50/Ile50* are located in the flap region. Thus, it can be inferred that the

278

dissociation rate constants of the 37 HIV-1 PIs are mainly determined by the

279

interactions with the active-site and flap-region residues in the first half (SMD time

SMD-Based Interaction-Energy Fingerprints Can Predict Accurately

Recommend Documents