Multiscale Simulation of ReceptorâDrug Association Kinetics

Article pubs.acs.org/JCTC

Multiscale Simulation of Receptor−Drug Association Kinetics: Application to Neuraminidase Inhibitors Fabian Zeller, Manuel P. Luitz,† Rainer Bomblies, and Martin Zacharias†,* Physik-Department T38, Technische Universität München, James-Franck-Str. 1, 85748 Garching, Germany S Supporting Information *

ABSTRACT: A detailed understanding of the drug−receptor association process is of fundamental importance for drug design. Due to the long time scales of typical binding kinetics, the atomistic simulation of the ligand traveling from bulk solution into the binding site is still computationally challenging. In this work, we apply a multiscale approach of combined Molecular Dynamics (MD) and Brownian Dynamics (BD) simulations to investigate association pathway ensembles for the two prominent H1N1 neuraminidase inhibitors oseltamivir and zanamivir. Including knowledge of the approximate binding site location allows for the selective confinement of detailed but expensive MD simulations and application of less demanding BD simulations for the diffusion controlled part of the association pathway. We evaluate a binding criterion based on the residence time of the inhibitor in the binding pocket and compare it to geometric criteria that require prior knowledge about the binding mechanism. The method ranks the association rates of both inhibitors in qualitative agreement with experiment and yields reasonable absolute values depending, however, on the reaction criteria. The simulated association pathway ensembles reveal that, first, ligands are oriented in the electrostatic field of the receptor. Subsequently, a salt bridge is formed between the inhibitor’s carboxyl group and neuraminidase residue Arg368, followed by adopting the native binding mode. Unexpectedly, despite oseltamivir’s higher overall association rate, the rate into the intermediate salt-bridge state was found to be higher for zanamivir. The present methodology is intrinsically parallelizable and, although computationally demanding, allows systematic binding rate calculation on selected sets of potential drug molecules.

■

trypsin.7,8 However, using these approaches, large sampling efforts are spent on the diffusive search of the binding site. For these diffusion controlled regimes, Brownian dynamics (BD)9 simulations are an efficient representation and have been previously used to simulate entire association pathways,10−16 protein−protein association kinetics, the influence of electrostatic guiding,17 and protein crowding.18 However, BD simulations often ignore or simplify the flexibility of the receptor and ligand. Also, effects due to explicit solvent molecules are neglected and short-range interactions between receptor and ligand may only be approximately included. It typically requires the prior definition and optimization of reaction criteria to achieve good agreement with experiment.17,19 Advanced sampling atomistic explicit solvent MD methods have been developed to efficiently sample putative binding sites for drug molecules on protein molecules based on metadynamics20,21 or Hamiltonian replica exchange methods.3,22 However, these approaches do not allow a direct extraction of binding kinetics. Inspired by previous work on calculating association rates,13,23−26 we present an extended multiscale approach that

INTRODUCTION Until recently, determining and optimizing equilibrium binding affinities to increase the inhibition efficacy of potential drug molecules has been at the focus of most drug design efforts.1 There is, however, growing evidence that drug efficacy is also determined by binding kinetics.2 Thus, the detailed association pathways play an increasingly important role in drug design. Knowledge of the complete pathway ensemble may additionally reveal alternative binding modes which remain disregarded when optimizing the affinity of a single binding pose.3,4 The atomistic simulation of complete association pathways is, however, still computationally demanding. In contrast to the calculation of equilibrium quantities, which can be obtained from enhanced sampling methods, the determination of kinetics requires simulation at time scales of the physical association process. Ultralong simulations of milliseconds on specialized hardware at atomistic detail showed that ligand− receptor complex formation can be covered within molecular dynamics (MD) simulations without including prior knowledge of binding sites and reveal the principles of binding mechanisms.5,6 As an alternative, the application of Markov state models allows for the combination of independent short MD simulations. For a charged enzyme inhibitor, this approach has been used to investigate the association of benzamidine to © 2017 American Chemical Society

Received: June 17, 2017 Published: August 18, 2017 5097

DOI: 10.1021/acs.jctc.7b00631 J. Chem. Theory Comput. 2017, 13, 5097−5105

Article

Journal of Chemical Theory and Computation

binding processes,10−16 our method accounts for full flexibility of the receptor and ligand and includes explicit solvent molecules in the region where the receptor and ligand commence to strongly interact.

combines Molecular Dynamics (MD) and Brownian Dynamics (BD) simulations to efficiently investigate the association process of drug molecule candidates to a target receptor binding site on readily available hardware. The approach limits the computationally expensive atomistic MD simulations to the vicinity of a binding site of interest which is often already known (Figure 1, blue area) and connects the MD regime to

■

MATERIALS AND METHODS Molecular Dynamics Simulations, System Preparation and Equilibration. To avoid any bias toward a bound conformation of the binding site, a neuraminidase apo structure was used for the association simulations (monomeric apo H1N1 neuraminidase from PDB: 4B7M30) which contains three point mutations. The wild type amino acid sequence of H1N1 Cal07 was restored by in silico reverse mutating residues Ile106Val, Arg223Ile, and Asp248Asn in the starting structure. Residual phosphate ions and monosaccharides were removed from the structure, but the crystallographic water and one CA2+ ion bound close to the receptor site was conserved. To reduce the computational overhead only one monomer of the tetrameric neuraminidase complex was used. The protonation states of neuraminidase at neutral pH were calculated with the karlsberg+ software package31,32 for both complexes and kept constant during simulations at neutral pH (see Supporting Information Table S5). For comparison with the bound conformation, reference crystal structures of oseltamivir and zanamivir bound to H1N1 neuraminidase were taken from PDBs 4B7Q and 4B7R.30 The MD simulations were carried out with the pmemd.cuda33 module of the Amber14 package34 at atomistic detail,35,36 including explicit solvent (at neutral pH). The parmff99SB-ILDN force field37 was used for the protein and topologies. For the drug molecules, oseltamivir/zanamivir (structures indicated in Supporting Information Figure S3) were parametrized using the GAFF force field38 as part of the antechamber package.39 Calcium ion parameters were taken from Bradbrook et al.40 The protein was solvated in 8214 (zanamivir) and 8150 (oseltamivir) molecules of the TIP3P41 explicit water model in an octahedral periodic box. Sodium and chloride ions were added in order to render the overall charge neutral. Short range van der Waals and electrostatic interactions were cut off beyond 9 Å; the Particle Mesh Ewald42 (PME) method was used for long-range interactions. The Berendsen thermostat and barostat43 were used to control the temperature (310 K) and pressure (1.01 bar) of the simulations. For equilibration, energy minimization was performed for 2000 steps with the steepest descent algorithm. Then the system was heated and equilibrated in the NPT ensemble for 2 ns with all solute heavy atom positions restrained by a harmonic potential. Drug Molecule Association Simulations. During drug association simulations the hydrogen mass repartitioning scheme36 and the SHAKE method35 were applied to the nonsolvent hydrogen atoms. This allowed using an MD integration time step of 4 fs. To obtain a Canonical ensemble the simulations were performed in the NVT ensemble using the Andersen thermostat44 with a collision frequency of 2 ps−1 and a randomization of the velocities to T = 310 k every 1000 steps. For the association simulations, the systems were simulated for a total of 85.7 μs (50.0 μs for oseltamivir and 35.7 μs for zanamivir). Bound Mode Determination via Clustering. Stable bound modes were identified through clustering with respect to the RMSD (of the ligand) between sampled states. For each ligand, all trajectories that resulted in a bound drug molecule were concatenated. The concatenated trajectory with a time step of 1 ns was superimposed with respect to the rigid, β-sheet

Figure 1. Schematic depiction of the simulation setup. The isotropic rate from infinity into a sphere of radius 60 Å (dotted circle) can be obtained by continuum theory. The diffusion and long-range electrostatics dominated part of the association pathways from the sphere into the encounter surface (ES, red) was treated with Brownian dynamics simulations. Close range interactions were covered by MD simulations with explicit water representation, starting from a separately generated atomistic ensemble at the ES. The propagation was terminated either when the ligand left the MD region (blue segment, defined as cone extending up to 32 Å form protein center-ofmass; see Supporting Information, part 1) or after a minimal residence time of 2 μs. MD trajectories in which the ligand left the MD region were continued as an ensemble of BD trajectories to determine their probability to re-enter the ES.

less demanding BD simulations in the diffusion controlled regime far away from the receptor (Figure 1, dashed sphere up to red encounter surface). We evaluate the method on two prominent, clinically relevant inhibitors (oseltamivir and zanamivir) of an important drug target, influenza H1N1 neuraminidase. Detailed association pathways for oseltamivir and zanamivir were obtained, revealing a single, common intermediate state. A higher overall binding rate for oseltamivir is predicted, coinciding with experimental findings. Surprisingly, however, the rate into the intermediate state is higher for zanamivir. The confinement of the MD simulations to the vicinity of a binding site allows for long simulation of the relevant parts of the association pathway. We also compare the application of a more generic binding criterion in the form of a minimal residence time with a geometric criterion based on the similarity to a known bound structure. A residence time criterion eliminates the necessity of reaction criteria to be fitted to prior knowledge about the binding event, e.g. a X-ray structure, as used in previous studies which relied on geometric definitions of the bound state.12,13,19,23,27−29 The method does not bias sampling to certain regions of the association pathway space, nor does it require identification of transition states. In contrast to e.g. ultralong MD simulations or Markov state model approaches,7,8 the major part of the sampling effort is not spent on the diffusive search of the binding site. Unlike the sole application of BD simulations for the investigation of 5098


Article


Table 1. Rate from Infinite Distance into Encounter Surface (ES) Obtained from Continuum Theory and BD Simulations MD BD (k∞→ ES ), Number of Binding MD Trajectories over Total Number of MD Trajectories Started at the ES (αES → bound ), Total calc. exp. 30 a Calculated On-Rate (kon ) and Experimental On-Rate (kon ) for Both Ligands Molecule

kBD ∞→ES [1/μMs]

αMD ES→bound, 2 μs

kcalc. on [1/μMs]

αMD ES→bound, 1 Å

kcalc. on [1/μMs]

αMD ES→bound, 2 Å

kcalc. on [1/μMs]

kexp. on [1/μMs]

Oseltamivir Zanamivir

794 [776,822] 737 [519,754]

19/676 4/606

25.5 5.5

20/676 4/606

27.9 5.5

22/676 25/606

30.7 35.5

3.1 1.3

The calculated rates are given for different definitions of the bound state: A minimal residence time of 2 μs and a minimal RMSD threshold of 1 and BD 50 2 Å. The statistical error of k∞→ ES was obtained by the bootstrap method in the form of 95% confidence intervals. A statistical error estimate of MD αES → bound and thus the total on-rate kon is not meaningful due to the small number of sampled events. a

Brownian Dynamics Simulations. The MD simulations starting at the vicinity of the binding site were complemented by Brownian Dynamics (BD) simulations starting outside of this region. For each drug molecule, 106 BD trajectories were started from a spherical shell of radius b = 60 Å with the receptor at the center to determine the probability of arriving at the ES before diffusing to distances larger than q = 100 Å (see Supporting Information Figure S2). Note that test calculations indicated that beyond the spherical shell of radius b = 60 Å the interaction of the protein and neutral ligand can be neglected (for a charged ligand this may not necessarily be the case). Also, 106 BD trajectories were started from the MD configurations terminated outside of the MD area to determine the probability to re-enter the vicinity of the binding site before diffusing to distances larger than b. From the isotropic rate into a sphere of radius b, the results of the BD simulations, and the fraction of binding MD trajectories, the on-rate was calculated using the method of Northrup et al.13 The BD simulations were performed using an in-house program based on the Ermak− McCammon algorithm.47 In the BD simulations, the drug molecule and the receptor are treated as rigid structures that undergo translational and rotational diffusion13 in an implicit continuum solvent.9 Electrostatic forces on the drug molecule in the field of the receptor were taken into account by solving the Poisson−Boltzmann equation with the APBS software48 in order to calculate the electrostatic potential of the protein (dielectric constant ϵ = 2.0) embedded in an aqueous environment (dielectric continuum with ϵ = 80.0) at a temperature of 310 k. The influence of the drug molecule on the electrostatic field was neglected. Changes in the ion concentration (0 to 0.2 M) in the FDPB resulted in negligible effects on the rates in the BD regime (presumably because the ligands are overall charge neutral). Hydrodynamic interactions between the protein and the drug molecule were neglected (further details are given in Supporting Information, section 2). In order to prevent steric clashes between ligand and receptor molecules a configuration was defined as overlapping, when the shortest distance between any atom of the ligand and any atom of the receptor was below a collision radius of 2 Å. BD propagation steps which resulted in such overlapping configurations where then rejected. The diffusion coefficients were calculated according to the Stokes−Einstein equation, using the radius of gyration of the molecules as resulting from the MD simulations (see Supporting Information Table S4). The propagation time step was 100 ps. Despite the large number of simulated trajectories, due to the simplicity of the model, the computational cost of the BD simulations is negligible in comparison to the MD simulations.

forming protein backbone part of neuraminidase. For each snapshot the RMSD of the ligand with respect to all other snapshots was calculated. The snapshot with the largest number of neighbors, i.e. snapshots within an RMSD 1.5 Å, formed the largest cluster and, with all its neighbors, was removed from the pool of structures. This process was repeated until all structures were assigned to a cluster, so that every structure was a member of exactly one cluster. The implementation in GROMACS v4.6.545 was used, based on the approach by Daura et al.46 Association Pathways and On-Rate Calculation. The diffusion dominated part of the association pathway (simulated using the Brownian Dynamics method; see below) was connected to the computationally expensive Molecular Dynamics (MD) simulations at an encounter surface (see Figure 1, red circle) in the vicinity of the binding site. The surface was defined by a 12 Å distance of the center of mass of the Cα atoms of Asn248 and Pro431. This point is located approximately at the center of the entrance to the neuraminidase binding site, and the sterically accessible regime is defined by a cone (illustrated in Figure 1; see also Supporting Information section 1 and Tables S1, S2 and Figure S1 for a detailed definition of the encounter surface). At the core of this study are atomistic MD simulations with an explicit solvent representation starting from the encounter surface. The selected encounter surface is close enough to the binding region to allow studying ligand binding efficiently using atomistic MD simulations in explicit solvent. At the same time, for the starting positions at the encounter surface, the placement and orientation of the drug molecules are sterically unrestricted. To obtain an unbiased ensemble of configurations at the encounter surface (ES) starting structures for the MD simulations were generated by placing the drug molecules in the MD region at distances of 31 Å, which is several angstroms outside the ES. Trajectories were propagated from these starting structures using a nondeterministic Langevin dynamics integrator with a collision frequency of 5 ps−1. As long as the respective drug molecule resided within the MD region, the configurations from these trajectories were saved every 10 ps. This procedure was repeated until ∼5000 configurations had been generated for both inhibitors. The generated configurations were propagated until leaving the MD region or entering the ES. This resulted in an ensemble of ∼600 configurations for each drug molecule with the drug molecule located at the ES. The total simulation time required for the generation of the ES ensembles sums up to ∼100 ns. Further association simulations were started from this ensemble. Each of these drug association simulations was propagated until leaving the MD area (Figure 1, blue area) or after a maximum simulation time of 2 μs. Trajectories staying within the MD area for 2 μs were considered to represent bound inhibitors. 5099


Article


Figure 2. (Left) Number of trajectories within binding region (MD regime) vs simulation time for zanamivir and oseltamivir association simulations. (Right) Number of successful trajectories reaching a bound state defined by an RMSD (ligand) threshold vs RMSD for the zanamivir and oseltamivir inhibitor case.

■

RESULTS A combined BD and MD approach was employed to calculate the association rates for oseltamivir and zanamivir inhibitors (chemical structures given in Supporting Information Figure S3) to the neuraminidase target enzyme. Knowledge of the approximate location of the binding site was included in the form of a geometrically defined encounter surface (ES) (see Figure 1). The ES is designed to completely extend over the entire binding site entrance, such that it has to be crossed by approaching inhibitors in order to bind. Far away from the binding site, up to the ES, the association pathways are represented by coarse BD simulations. Starting from the ES, the association pathways are represented by atomistic explicit solvent MD simulations. Thus, the ES is ideally located in the region where specific short-ranged interactions between inhibitors and binding site start to play a role. Here, the ES was defined by a distance of 12 Å between the inhibitors and the approximate center of the binding site entrance (see Methods and Supporting Information). Exemplary trajectory snapshots of inhibitors at the ES are illustrated in Figure 3A. The steady-state rate of oseltamivir and zanamivir from BD infinite distances into the ES k∞→ ES (Table 1) was calculated from the isotropic continuum rate into a sphere surface of radius 60 Å and the fraction of BD trajectories starting at this sphere surface and, subsequently, reaching the ES (see Methods). As expected for the two inhibitors with similar diffusion constants, at distances outside the ES where the binding site does not directly interact with the inhibitors, the rates are roughly equal (see Table 1). Note that for this part the rates were also minimally affected by changing salt concentration or other parameters in the Poisson−Boltzmann electrostatic calculations to obtain forces on the (neutral) drug molecules (see Methods). The transitions from the encounter surface into the actual binding site were studied with Molecular Dynamics (MD) simulations. MD trajectories were started with an inhibitor at

the encounter surface and terminated when it left the MD area (Figure 1, blue area). Several options to define ligand bound states were considered. Frequently, the root-mean-square deviation (RMSD) from a known bound location at the receptor binding site is used as a criterion for defining a bound state (RMSD is here the RMSD of the ligand heavy atoms after best superposition with respect to the rigid, β-sheet forming protein backbone part of neuraminidase).19,28,29 However, in this case the results may depend on the RMSD threshold and possible force field artifacts that may favor a bound state (with appropriate contacts) that differs from the experimentally determined structure. In experiments the RMSD of the ligand from a bound state cannot be obtained and a bound state is typically detected by an accessible physicochemical property that changes upon complex formation (and stays usually constant for the rest of the measurement time). We compared definitions based on a whole range of possible RMSD thresholds with a definition based on residence time (in principle equivalent to the experimental setup albeit on a smaller time scale given by the maximum simulation time). In this case trajectories were considered bound when the inhibitor resided in the MD region for at least 2 μs corresponding to the maximum MD simulation time (see Methods). This binding criterion has the additional advantage that it is independent of any geometric definition of the bound state which can be very sensitive to a chosen threshold or force field accuracy. However, as a drawback it can be affected by the maximum simulation duration. In order to check the convergence we plot the number of trajectories that stay within the binding region vs simulation time. This should converge to an asymptote reflecting the true binding kinetics (Figure 2). The majority of MD trajectories leaves the MD region after less than 0.1 μs, and only a small fraction enters the binding sites and required longer simulation times. In the case of oseltamivir, the number of trajectories that stays in the binding region converges relatively rapidly to around 20 reaching 19 after the maximum 5100


Article

Journal of Chemical Theory and Computation simulation length of 2 μs (Figure 2). If one instead uses an RMSD threshold criterion for binding, the number of successful binding trajectories depends strongly on the choice of the threshold. Indeed, for an RMSD threshold of 1−2 Å the number of successful trajectories is very close to the result based on residence time in the MD regime (see above) but a choice of 3 Å gives already a significantly larger number of binding events. In the case of zanamivir, the result is more complex; after 2 μs, only four drug molecules transition into stable binding modes and remain in the MD regime. When using an RMSD based binding criterion, the resulting successful binding trajectories vary considerably even for the interval between an RMSD of 1−2 Å (Figure 2). As we demonstrate below by simulations starting from the bound form, the strong RMSD threshold dependence and loss of bound states is most likely due to a force field artifact and not a principle problem of our simulation approach. The result demonstrates that if the drug molecule approaches a stable bound state within the maximum MD simulation time (as is the case for oseltamivir), the evaluation of successful trajectories based on the maximum residence time criterion can reach a converged asymptote and agrees also with a criterion based on deviation from the known native binding mode. However, the latter criterion requires precise knowledge of the native binding geometry not always accessible in the case of drug design/screening efforts. MD trajectories terminated at the border of the MD region were continued as an ensemble of BD trajectories to determine the probability of rearriving at the encounter surface. The total on-rate kon (Table 1) was calculated from the BD simulations MD and the fraction of bound MD trajectories αES → bound by the 13 method of Northrup (see Methods, SI). The calculated total association rates are based on the number of trajectories reaching the bound state and, for the minimal residence time definition, coincide with the experimental ranking of the inhibitors. However, the absolute values are overestimated by a factor of 4−8. Exemplary inhibitor configurations serving as starting points for the MD simulations at the encounter surface are shown in Figure 3A. By analyzing all binding trajectories of both oseltamivir and zanamivir, the formation of the salt bridge between the carboxyl group of the drug molecule and Arg368 of neuraminidase (distance smaller than 3 Å, Figure 3B) was identified as a common first transient interaction with the binding site. This key interaction was observed in a variety of orientations of the ligand relative to the binding site and was only restricted sterically by the shape of the binding entrance region. The probability of the MD trajectories to form the carboxyl salt bridge was found to be higher for zanamivir (101 of 606) than for oseltamivir (64 of 676), predicting a higher rate into the salt-bridge ensemble for zanamivir. This can be explained by the orientational bias caused by the higher electrostatic dipole moment of zanamivir (28.3 D) interacting with the electrostatic field of the positively charged binding site in comparison to oseltamivir (23.0 D) (see Figure 4). Although zanamivir adopted the intermediate salt-bridge state more frequently, it reached final bound modes with the same (RMSD threshold = 2 Å) or a lower probability (binding based on residence time) than oseltamivir, agreeing with the experimentally determined ratio of on-rates. Thus, the subsequent specific short-range interactions and detailed binding site rearrangements must be the reason for a higher overall on-rate for oseltamivir compared to zanamivir. For the

Figure 3. Predominant association pathways for oseltamivir (orange sticks) and zanamivir (green sticks) to the binding site (light blue area) of neuraminidase (white). MD simulations were started from ligand positions at the ES (A, exemplary ligand positions). In all trajectories that resulted in binding, the ligand’s carboxyl group first formed a key contact with Arg368 (B). The largest clusters from bound trajectories (C) are depicted with population percentage and ligand RMSD from X-ray structure (blue numbers). Colored arrows indicate transition probabilities from ES (A) via initial salt bridge state (B) to bound modes (C). Numbers associated with these arrows represent the number of trajectories that started from or reached the respective state. Bound structures from X-ray crystallography30 are shown on the left. Important binding site residues of neuraminidase are shown in white stick representation.

transition from the salt-bridge ensemble into the stable modes, no distinct predominant intermediate states could be identified. The 150-loop (residues 146−152), known to be very flexible,49 remained in a closed conformation in all oseltamivir binding trajectories. For zanamivir in two of the 2 μs trajectories, the flexible 150-loop opened up (shown in Figure 5). To analyze the most stable binding modes indicated by the simulations, the bound trajectories were clustered with respect to ligand RMSD (RMSD of ligand after superpositioning the rigid β-sheet part of the receptor onto the native complex). In Figure 3C the most populated clusters are shown with their respective population. For oseltamivir, in nearly 80% of the frames the inhibitor adopted a binding geometry very similar to the crystal structure, with an RMSD of only 1.0 Å (see Co1). This pose was observed to exchange with two secondary binding modes. The diethane group was frequently observed to flip out of the corresponding binding cavity (Co2) and vice versa. For the reverse transition (Co2 to Co1), the opening of the Asn295 side chain appears to be crucial. Also, the carboxyl group detaches and the salt bridge occasionally opens (Co3), 5101


Article


Indeed, simulations of zanamivir starting from the crystal structure resulted in similar deviations of the sampled ligand placements from the structure in the crystal. In addition, ligand binding states close to the configurations found in our association simulations (started from the apo structure of the enzyme) were also sampled (see Figure S4 in SI). This result indicates that the deviation of some of the sampled zanamivir bound configurations from the native state is not due to insufficient sampling (trapping in configuration on the path to the bound state) but also occurs rapidly when starting from the bound structure. It is important to emphasize that this does not invalidate our approach; rather, it indicates to interpret results with care if the ligand-bound states start to differ appreciably from the native bound state.

■

Figure 4. Illustration of electrostatic properties of neuraminidase and inhibitors. The surface potential (color coded, blue positive, red negative) and electrostatic field lines of neuraminidase were obtained from Poisson−Boltzmann calculations.48 The field lines guide the two ligands oseltamivir (orange) and zanamivir (green) into the positively charged binding site according to their respective dipole moment represented as arrows.

DISCUSSION During association processes, the initial approach of a ligand toward the receptor is governed by a diffusive random walk and possibly guided by long-range electrostatic interactions. This regime can be accurately represented by BD simulations,9 treating molecules as rigid bodies with their translational and rotational movement modeled as a diffusion process under the influence of electrostatic forces. The solvent is represented by a stochastic force, and effects of the molecular solvent structure are neglected. These simplifying assumptions allow far longer simulation times compared to MD, but break down when the ligand comes into contact with the binding site. In order to extract association rates from such approaches, in previous studies, geometric reaction criteria were applied.13,19,23,29 However, as resulting kinetics are highly sensitive to the choice of these criteria, the predictive value is limited. For the binding site vicinity, a more detailed model with atomistic resolution and an explicit representation of the solvent, shown to significantly influence ligand binding mechanisms,51,52 is advantageous. Atomistic explicit solvent MD simulations provide such a high level of detail at the cost of increased computational demands. Multiscale approaches have been proposed that combine the BD simulations in the diffusive regime with expensive but more accurate MD simulations for close range interactions.23,24 In combination with milestoning such methodology has recently been applied to study the binding kinetics of benzamidine to trypsin in the SEEKR approach (requiring 19 μs simulation time).26 The same system with a charged benzamidine ligand approaching an enzyme active site has also been used in multiple MD association simulations,7 using a total simulation time of around 50 μs, which is in the same order as the simulation time we employed for a single ligand. In the present study, we employed a combination of BD and explicit solvent MD simulations to investigate the association of two clinically relevant overall neutral inhibitors, oseltamivir and zanamivir, to H1N1 neuraminidase. It is important to note that the association rate constant for these ligands is roughly 10 times smaller than for the benzamidine/trypsin test system (and strong electrostatic guiding due to the charge complementarity between benzamidine and trypsin active site) indicating a considerable saving of computer time of our BD/MD compared to an approach based entirely on MD simulations. We also evaluated the option to obtain the association rate from the asymptote of the number of drug trajectories that do not leave the binding region within the maximum simulation time. The advantage of this approach is that only knowledge of the rough location of the binding site is required and all simulations can be started

Figure 5. Opening of the 150-loop of neuraminidase with bound zanamivir. The closed loop conformation of the X-ray crystal structure30 (light red cartoon) is shown in comparison with a sampled open loop (blue cartoon) conformation. The bound X-ray crystal structure of zanamivir (green) and contacting side chains of the loop are depicted as sticks; the binding site is colored in light blue.

with all other binding features remaining closed. The 150-loop remained closed for all binding modes similar to the X-ray structure. In the zanamivir simulations also binding geometries very similar to the crystal structure were sampled. Nevertheless, the most frequently observed modes were deviating from the crystal structure pose by more than 3.7 Å (Figure 3). Despite the larger deviation, important key contacts between zanamivir and the enzyme similar to the crystal structure were formed. In modes Cz1 and Cz5 (Figure 3), the carboxyl group of zanamivir formed a stable salt bridge with Arg118 and the guanidinium group formed a salt bridge with Asp151. In modes Cz2 and Cz3, Arg118 and Asp151 were displaced as a result of the partially opened 150-loop, not allowing interactions with zanamivir. The carboxyl Arg368 salt-bridge was formed. Mode Cz4 shows a closed carboxyl-Arg368 salt bridge with the glycerol group of zanamivir pointing in the direction of the 150loop. We suspect that for zanamivir the significant dependence of the calculated association rates on the RMSD threshold for binding and the loss of bound states even after long simulation times (see above) might be caused by force field artifacts. 5102


Article


ensemble of receptor starting configurations or consider a BD approach with a flexible receptor structure. It is important to emphasize that the assumption that the drug molecules can reach a stable native bound state within the maximum explicit MD simulation time is central to our approach when estimating the number of successful trajectories from the number of trajectories that resided in the binding region. In such case with a large number of starting configurations, it is possible to extract the association rate without any detailed definition of binding criteria and only an approximate knowledge of the binding site location. In the case of activation-controlled systems, without knowing the binding site, the application of the present approach might be difficult because the determination of an appropriate maximum simulation time may require performing several test simulations and a convergence analysis. The present methodology of efficiently combining BD and MD simulations is optimal in the case of diffusion controlled binding reactions. However, in cases where the final association of the ligand toward the native binding mode involves crossing large energy barriers or large scale conformational changes and dominates the association process (activation control), it might be less useful. In such cases the computational savings of using the BD approach for the initial diffusive encounter may play only a minor role at the cost of much longer explicit solvent MD simulations required to reach the final bound configuration. Even in the present case although the majority of MD trajectories starting at the encounter surface are discontinued after leaving the binding region, the approach is still computationally demanding mostly due to the small number of long simulations that enter the binding region. However, atomistic resolution of the association pathway in the vicinity of the binding site is of significant importance to account for side chain rearrangements and solvent molecule effects, while in the bulk state a coarse grained representation is sufficient. With steadily increasing hardware performance applied to selected sets of drug candidates, this method may complement available techniques for calculating equilibrium quantities in the drug discovery process by an efficient and highly parallelizable approach for investigating the dynamics of binding processes. From the obtained kinetic onrates, kinetic off-rates or the average thermodynamic residence time can also be obtained when the equilibrium binding affinity is known. Efficient free energy methods are available to calculate such absolute binding free energies.53,54

from the apo form of the receptor to avoid a bias toward the bound form. This avoids geometric definitions requiring previous assumptions about the binding, and thus, in contrast to previous approaches, no fitting of a binding criterion to experimental data is required.13,19,28,29 Indeed, for oseltamivir this criterion gave very good agreement with an RMSD based binding criterion if an RMSD threshold of 1−2 Å was used. For the zanamivir ligand the sampled placements in the binding region showed larger deviations from the native structure (weaker binding presumably due to force field artifacts), and therefore, the agreement between rates obtained with both binding criteria showed larger differences. Nevertheless, the two binding criteria (residence time and ligand RMSD) are consistent when selecting a ligand RMSD of 1 Å resulting in the same number of bound events for both ligands (∼20 for oseltamivir, ∼4 for zanamivir). This consistency check with respect to the established RMSD binding criterion further confirms the validity of our approach using only the residence time of the ligand in the binding site as a binding criterion for the present cases. BD simulations were used to estimate the rate into the ES. From the ES, MD simulations were started and run only in the vicinity of the known binding site. As expected for two electrically neutral inhibitors with similar diffusion constants, BD simulations yielded similar association rates into the encounter surface (Table 1). Due to the change to an atomistic water representation, BD simulations were not directly continued at MD resolution. Instead, MD simulations were started from an independently generated ensemble at the ES (red sphere shell in Figure 1). The generation of an MD based ensemble in the presence of explicit solvent results in a more realistic probability distribution of drug molecule starting placements at the encounter surface compared to an ensemble based on the more approximate BD simulations. Our MD simulations reveal one single, essential intermediate state in which a salt bridge between the carboxyl group of the inhibitors and Arg368 of neuraminidase is formed. From the salt-bridge state, further association proceeds as a broad ensemble of pathways and is governed by complex short-ranged interactions and rearrangements of the binding site. Interestingly, for the association processes of other receptor-inhibitor complexes, simulation studies indicated multiple distinct intermediate states.8 Although zanamivir has a higher association rate toward the salt-bridge state due to the guiding effect of its higher dipole moment, oseltamivir has a higher overall on-rate (using the residence time criterion for indicating bound states). The probability of oseltamivir reaching a bound state from the intermediate state is about 10 times higher than that of zanamivir. This relative result is in good agreement with experiment and highlights the importance of full atomistic coverage of the final steps of association. In the case of zanamivir, large scale rearrangements of the 150-loop, also found to be flexible in experiments and simulation,30,49 affected key binding features (Figure 5). In contrast, loop flexibility seems to play only a minor role in oseltamivir binding. It underpins that relevant conformational rearrangements need to be accessed by the simulation when extracting binding kinetics, as they can specifically influence the binding of different inhibitors. The closed loop conformation is overrepresented in the starting ensemble, as all simulations start from the same crystal structure receptor conformation. In future studies, one could consider using an equilibrated

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.7b00631. Detailed information on the definition of the encounter surface, association rate calculations, and dynamics of zanamivir in complex with neuramidase (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Martin Zacharias: 0000-0001-5163-2663 Present Address †

Center for Integrated Protein Science Munich, 81377 München, Germany.

5103


Article

Journal of Chemical Theory and Computation Notes

(18) Mereghetti, P.; Gabdoulline, R. R.; Wade, R. C. Brownian dynamics simulation of protein solutions: structural and dynamical properties. Biophys. J. 2010, 99, 3782−3791. (19) Sung, J. C.; Wynsberghe, A. W. V.; Amaro, R. E.; Li, W. W.; McCammon, J. A. Role of Secondary Sialic Acid Binding Sites in Influenza N1 Neuraminidase. J. Am. Chem. Soc. 2010, 132, 2883− 2885. (20) Söderhjelm, P.; Tribello, G. A.; Parrinello, M. Locating binding poses in protein-ligand systems using reconnaissance metadynamics. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, 5170−5175. (21) Cavalli, A.; Spitaleri, A.; Saladino, G.; Gervasio, F. L. Investigating drug−target association and dissociation mechanisms using metadynamics-based algorithms. Acc. Chem. Res. 2015, 48, 277− 285. (22) Ostermeir, K.; Zacharias, M. Accelerated flexible protein-ligand docking using Hamiltonian replica exchange with a repulsive biasing potential. PLoS One 2017, 12, e0172072. (23) Luty, B. A.; El Amrani, S.; McCammon, J. A. Simulation of the bimolecular reaction between superoxide and superoxide dismutase: synthesis of the encounter and reaction steps. J. Am. Chem. Soc. 1993, 115, 11874−11877. (24) Votapka, L. W.; Amaro, R. E. Multiscale estimation of binding kinetics using molecular dynamics, brownian dynamics, and milestoning. J. Biomol. Struct. Dyn. 2015, 33, 26−27. (25) Votapka, L. W.; Amaro, R. E. Multiscale estimation of binding kinetics using Brownian dynamics, molecular dynamics and milestoning. PLoS Comput. Biol. 2015, 11, No. e1004381. (26) Votapka, L. W.; Jagger, B. R.; Heyneman, A. L.; Amaro, R. E. SEEKR: Simulation Enabled Estimation of Kinetic Rates, A Computational Tool to Estimate Molecular Kinetics and Its Application to Trypsin−Benzamidine Binding. J. Phys. Chem. B 2017, 121, 3597− 3606. (27) Alsallaq, R.; Zhou, H.-X. Electrostatic rate enhancement and transient complex of protein−protein association. Proteins: Struct., Funct., Genet. 2008, 71, 320−335. (28) Le, L.; Lee, E. H.; Hardy, D. J.; Truong, T. N.; Schulten, K. Molecular Dynamics Simulations Suggest that Electrostatic Funnel Directs Binding of Tamiflu to Influenza N1 Neuraminidases. PLoS Comput. Biol. 2010, 6, e1000939. (29) Saglam, A. S.; Chong, L. T. Highly Efficient Computation of the Basal kon using Direct Simulation of Proteinâ−Protein Association with Flexible Molecular Models. J. Phys. Chem. B 2016, 120, 117−122. (30) Van Der Vries, E.; Collins, P. J.; Vachieri, S. G.; Xiong, X.; Liu, J.; Walker, P. A.; Haire, L. F.; Hay, A. J.; Schutten, M.; Osterhaus, A. D. H1N1 2009 pandemic influenza virus: resistance of the I223R neuraminidase mutant explained by kinetic and structural analysis. PLoS Pathog. 2012, 8, e1002914. (31) Rabenstein, B.; Knapp, E.-W. Calculated pH-dependent population and protonation of carbon-monoxy-myoglobin conformers. Biophys. J. 2001, 80, 1141−1150. (32) Kieseritzky, G.; Knapp, E.-W. Optimizing pKa computation in proteins with pH adapted conformations. Proteins: Struct., Funct., Genet. 2008, 71, 1335−1348. (33) Salomon-Ferrer, R.; Götz, A. W.; Poole, D.; Le Grand, S.; Walker, R. C. Routine microsecond molecular dynamics simulations with AMBER on GPUs. 2. Explicit solvent particle mesh Ewald. J. Chem. Theory Comput. 2013, 9, 3878−3888. (34) Case, D.; Babin, V.; Berryman, J.; Betz, R.; Cai, Q.; Cerutti, D.; Cheatham, T., III; Darden, T.; Duke, R.; Gohlke, H.; Goetz, A.; Gusarov, S.; Homeyer, N.; Janowski, P.; Kaus, J.; Kolossváry, I.; Kovalenko, A.; Lee, T.; LeGrand, S.; Luchko, T.; Luo, R.; Madej, B.; Merz, K.; Paesani, F.; Roe, D.; Roitberg, A.; Sagui, C.; Salomon-Ferrer, R.; Seabra, G.; Simmerling, C.; Smith, W.; Swails, J.; Walker, R.; Wang, J.; Wolf, R.; Wu, X.; Kollman, P. AMBER 14, 2014. (35) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 1977, 23, 327−341.

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank the DFG (German Research Foundation) grant SFB (Collaborative Research Center) 1035 (project B02) for financial support and Center of Integrated Protein Science Munich (CIPSM) for financial support. Computer resources for this project have been provided by the Gauss Centre for Supercomputing/Leibniz Supercomputing Centre under grant pr48po.

■

REFERENCES

(1) Copeland, R. A.; Pompliano, D. L.; Meek, T. D. Drug-target residence time and its implications for lead optimization. Nat. Rev. Drug Discovery 2006, 5, 730−739. (2) Pan, A. C.; Borhani, D. W.; Dror, R. O.; Shaw, D. E. Molecular determinants of drug−receptor binding kinetics. Drug Discovery Today 2013, 18, 667−673. (3) Luitz, M. P.; Zacharias, M. Protein−Ligand Docking Using Hamiltonian Replica Exchange Simulations with Soft Core Potentials. J. Chem. Inf. Model. 2014, 54, 1669−1675. (4) Fujitani, H.; Tanida, Y.; Ito, M.; Jayachandran, G.; Snow, C. D.; Shirts, M. R.; Sorin, E. J.; Pande, V. S. Direct calculation of the binding free energies of FKBP ligands. J. Chem. Phys. 2005, 123, 084108. (5) Dror, R. O.; Pan, A. C.; Arlow, D. H.; Borhani, D. W.; Maragakis, P.; Shan, Y.; Xu, H.; Shaw, D. E. Pathway and mechanism of drug binding to G-protein-coupled receptors. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 13118−13123. (6) Dror, R. O.; Green, H. F.; Valant, C.; Borhani, D. W.; Valcourt, J. R.; Pan, A. C.; Arlow, D. H.; Canals, M.; Lane, J. R.; Rahmani, R.; Baell, J. B.; Sexton, P. M.; Christopoulos, A.; Shaw, D. E. Structural basis for modulation of a G-protein-coupled receptor by allosteric drugs. Nature 2013, 503, 295−299. (7) Buch, I.; Giorgino, T.; De Fabritiis, G. Complete reconstruction of an enzyme-inhibitor binding process by molecular dynamics simulations. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 10184−10189. (8) Plattner, N.; Noe, F. Protein conformational plasticity and complex ligand-binding kinetics explored by atomistic simulations and Markov models. Nat. Commun. 2015, 6, 7653. (9) Ermak, D. L.; McCammon, J. A. Brownian dynamics with hydrodynamic interactions. J. Chem. Phys. 1978, 69, 1352−1360. (10) Zhou, H. X. Kinetics of diffusion-influenced reactions studied by Brownian dynamics. J. Phys. Chem. 1990, 94, 8794−8800. (11) Zhou, H.-X. Brownian dynamics study of the influences of electrostatic interaction and diffusion on protein-protein association kinetics. Biophys. J. 1993, 64, 1711. (12) Qin, S.; Pang, X.; Zhou, H.-X. Automated prediction of protein association rate constants. Structure 2011, 19, 1744−1751. (13) Northrup, S. H.; Allison, S. A.; McCammon, J. A. Brownian dynamics simulation of diffusion-influenced bimolecular reactions. J. Chem. Phys. 1984, 80, 1517−1524. (14) Yu, X.; Martinez, M.; Gable, A. L.; Fuller, J. C.; Bruce, N. J.; Richter, S.; Wade, R. C. webSDA: A web server to simulate macromolecular diffusional association. Nucleic Acids Res. 2015, 43, 220−224. (15) Wade, R. C.; Luty, B. A.; Demchuk, E.; Madura, J. D.; Davis, M. E.; Briggs, J. M.; McCammon, J. A. Simulation of enzyme−substrate encounter with gated active sites. Nat. Struct. Mol. Biol. 1994, 1, 65− 69. (16) Spaar, A.; Dammer, C.; Gabdoulline, R. R.; Wade, R. C.; Helms, V. Diffusional encounter of barnase and barstar. Biophys. J. 2006, 90, 1913−1924. (17) Gabdoulline, R. R.; Wade, R. C. Protein-protein association: investigation of factors influencing association rates by Brownian dynamics simulations. J. Mol. Biol. 2001, 306, 1139−1155. 5104


Article

Journal of Chemical Theory and Computation (36) Hopkins, C. W.; Le Grand, S.; Walker, R. C.; Roitberg, A. E. Long-Time-Step Molecular Dynamics through Hydrogen Mass Repartitioning. J. Chem. Theory Comput. 2015, 11, 1864−1874. (37) Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved side-chain torsion potentials for the Amber ff99SB protein force field. Proteins: Struct., Funct., Genet. 2010, 78, 1950−1958. (38) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general amber force field. J. Comput. Chem. 2004, 25, 1157−1174. (39) Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A. Automatic atom type and bond type perception in molecular mechanical calculations. J. Mol. Graphics Modell. 2006, 25, 247−260. (40) Bradbrook, G. M.; Gleichmann, T.; Harrop, S. J.; Habash, J.; Raftery, J.; Kalb, J.; Yariv, J.; Hillier, I. H.; Helliwell, J. R. X-Ray and molecular dynamics studies of concanavalin-A glucoside and mannoside complexes Relating structure to thermodynamics of binding. J. Chem. Soc., Faraday Trans. 1998, 94, 1603−1611. (41) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926−935. (42) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An N Log (N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089−10092. (43) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81, 3684−3690. (44) Andersen, H. C. Molecular dynamics simulations at constant pressure and/or temperature. J. Chem. Phys. 1980, 72, 2384−2393. (45) Hess, B.; Kutzner, C.; Van Der Spoel, D.; Lindahl, E. GROMACS 4: algorithms for highly efficient, load-balanced, and scalable molecular simulation. J. Chem. Theory Comput. 2008, 4, 435− 447. (46) Daura, X.; Gademann, K.; Jaun, B.; Seebach, D.; van Gunsteren, W. F.; Mark, A. E. Peptide folding: when simulation meets experiment. Angew. Chem., Int. Ed. 1999, 38, 236−240. (47) Ermak, D. L.; McCammon, J. Brownian dynamics with hydrodynamic interactions. J. Chem. Phys. 1978, 69, 1352−1360. (48) Baker, N. A.; Sept, D.; Joseph, S.; Holst, M. J.; McCammon, J. A. Electrostatics of Nanosystems: Application to Microtubules and the Ribosome. Proc. Natl. Acad. Sci. U. S. A. 2001, 98, 10037−10041. (49) Amaro, R. E.; Minh, D. D.; Cheng, L. S.; Lindstrom, W. M.; Olson, A. J.; Lin, J.-H.; Li, W. W.; McCammon, J. A. Remarkable loop flexibility in avian influenza N1 and its implications for antiviral drug design. J. Am. Chem. Soc. 2007, 129, 7764−7765. (50) Kunsch, H. R. The Jackknife and the Bootstrap for General Stationary Observations. Ann. Statist. 1989, 17, 1217−1241. (51) Klebe, G. Applying thermodynamic profiling in lead finding and optimization. Nat. Rev. Drug Discovery 2015, 14, 95−110. (52) Ladbury, J. E. Just add water! The effect of water on the specificity of protein-ligand binding sites and its potential application to drug design. Chem. Biol. 1996, 3, 973−980. (53) Woo, H.-J.; Roux, B. Calculation of absolute protein−ligand binding free energy from computer simulations. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 6825−6830. (54) Lee, M. S.; Olson, M. A. Calculation of absolute protein-ligand binding affinity using path and endpoint approaches. Biophys. J. 2006, 90, 864−877.

5105


Multiscale Simulation of ReceptorâDrug Association Kinetics

Recommend Documents

Multiscale Simulation of ReceptorâDrug Association Kinetics