Characterizing Drug–Target Residence Time with Metadynamics: How

Article pubs.acs.org/jcim

Characterizing Drug−Target Residence Time with Metadynamics: How To Achieve Dissociation Rate Efficiently without Losing Accuracy against Time-Consuming Approaches Huiyong Sun,†,‡ Youyong Li,§ Mingyun Shen,∥ Dan Li,† Yu Kang,† and Tingjun Hou*,†,‡ †

College of Pharmaceutical Sciences, Zhejiang University, Hangzhou, Zhejiang 310058, P. R. China State Key Lab of CAD&CG, Zhejiang University, Hangzhou, Zhejiang 310058, P. R. China § Institute of Functional Nano and Soft Materials (FUNSOM), Soochow University, Suzhou, Jiangsu 215123, P. R. China ∥ Yangtze River Pharmaceutical Group, Taizhou, Jiangsu 225300, P. R. China ‡

S Supporting Information *

ABSTRACT: Drug−target residence time plays a vital role in drug efficacy. However, there is still no effective strategy to predict drug residence time. Here, we propose to use the optimized (or minimized) structures derived from holo-state proteins to calculate drug residence time, which could give a comparable or even better prediction accuracy compared with those calculated utilizing a large number of molecular dynamics (MD) structures based on the Poisson process. Besides, in addition to the Poisson process, one may use fewer samples for predicting residence time due to the reason that, in a large extent, the calculated drug residence time is stable and independent of the number of samples used for the prediction. With remarkably reduced computational load, the proposed strategy may be promising for large-scale drug residence time prediction, such as post-processing in virtual screening (VS) and lead compound optimization.

■

INTRODUCTION The past decade has witnessed a great success in drug discovery, where a very important conception has been introduced by Copeland et al., namely, drug−target residence time,1 which is thought to play critical role in drug efficacy in open systems, such as human bodies.2,3 Drug−target residence time emphasizes how long a drug stays in its target rather than how tightly the drug binds to it. Here, 1/koff is often used to characterize drug−target residence time, where koff is the dissociation rate constant of a reaction and usually used as the aim for various predictions. Accurate prediction of drug residence time is notoriously difficult since the typical time scale for drug dissociation ranges from minutes to days,1,4,5 which is unrealistic for traditional molecular dynamics (MD) simulations even conducted with special purpose computers, such as Anton, which can extend the MD time scale to millisecond (ms).6−8 Neither is it feasible for the classical enhanced sampling approaches, such as steered molecular dynamics (SMD) in conjunction with Jarzynski equality,9 umbrella sampling (US), 10 and adaptive biasing force (ABF),11,12 which mostly work with time-independent frameworks and can only estimate drug−target residence time from convergent free energy profiles,13−15 thereafter making the prediction of drug residence time very time consuming as well. Nevertheless, up to date, numerous works have been done to accurately predict drug−target residence time (or dissociation rate). For instance, initiating from different microstates, Buch et © 2017 American Chemical Society

al. conducted 495 conventional MD simulations (with 100 ns for each) to build the Markov state model (MSM) for the benzamidine−trypsin system to predict the transition-state kinetics (i.e., the binding/unbinding process) of the ligand, and they found several key metastable states for drug association and dissociation, making the predicted kinetics parameters closely consistent with the experimental kon and koff values.16 Barril and co-workers proposed to use water-shielded hydrogen bonds to predict the binding kinetics which determines the association/dissociation processes of small molecular drugs.17 Moreover, as mentioned above, the association/dissociation rate constant for a ligand to its target can be easily obtained from convergent free energy surfaces as well. For example, through two-dimensional US simulations, Tiwary and colleagues constructed the convergent free energy landscape to predict the influence of water molecules to the association and dissociation processes of the ligand−receptor interaction.14 Moreover, using an advanced molecular docking scheme, Bai et al. investigated the whole conformational space of drug (−)-Huperzine A to its target acetylcholinesterase (AChE), and they could estimate both the binding affinity and the kinetics-associated parameters (i.e., kon and koff values) simultaneously.18 Besides, as the most rapid free energy landscape characterization methodology, various versions of Received: February 9, 2017 Published: July 27, 2017 1895

DOI: 10.1021/acs.jcim.7b00075 J. Chem. Inf. Model. 2017, 57, 1895−1906

Article

Journal of Chemical Information and Modeling

Figure 1. Structures (panels A1−H1), typical free energy landscapes (panels A2−H2, based on the minimized structures), and drug dissociation trajectories (panels A3−H3) of the investigated systems. The favorable dissociation pathway in each system is plotted in the black dotted line in the corresponding free energy landscape (panels A2−H2), and the detected dissociation time is shown in the red dotted line in panels A3−H3.

metadynamics19−21 have been employed for the predictions of drug kinetics-associated parameters.22,23 Gaspari and co-workers used a metadynamics-based convergent free energy surface (FES) to determine the association and dissociation rate

constants and the binding pathways of two drugs (parasubstituted n-alkyl and hydroxy-ethylene-benzenesulfonamides) to their target hCAII, which could give consistent kinetics parameters compared with their surface plasmon resonance 1896

DOI: 10.1021/acs.jcim.7b00075 J. Chem. Inf. Model. 2017, 57, 1895−1906

Article

Journal of Chemical Information and Modeling (SPR) measurements.24 Recently, Cavalli’s group used a newly developed approach, called scaled-MD,25 to determine the kinetic parameters of a set of ligands targeting GK1 for the treatment of type 2 diabetes mellitus, which showed good agreement between the experimental data as well.26 Furthermore, Bortolato et al. employed the maximum biasing potential derived from metadynamics to approximately represent the drug residence time.27 Although the approach can efficiently discriminate the long-residence time ligands from the shortresidence time ligands, the calculated absolute residence times are mostly several orders of magnitude lower than the experimental data, which we discuss in the following part of the paper. Besides, our group also used a convergent free energy landscape derived from funnel-based well-tempered metadynamics (FWTM) to investigate the crizotinib resistance mechanism in the G2032R mutant of ROS1, and we found that the attenuated drug residence time (1/koff) is the main reason for the serious drug resistance phenomenon in the G2032R mutated ROS1.28 Although the pioneering works enabled us to successfully and accurately predict the kinetics-associated parameters (such as ligands binding/unbinding to their targets), it is indeed very time consuming to determine the kinetics-associated parameters for a reaction (which usually needs tens of microseconds for conventional MD simulations and several microseconds for classical enhanced sampling methodologies) even with the most rapid FES characterization method, such as metadynamics29,30 (which typically needs hundreds of nanoseconds to converge the free energy landscape for complicated systems, i.e., drug− target interaction31−34), thereafter hindering the methodologies being implemented to large-scale tasks, for instance, postprocessing in virtual screening. Fortunately, Tiwary and Parrinello introduced a simple yet very powerful metadynamics-based method to predict the transition rates between different metastable states recently,35 which can estimate the drug residence time based on the first dissociation event rather than characterizing whole convergent FESs.36 Although Salvalaglio et al. showed that accurate prediction of transition time needs a large number of independent simulations that distribute according to Poisson distribution,37 it is still unclear whether the simplest case, namely, using the single optimized structure derived from holo-state proteins, is capable of generating a comparable prediction. If so, the computational cost of the drug residence time prediction will be significantly reduced to facilitate the application in large-scale predictions. Therefore, herein, we systematically investigated the performance of the Tiwary algorithm using both the optimized structures and independent MD structures. The Bortolato algorithm is also added to give a comparison. The proposed strategy may furnish useful information for large-scale prediction of drug residence time.

drug residence time calculations. The missing residues such as the A-loop in EGFR were modeled with the loop module in Sybyl-X1.2. The protonation state of the proteins, such as histidines, were determined by PROPKA version 3.1.44 To obtain the simulation parameters for the small molecules, the Hartree−Fock method based on the 6-31G basis set was first used for the structure optimization, followed by the electrostatic potential calculation. The atomic charges were fitted by the restrained electrostatic potential technique (RESP)45 embedded in the Amber14 simulation package. The structural optimization was accomplished by Gaussian09.46 The ff99SB-ILDN force field47 and general Amber force field (gaff)48 were used for the proteins and drugs, respectively. Na+ and Cl− were used to neutralize the unbalanced charges of the investigated systems. A TIP3P49 water box was finally added to the drug−target complexes with the water molecules extended 10 Å out of the boundary of the solute. Molecular Mechanics Minimization. To remove the unfavorable contacts of the repaired side chains of the missing residues, molecular mechanism minimization was first used to optimize the structure of the ligand−receptor complexes. Here, the commonly used four steps of minimization strategy was employed for the investigated systems:50−52 (1) All the hydrogen atoms were relaxed for 5000 steps with all the heavy atoms restrained. (2) The heavy atoms in solvent (including the O atom in water molecules and ions) were released to be free (5000 steps minimization). (3) Only the backbone atoms in the protein were restrained with other atoms set free (5000 steps minimization). (4) The system was minimized for 20,000 steps without any restraint on the atoms. In the optimization process, the cutoff value of the short-range interactions (including van der Waals and short-range electrostatic interactions) was set to 10 Å, and the long-range interactions (long-range electrostatic interaction) were handled by the particle mesh Ewald algorithm (PME).53 Conventional MD Simulation. In the process of MD simulation, to prevent the too fast vibration of the hydrogen atoms in the system, all covalent bonds involving hydrogen atoms were restrained with a SHAKE algorithm,54 and the time step was set to 2 fs. On the basis of the optimized structures, the systems were at first gradually heated from 0 to 310 K within 1 ns simulation time in an NVT ensemble (with the backbone in proteins restrained with 5 kcal/mol Å2). Then, the systems were relaxed with the restraint gradually decreased from 5 to 0 kcal/mol Å2 within another 1 ns in an NPT ensemble (under the targeted temperature of 310 K and pressure of 1 atm). Here, the temperature was maintained by Langevin dynamics, and the pressure was controlled by the Poisson Piston algorithm.55 Finally, a 50 ns production run was performed for each system with the collection interval set to 20 ps (namely, 2500 frames for each system). All the simulations were performed with NAMD, version 2.9.56 Metadynamics Simulation. Metadynamics29,30 has been widely used in the scope of drug−target interaction, where the binding/unbinding process of the drugs can be observed repeatedly.28,31,57,58 Through periodically adding historydependent biasing potential on the selected reaction coordinates, the free energy surface can be flattened to accelerate the occurrence of the rare events, such as the dissociation of drugs from their targets, which may occur from minutes to days and is notoriously difficult to observe in conventional MD simulations. Equation 1 shows the most widely used metadynamics formalism, namely, well-temped

■

METHODS Simulation Systems Setup. Six systems derived from Dahl’s work,4 including HIV protease complexed with lopinavir (1MUI38), atazanavir (2AQU39), and indinavir (3WSJ40); renin complexed with aliskiren (2V0Z41); EGFR complexed with lapatinib (3BBT42); and neuraminidase complexed with oseltamivir (4HZX43), were employed for the drug residence time calculations. Specially, two conformations of lopinavir (Figure 1A and B) and atazanavir (Figure 1C and D) were found in the crystallized structures of HIV protease (1MUI and 2AQU). Thereafter, the two conformations were both used for 1897

DOI: 10.1021/acs.jcim.7b00075 J. Chem. Inf. Model. 2017, 57, 1895−1906

Article

Journal of Chemical Information and Modeling metadynamics,59 which was also used here for the investigation of the dissociation events.

initial hill height (ω) was set to 1 kcal/mol ps, and the deposition time (τ) was set to 2 ps to give enough dissociation time for the release of drugs without adding biasing potential on the dissociation boundary.35 The bias-factor (ΔT) was set to 3100 corresponding to 10 times the simulation temperature (310 K). Residence Time Calculation. As mentioned above, it is notoriously difficult to calculate drug residence time from conventional MD simulation and classical enhanced sampling methodologies (such as umbrella sampling). Fortunately, several metadynamics-based algorithms have been proposed to calculate drug residence time, such as the Tiwary algorithm based on transition state theory35 and Bortolato algorithm based on maximum biasing potential.27 Thus, here we used the two methods for residence time calculation. In the Tiwary algorithm, they deduced the residence time (τ) through transition state theory, where the residence time can be represented by the product of the metadynamics-associated transition time (τM) and the acceleration factor (α) as shown in eq 2. The calculation of the acceleration factor is shown in eq 3.

⎛ (s − s(t ′))2 ⎞ ωe−V (s , t ′)/ΔT exp⎜ − ⎟ 2σ 2 ⎠ ⎝ t ′= τ ,2τ ,··· t

V (s ,t ) =

(t ′ < t )

∑

(1)

In eq 1, V(s,t) is the history-dependent Gaussian potential added on the RC at position s and cumulative time t, and ωe−V(s, t′)/ΔT denotes the hill height of the biasing potential, which will be exponentially scaled down with the increase in the concurrent biasing potential V(s,t′) at each time interval τ. Here, ΔT is called the bias-factor and is used for controlling the decreasing rate of the biasing potential. After sufficient simulation time, the FES of the selected RCs can be flattened and sampled with the same probability. Nevertheless, since we only need to calculate the acceleration factor (α) here, which is associated with the dissociation time of the drugs and is detailed below (eq 3), it is not necessary to converge the whole free energy landscape. It is well known that the slowest degree of freedom usually fulfills the bottleneck of a reaction, such as the dissociation of the drugs from their active pockets, which is a typical rare event and can be hardly observed in traditional MD simulations. To solve this problem, therefore, the center of mass (CoM) distance between the drug and the corresponding active pocket (heavy atoms in residues within 5 Å of the drug) is employed as the main RC (RC1) for the dissociation event simulation. Recently, Salvalaglio et al. have shown that it is vital to use multidimensional reaction coordinates to accurately characterize the residence time of the biased molecule.37 Thereafter, we also employed the conformational change of the active pocket (atoms same as RC1, characterized by the RMSD change of the active pocket) as the second RC for the metadynamics simulation (RC2), which considers the induced-fit phenomenon of the drug binding and unbinding processes (Figure 1A2− 1G2) and is also known as one of the slowest degree of freedoms in drug association/dissociation events.60 To sufficiently assess the performance of metadynamics upon drug residence time calculation, the minimized (optimized) structure and 48 structures derived from the corresponding last 48 ns MD trajectory (that is, after 2 ns equilibrium, the structures were collected each 1 ns from the last 48 ns MD trajectory, where a total of 50 ns of MD simulation is performed for each system) for each drug−target complex were used for the dissociation process simulation. Each structure was at first overlapped upon the minimized structure and immersed into a TIP3P water box with the box boundary extended 12 Å out of the solute, which is sufficient for drug dissociation since the maximum radius of the drug (atazanavir, radius of ∼8.5 Å) is less than 10 Å. Then, each system derived from the corresponding MD trajectory was minimized for 20,000 steps with the heavy atoms in the protein and ligand restrained with 5 kcal/mol Å2. Afterward, all the systems were submitted to a short time NVT simulation (1 ns) to equilibrate the systems (with 2 kcal/mol Å2 restraint on the heavy atoms of the drug− target complexes). Finally, the last structures restarted from the NVT simulation were used for metadynamics simulations with only the C-terminal and N-terminal residues restrained with 2 kcal/mol Å2 to prevent the drifting of the systems. The bin size was set to 0.1 Å for both the RCs. The Gaussian width (σ) was set to 0.8 Å (0.1 × 8 as implemented in NAMD code). The

τ = τM × α(t ) α (t ) ≈

(2)

Z0 = ⟨e β[V (s , t )]⟩M = (1/t ) ZM

∫0

t

dt ′e βV (s , t ′)

(3)

where V(s,t′) is the temporal biasing potential deposited in the basin region of the RC (here representing the binding pocket of the target) and was summed up at the initial position of the drug (bound state) and stopped when the drug dissociated 10 Å away from the initial position (Figure 1A3−G3). Here, β represents the inverse temperature (β = 1/kBT), where kB and T denote the Boltzmann constant and simulation temperature, respectively. We used the maximum biasing potential in each time step to calculate the residence time for the reason that the ligand will stay in the binding pocket until sufficient biasing potential is deposited in the binding basin. In the Bortolato algorithm, they proposed to characterize drug residence time with the maximum biasing potential derived from metadynamics simulations,27 which, in principle, can be represented into dissociation rate constant (koff) using Eyring’s equation koff =

kBT −ΔGoff /RT e h

(4)

where ΔGoff is the maximum biasing potential derived from metadynamics simulation, and kB, h, and R represent the Boltzmann constant, Planck’s constant, and ideal gas constant, respectively. All the maximum biasing potentials were determined as soon as the drug dissociated from the binding pocket (with 10 Å cutoff). The resulted koff values calculated by Eyring’s equation were then used to compare with the corresponding koff values determined by Tiwary’s algorithm. Mean Residence Time Calculation with Poisson Distribution. To calculate the mean residence time, one should derive it from several independent MD simulations. However, the average residence time may be susceptibly distorted with even one sample dropped into the long-time tail of the distribution.14,36,37 Alternatively, considering the dynamics is Markovian in nature, whose current state only depends on the last one, the occurrence of drug dissociation (known as rare events) should follow Poisson distribution, which can give the average transition time by avoiding incorporating the long-time tail influence into the statistics.37 1898

DOI: 10.1021/acs.jcim.7b00075 J. Chem. Inf. Model. 2017, 57, 1895−1906

Article

Journal of Chemical Information and Modeling

Figure 2. KS test between random samples and all samples, where the black line in each panel is the empirical CDF derived from the metadynamics simulation, and the red lines are the theoretical CDFs generated according to Poisson distribution.

theoretical CDF (TCDF, which is fitted from the ECDF with minimum variance, such as red curves in Figure 2) with the two-sample Kolmogorov−Smirnov (KS) test.37 Thereafter, we used CDFs for the residence time calculation and the simulation quality inspection. Moreover, to give a comparison, we also employed the medium value and the geometric average of the independent MD simulations for the residence time characterization.

To characterize the Poisson processes of rare events, the cumulative distribution function (CDF) is usually used to describe the probability of observing at least one time of the purposed event. Equation5 shows the CDF of a Poisson process to observe at least one dissociation/transition event within time t (here, t is replaced by log koff, thus the reverse cumulative curve is shown in Figure 2), where τ represents the average dissociation/transition time of the several independent MD simulations. With the simulation time exponentially growing, the probability of occurrence of a rare event will increase with an S-shaped curve up to 1 (or vice versa), implying that the event must happen within the time scale. ⎛ t⎞ Pn ≥ 1 = 1 − P0 = 1 − exp⎜ − ⎟ ⎝ τ⎠

■

RESULTS AND DISCUSSION Typical Simulation Results for Residence Time Calculation. As mentioned above, to calculate drug residence time, appropriate collective variables (CVs) should be selected at first. Here, we used the dissociation distance of drugs (characterized by drug−target distance) and conformational changes of protein pockets (characterized by RMSD changes upon the binding/unbinding process of drugs) as RC1 and RC2, respectively. Figure 1 illustrates the bound-state crystal structure (Figure 1A1−H1) and the corresponding typical metadynamics result (free energy landscapes in Figure 1A2−

(5)

Besides, it is an advantage to use CDF for the drug residence time characterization due to the reason that one may easily detect whether the simulation processes are Poisson-like by comparing the empirical CDF (ECDF, which is derived from the simulation result, such as black curves in Figure 2) and the 1899

DOI: 10.1021/acs.jcim.7b00075 J. Chem. Inf. Model. 2017, 57, 1895−1906

Article

Journal of Chemical Information and Modeling Table 1. Summarization of Calculated Residence Time Based on Various Strategies Target

Renin

EGFR

Neuraminidase

Drug

Lopinavirf

HIV Protease Atazanavirf

Indinavir

Aliskiren

Lapatinib

Oseltamivir

UTEh

Residence time (minimization, tmin, s−1)a

3.57 × 10−4 ±5.35 × 10−5g (−0.23)h

3.12 × 10−4 ±2.47 × 10−5 (−0.35)

6.89 × 10−5 ±4.67 × 10−6 (−1.37)

1.94 × 10−5 ±2.21 × 10−6 (−0.75)

4.09 × 10−9 ±2.20 × 10−10 (−3.98)

1.77 × 10−3 ±6.96 × 10−5 (0.85)

7.53

Residence time MD (median, tmed, s−1)b

5.16 × 10−4 (−0.10)

2.97 × 10−4 (−0.37)

6.97 × 10−3 (0.64)

3.86 × 10−6 (−1.46)

7.16 × 10−8 (−2.74)

0.91 (3.56)

8.87

Residence time MD (geometric average, tgeo, s−1)c

6.03 × 10−4 (−0.03)

2.48 × 10−4 (−0.44)

4.52 × 10−3 (0.45)

3.73 × 10−6 (−1.47)

4.49 × 10−8 (−2.94)

1.24 (3.70)

9.03

Residence time MD (Poisson distribution, τTiwary, s−1)d

7.46 × 10−4 (0.07)

5.57 × 10−4 (−0.09)

1.11 × 10−2 (0.84)

8.18 × 10−6 (−1.13)

4.77 × 10−7 (−1.92)

2.06 (3.92)

7.97

Residence time MD (maximum biasing potential, τBortolato, s−1)e

6.24 (3.98)

2.75 (3.60)

54.95 (4.54)

5.01 × 10−2 (2.66)

6.31 × 10−4 (1.21)

2.57 × 103 (7.01)

23.0

Experimental residence time (Koff, s−1)

6.50 × 10−4

6.90 × 10−4

1.60 × 10−3

1.10 × 10−4

3.90 × 10−5

2.50 × 10−4

−

PDB code

1MUI38

2AQU39

3WSJ40

2V0Z41

3BBT42

4HZX43

−

a

b

Residence time estimated based on minimized structures using Tiwary algorithm. Residence time estimated based on the median value of the 48 (or 96) MD structures using Tiwary algorithm. cResidence time estimated based on the geometric average of the 48 (or 96) MD structures using Tiwary algorithm. dResidence time estimated based on the Poisson distribution characterized transition time of the 48 (or 96) MD structures using Tiwary algorithm. eResidence time estimated based on the Poisson distribution characterized transition time of the 48 (or 96) MD structures using Bortolato algorithm. fProperties calculated based on two conformations of the crystallized drug. gStandard deviations derived from 10−20 Å of the reaction coordinate. hFold change between the predicted residence time and the corresponding experimental data (reported based on log10). i Unsigned total errors of the predicted residence time for the investigated systems.

final position of the drugs) is used to determine the release event of the drugs for the reason that all the investigated drugs have a radius less than 10 Å. Furthermore, it is reasonable because the calculated residence time is stable and independent of the cutoff values used in the determination of a release event (10−20 Å cutoff values were used for the validation, Figure S1 in the Supporting Information). As shown in Figure 1A3−H3, all the release events (regions around the red dotted lines) correspond to very steep dissociation curves, where the drugs unbind from the targets within a very short simulation time (at the level of several picoseconds) and obey the demand that there should be no biasing potential deposited in the transition region (which is also validated by the cutoff-independent residence times in Figure S1).35 Comparison of Residence Time Predicted by Optimized Structures and MD Trajectories Based on Tiwary Algorithm. Table 1 summarizes the predicted drug residence time (or dissociation rate, s−1) based on various calculation protocols (the detailed information can be found in Table S1). Due to the exponentially growing form of the accelerating factor (α, eq 3) that may show very long tails (unusual samples) far away from the normally sampled region, one should avoid calculating drug residence time by averaging a series of residence times derived from the MD structures directly. To solve this problem, we calculated the average drug residence time based on the median value (tmed), geometric average (tgeo), and Poisson distribution characterized transition time (τTiwary) (Table 1), which can effectively attenuate the influence of the long tails of the unusual samples in the data set.37 Thus, we could compare these values (namely, tmin, tmed,

H2 and dissociation trajectories in Figure 1A3−H3) based on the optimized structure for each drug−target complex, where two complexes, lopinavir (panels A1 and B1 in Figure 1)- and atazanavir (panels C1 and D1 in Figure 1)-bound HIV proteases, illustrate two opposite binding conformations of the drugs in the crystallized structures and were investigated individually in the following study. As mentioned before, to calculate the residence time, one only needs to study the first dissociation event rather than gaining the whole convergent free energy landscape. Therefore, the illustrated FESs are much more unconvergent and not suitable for the estimation of ligand−receptor binding free energy. Nevertheless, the depositing time-dependent biasing potential in the active basin (dark blue region in each FES) could accurately predict the drug residence time due to the exponential distribution of the biasing potentials (eq 3) that will be dominated by the potentials in the deep basin, and thus, a precise estimation of unbinding boundary is actually not necessary.35 The black dotted lines in each FES (Figure 1A2−H2) show the favorable unbinding pathways upon the dissociation of the drugs, where each drug−target complex selects an up and down dissociation pathway with the RMSD of the apo-state protein higher than the corresponding holo one, indicating that high-dimensional FES is necessary to avoid too high free energy barriers in a rare event and is consistent with the result of Salvalaglio.37 Panels A3−H3 in Figure 1 exhibit the time evolution of the first dissociation event of the drugs, where the red dotted lines represent the corresponding simulation time against the occurrence of the dissociation event for each complex. Here, a 10 Å cutoff (difference between the original position and the 1900

DOI: 10.1021/acs.jcim.7b00075 J. Chem. Inf. Model. 2017, 57, 1895−1906

Article

Journal of Chemical Information and Modeling tgeo, and τTiwary) directly to propose a better strategy for the prediction of drug residence time. To facilitate the following analysis, the residence time difference (Δt) between the predicted value (tpred) and the experimental data (texp) was adjusted from the exponential form to the log10 form (Δt = log10

t pred texp

(2AQU), whereas different residence times are shown for the two binding modes of lopinavir (1MUI). Nevertheless, the mean residence times (Poisson distribution characterized transition time based on 96 samples) derived from the two conformations of the ligands are very consistent with the experimental data, indicating that all the crystallized conformations of the ligands should be taken into consideration when calculating residence time. Although a relatively small number of systems were used in this comparison, the conclusion is consistent with our previous work on the assessment of the MM/GBSA method that, compared with those predicted based on long-time MD structures, the minimized structures can usually give a better correlation between the experimental data (for both crystallized structures and docking structures).62−64 The deeper reason may also be attributed to the error of the force field, where any error of the MD structures may significantly affect the predicted residence time. Comparison of Residence Time Predicted by Tiwary and Bortolato Algorithms. As has been mentioned above, recently Bortolato and co-workers proposed to characterize drug residence time with the maximum biasing potential derived from metadynamics simulations.27 Thus, we give a comparison of the two methods. As shown in Figure 3, a very

, where each Δt is reported in the bracket

below the residence time in Table 1). Moreover, the unsigned total error (UTE) for each prediction strategy is reported in Table 1 as well. It can be found in Table 1 that the UTE based on the minimized structures is the lowest among the four strategies (UTEmin = 7.73 unit, UTEmed = 8.87 unit, UTEgeo = 9.03 unit, and UTEτ_Tiwary = 7.97 unit), implying that the minimized structures are favorable in the prediction of drug residence time. Besides, the overall accuracy of the predicted residence time derived from the Poisson distribution (τTiwary) is also better than those based on the median values and geometric averages, suggesting that it is necessary to estimate drug residence time from the Poisson distribution when using a large number of MD structures. Detailed analysis shows that four out of six optimized systems (1MUI, 2AQU, 2V0Z, and 4HZX) exhibit |Δt| less than 1 (namely, the error of the predicted value is within 10 folds of the experimental data), whereas three systems (1MUI, 2AQU, and 3WSJ) based on the MD trajectories show |Δt| less than 1, meaning that the strategy of using minimized structures for residence time estimation may be very promising to be used in large scale VS. Nevertheless, the predicting error of the system of 3WSJ for the minimized structures and 2V0Z for the MD structures are still within 100 fold of the experimental data (|Δt|