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
ACS Paragon Plus Environment
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
Journal of Chemical Information and Modeling
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
ACS Paragon Plus Environment
17
etc.. However, these
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
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
ACS Paragon Plus Environment
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
Journal of Chemical Information and Modeling
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
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 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
ACS Paragon Plus Environment
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
Journal of Chemical Information and Modeling
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
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
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
ACS Paragon Plus Environment
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
Journal of Chemical Information and Modeling
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