Oxime Unbinding from AChE Gorge - ACS Publications - American

Subscriber access provided by CMU Libraries - http://library.cmich.edu

Article

Protein-Drug Interactions with Effective Polarization in Polarizable Water: Oxime Unbinding from AChE Gorge Arup K Pathak, and Tusar Bandyopadhyay J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.5b08930 • Publication Date (Web): 15 Oct 2015 Downloaded from http://pubs.acs.org on October 19, 2015

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B 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 33

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

The Journal of Physical Chemistry

Protein-Drug Interactions with Effective Polarization in Polarizable Water: Oxime Unbinding from AChE Gorge Arup K. Pathak and Tusar Bandyopadhyay∗ Theoretical Chemistry Section, Bhabha Atomic Research Centre, Mumbai 400 085, INDIA E-mail: [email protected]

Phone: +91 22 2559 0300. Fax: +91 22 2550 5151

(Preprint submitted to J. Phys. Chem. B Monday 12th October, 2015)

∗

To whom correspondence should be addressed

1 ACS Paragon Plus Environment


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

Abstract Despite the fact that polarizability of water is different in the bulk and in protein, simulations of protein-ligand complexes are mostly carried out in nonpolarizable water media. We present oxime (HI-6) unbinding from the active site gorge of AChE, known to be strongly influenced by intermolecular cation-π, hydrogen bridge (HB) and water bridge (WB) interactions, by molecular simulations with effective polarization in polarizable mean-field model of TIP3P water. Enabled by the recent availability of method of obtaining microkinetic of rare events, we set out to investigate the rate constants of unbinding transitions from one basin to the other through a combination of metadynamics and hyperdynamics simulations. The results underpin the importance of electronic polarization effects on the pathways, potential of mean force, rate constants, HB and WB dynamics of unbinding transitions of a drug molecule ligated to protein interior. The method is also applicable to unravel the binding mechanisms.

Keywords: Polarizable mean-field model of water, Effective polarization, Protein-drug complexes, Cation-π interaction, Metadynamics, Hyperdynamics, Kinetics, Hydrogen bond and water bridge dynamics


Page 2 of 33

Page 3 of 33

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


1. INTRODUCTION Cation-π interaction, 1 water bridge (WB) 2 and hydrogen bridge (HB) 3 formations are ubiquitous in biological recognition. Interactions of an ionic ligand with protein molecule, rich in aromatic residues, is one such example. Electronic polarizability plays a crucial role in deciding the fate of these recognition mechanisms. Molecular level predictions of the pathways and their associated microscopic rate constants during unbinding/binding transitions of protein-ligand complexes, 4 while these recognition features are at play, are essential for rational drug design and lead optimization. Molecular simulations could be a standard tool in this endeavour. However, current molecular mechanical models adopt two-body additive terms and thus are unable to treat electronic polarizability explicitly. Whereas, nonadditive force fields in polarizable models increases the computational cost significantly; moreover, consistency in the treatment of intra- and intermolecular polarizations is still an unresolved issue. 5 Thus to achieve a reasonable level of configurational sampling (that nonpolarizable force fields are good at) and to take account of the electronic polarizability (that polarizable force fields are built for), one has to look for optimal economical way to include polarizable molecular interactions in equivalent nonpolarizable models 6,7 for large biological simulations, like the present AChE (acetylcholinesterase)-oxime complex in explicit water that comprised of nearly 0.1 million atoms. Serine hydrolase, AChE is a pivotal enzyme that terminates nerve impulses by hydrolysing the cationic neurotransmitter acetylcholine at cholinergic synapses in a nearly diffusion controlled limit. The AChE monomer has an ellipsoidal shape with dimensions of ∼45 Å×60 Å×60 Å with a deep and narrow gorge; ∼20 Å long perforating almost half way into the centre of the monomer. 8 The base of the gorge holds the catalytic active site (CAS), while peripheral anionic site (PAS) 9 comprising of aromatic residues Tyr 72, Tyr 124, Trp 286, and Tyr 341 and an acidic residue Asp 74 is located 20 Å away from it at the rim of the active gorge that binds cationic ligands. The region between the CAS and PAS forms the hydrophobic pocket for binding aryl substrates and active site ligands. 3 ACS Paragon Plus Environment


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

Figure 1. (a) The protein-drug interaction system in a relaxed conformation of the bound complex. The cation-π interacting residues and the residues forming water bridge/hydrogen bond with the drug molecule (shown in red) are highlighted in grey and green colour, respectively. (b) Geometry and orientation used to describe cation-π interaction in the protein drug complex. The distance vector, roxime − rAromatic from one of the nitrogen atom of HI-6 to the geometrical centre of the aromatic ring and the angle, θ normal to the plane of the aromatic ring with the distance vector are shown. a and b are two vectors pointing towards two adjacent carbon atoms of the aromatic ring starting from its geometric centre. d and h are the projections of the distance vector onto the aromatic ring and along the normal to the aromatic plane, respectively. The two pyridinium rings of the oxime, one directed toward the catalytic active site (CAS pyridinium, CASP) and the other toward the peripheral anionic site (PAS pyridinium, PASP) of mAChE gorge, are also indicated.


Page 4 of 33

Page 5 of 33

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


The catalytic triad (serine, glutamic acid, and histidine) of the CAS is a vulnerable target of the lethal organophosphorus (OP) compounds. 10 Prior to aging of the enzyme due to OP compounds, the inhibited enzyme can be reactivated by strong nucleophiles, such as oximes. The high-resolution three dimensional structure of AChE, 11 complexed with the oximes provided valuable insight into the structural aspects involved in the interactions in the bound state: the crystallographic structure of bispyridinium oxime, HI-6 in complex with mAChE have revealed that one of its two pyridinium rings binds reversibly with peripheral site (PAS pyridinium; PASP ring) and the other is directed towards the active catalytic site (CAS pyridinium; CASP ring), as shown in Figure 1. The 4-carboxylamide-pyridinium ring (the PAS pyridinium ring) of HI-6 has been found to be sandwiched by the indole ring of Trp 286 and the phenol ring of Tyr 124 driven by strong cation-π interactions between the PASP ring and the side chains of Trp 286 and Tyr 124. 11 While the PAS of AChE holds this pyridinium ring tightly, the central linker of HI-6 is directed into the active-site gorge such that the 2-hydroxy-iminomethylpyridinium ring (CASP ring) could directed towards the active catalytic site. In addition, several residues form water bridge (WB) and hydrogen bond (HB) with the drug molecule. 12 Not all oximes involving all the aforementioned interactions with AChE are equally potent to liberate the inhibited enzymes. In fact, strategy of oxime drug administration differs from one country to the other. So far the major thrust has been to design oximes with highest binding affinity. However, the mean lifetime of the protein-drug complex, as decided by the microkinetics of unbinding, plays a decisive role for sustained drug efficacy. 13 Thus, fundamental studies of a realistic model, elucidating the structure-kinetic relationship of the protein-oxime complex in explicit water is of paramount importance. Here we study the unbinding transition of the oxime from the active-site gorge of AChE by conducting molecular dynamics simulations in electronic continuum (MDEC) model 6,7 that has been successfully applied in a number of occasions. In continuum electrostatics, the system is well defined by average charge distribution. In



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

effective polarizable model, in order to implicitly bring in the electronic screening effects, charges of all ionized groups and ions should be scaled by the inverse of the electronic part √ of water dielectric constant (1/ el , el = 1.78 ). 6,14 There are rigorous proofs on what conditions polarizable Hamiltonians are reduced to nonpolarizable ones and how this scaling factor appears naturally. On-the-fly, these charges can even be locally adjusted self consistently to reflect the local environment (e.g. the protein interior), such as the case in mean-field polarizable (MFP) model of water. 6 In effect, the resulting effective polarizable model, implicitly includes the electronic screening effects into an otherwise nonpolarizable model system and hence remains computationally as efficient as the standard nonpolarizable force field calculations are. Aided by this physically sound and numerically efficient method that have allowed us to study the protein-drug interactions with effective polarization in polarizable water medium, we have further employed a combination of metadynamics (MtD) and hyperdynamics (HD) simulations 15 to unravel the unbinding transition in a greater detail. We demonstrate that the strength of cation-π, HB and WB interactions, the unbinding transition pathways, the potential of mean force (PMF), and their associated rates can remarkably differ from one another depending on the treatment of electronic screening effects in the exemplary case of AChE-HI6 (oxime) complex.

2. COMPUTATIONAL DETAILS 2.1. Cation-π Interaction The subtle issues of the quantification of the strength of cation-π interaction are well known. 16 Here we endeavour to map the said interaction qualitatively, 17 through geometry based criteria. As shown in Figure 1 (b), the vector roxime − rAromatic , was defined as the distance between either of the cationic nitrogen atoms of HI-6 and the geometrical centre of 6 ACS Paragon Plus Environment

Page 6 of 33

Page 7 of 33

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


the aromatic ring of protein residue. And the angle θ accounts for the inclination of the said distance vector with the normal to the plane of the aromatic ring to be given as, ⎤ ⎡ a × b · (roxime − rAromatic ) ⎦. θ = Cos−1 ⎣ a × b |roxime − rAromatic |

(1)

Here the two vectors, a and b connects the centroid of the aromatic ring with its two adjacent carbon atoms. To account for the electrostatic attraction of the cationic charge toward the πelectron quadruple of the aromatic ring, an angular [f (θ)] and a spatial [f (r)] Fermi functions were considered: f (θ) =

1 1 + exp |θ|−0.785 0.785

;

(2)

and f (θ) =

1 Aromatic |−1 1 + exp |roxime −r0.6

.

(3)

The metacoordinate s, signifying the cation-π interactions between all the aromatic residues and the cationic nitrogen atoms of the ligand is given as, s=

nAromatic

f (θi ) × f (ri ).

(4)

i=1

In the bound state of the complex, we have monitored all possible cation-π interaction functionals between the two pyridinium rings of the oxime with the aromatic rings of the gorge at CAS, PAS, and at its inner wall by projecting the distance vector between the geometrical centre of the aromatic ring and the pyridinium nitrogen atom onto and along the normal to the aromatic plane, denoted by d and h, respectively.

2.2. System Setup We solvated the bound protein-drug complex (PDB ID code 2GYU) in a cubical box of 100 Å×100 Å×100 Å dimensions. The solvated system is comprised of a total 109646 atoms. 7 ACS Paragon Plus Environment


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

We have used two sets of water models: one, the primitive nonpolarizable TIP3P 18 and the other is the MFP/TIP3P. 6 In the MFP model of water, on-the-fly, the charges on the water atoms are locally adjusted self consistently to reflect the local environment. The mean-field nature of the local polarization is achieved through a well established damping scheme, parameterized by a averaging time τ . The resolution of the MFP model is limited by the HB lifetime (which is typically ∼1 ps). To achieve accuracy while averaging in the MFP model, several HBs should be broken and formed in between two simulation steps. Following Ref. 6, we have set a 5 ps value for τ . Due to HB and WB formations, water diffusion in protein environment is significantly reduced compared to that in bulk. Hence water molecules diffusing out during the relaxation time can be thought to be placed in a similar polar environment. This value of τ has been vouched to be used in biological applications of MFP/TIP3P water. All simulations were performed using GROMCAS 4.0.7 19 with MFP/TIP3P water model implemented. The protein was simulated using the AMBER99 force field. 20 HI-6 was modelled using the general AMBER force field (GAFF). 21 Their atomic charges were calculated with GAMESS suite of ab initio program system 22 by employing B3LYP/6-311++G** level of theory. One of the goals of the present work is to qualitatively capture the cation-π interaction in the protein-ligand system in MFP/TIP3P water. In all calculations we have maintained a charge neutral AChE, the active-site gorge of which is rich with aromatic residues having π electron cloud. These π electron clouds undergo cation-π interactions with the cationic nitrogen atoms of HI-6 as shown in Figure 1. In order to include effective electronic polarization in a mean-field way, the charges on the cationic nitrogen atoms of the oxime are rescaled √ (reduced) by the inverse of the electronic part of water dielectric constant (1/ el , el = 1.78 ). To compensate for the reduced charges on the two pyridinium N+ atoms, the excess charge is delocalized equally on the 5 carbon atoms of the two rings, without disturbing the charge distributions in the protein. The present treatment can be criticized for missing the protein polarization effects. However, as shown in reference 14 that even with this simplification, the


Page 8 of 33

Page 9 of 33

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


model already captures the major changes in polarizable water properties. The structure of HI-6 is presented in Supporting Information (SI) Figure S1 and its atomic charges are provided in Table S1. We energy-minimized the system before equilibrating the solvent in two stages; first, for 2 ns with the protein heavy atoms position restrained and then for 5 ns with restraints on the atoms of the protein backbone. Subsequently, in each of the solvents all the restraints are removed and 10 ns NPT equilibrium runs are performed at 298 K and 1 atm pressure using Nose-Hover thermostat 23 (with a coupling constant, 0.1 ps) and isotropic pressure coupling. The periodic boundary condition and minimum image convention were used along with the LINCS 24 procedure to constrain all bond lengths involving hydrogen. Long range electrostatics were treated by particle-mesh Ewald electrostatics. 25 The equilibrated structures were then used as the starting configurations for the subsequent MtD simulations.

2.3. Well-Tempered Metadynamics To accelerate the unbinding, which is otherwise too slow to be accomplished by equilibrium molecular dynamics (MD), one can adopt the well known metadynamics methods. 26 The original MtD formalism is based on building a history-dependent biasing potential, constructed from relatively small Gaussian-shaped repulsive hills, which are applied to a small number of specifically designed descriptors of the model, to be known as collective variables (CV). This discourages the system revisiting the same configurational space, resulting into faster sampling of configurational space. To study the unbinding transitions we have used protein-ligand centre of mass (COM) separation as CV. Because of its inherent nonequilibrium characteristics, this method is however, devoid of true equilibrium information unless the "slow" build-up limit is adopted. The well-tempered version of MtD 27 has solved this problem theoretically by rescaling the Gaussian height with the bias accumulated over time at a fictitious higher temperature, T + ΔT . One can control the extent of configurational space exploration by suitably tuning ΔT parameter between ΔT → 0 for equilibrium MD 9 ACS Paragon Plus Environment


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

and ΔT → ∞ for standard Mtd. Here we have used well-tempered version of metadynamics in the NVT-ensemble with a bias factor=10, Gaussian hills height=0.20 kJ/mol, width=0.05 nm and initial deposition rate of the hills=5ps. MtD simulations are performed utilizing PLUMED 28 plugin to the GROMACS 19 suite. At the end of the well-tempered runs, the converged free energy profiles of the unbinding transitions along the chosen COM CV are computed directly, 27 and through a reweighting scheme 29 when other unbiased degrees of freedom are used as a descriptor (such as cation-π metacoordinate) of unbinding.

2.4. Hyperdynamics Hyperdynamics simulations 30 were initiated from the bottom of each of the basins as discovered through the converged MtD runs. The accumulated hill potential along the chosen COM separation CV during a MtD run upto a certain build-up time is used so as to ensure that no bias is added at the barrier region (cf. Figure 2). Then we have performed equilibrium MD simulation using the physical potential as offered by the force field parameters, increased by the truncated time-independent MtD biasing potential. The last step allows ergodicity to be satisfied, since the metastabilties along the unbinding pathway are removed by the elevated potential and the diffusive dynamics in the chosen CV space can be sampled with ease with a reasonable computational effort. The biasing potential was read from external file using PLUMED 28 plugin. For the two water media, we have employed an ensemble of 50 independent starting configurations in each of the two basins with different initial velocities. The HD runs were continued for 15 ns, totalling a simulation time of 1.5 μs. Total number of escape events in the forward unbinding process during this simulation time are noted. These escape events also include all the forward trajectories that has once recrossed the barrier from reverse direction. The average boost factor, facc is calculated from the non-zero bias potential, Vbias (r) at a given basin in the bound state and calculated as, 10 ACS Paragon Plus Environment

Page 10 of 33

Page 11 of 33

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


ntot

1 facc = eβVbias (r) = eβVbias [r(ti )] , ntot i=1

(5)

where r can be read as COM separation and ntot is the number of such HD steps during which the bound states are confined to their respective basins.

3. RESULTS AND DISCUSSIONS 3.1. Unbinding Free Energy The aromatic residues lining the gorge are very much conserved among species. Molecular dynamics (MD) studies have shown that structural fluctuations of these residues are important in ligand binding. 31 The ligand unbinding through the gorge can thus be rare and the system can often rests for long in local minima before it crosses a barrier where old interactions are weakened and newer interactions start building with significant entropy changes. Conventional MD, limited by the sampling time scale that it can achieve, becomes impractical to study such rare events. MtD is a powerful method to enhance the sampling of the rare events. 26 The well-tempered variant 27,32 of which can successfully avoid the convergence and accuracy issues associated with the nonequilibrium simulation protcol. The inherent feature of the MtD method that no a priori assumptions are made on the final state, has also made it a powerful path finding tool in the multidimensional free energy landscapes. 33 In order to generate a pathway between the crystallographic bind state and the unbound state of the oximes, we first attempted performing MtD run with a single CV, namely the COM separation between the protein and the oxime.



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

Figure 2. PMF of the unbinding transition employing COM separation between the protein and the oxime as the CV in well-tempered MtD simulation for (a) nonpolarizable and (b) polarizable model. The three basins are denoted by numbers coloured in blue, while the two barriers are marked by red coloured numbers. The bias potentials (in kJ/mol) for HD runs initiated from basins 1 and 2 are indicated by dashed lines.

Figure 2 shows the PMF profile of the unbinding transitions with various basins and barriers along the pathway as marked by numbers. Snapshots of structures corresponding to these basins and barriers are provided Figure 3 and Figure 4 for the two water media. Clearly free energy barriers associated with the inter-basin transitions and the transition pathways are markedly different depending on the inclusion of the water polarization effects. This result also demonstrate the distinctive role played by water molecules in the unbinding process. Note that outside the gorge at the unbound state, the orientations of HI-6 with respect to the protein are different: in nonpolarizable model calculation, HI-6 lie flat with the outer surface of the protein, while in polarizable model calculation it lies somewhat perpendicular to the surface of the protein. In both the cases the CAS pyridinium ring of HI-6 still interacts with the residues in the rim of the gorge leading to basin 3. This results into higher COM separation (≈ 1 nm) in the case of polarizable model calculation to arrive at an equivalent of basin 3 of nonpolarizable situation. 12 ACS Paragon Plus Environment

Page 12 of 33

Page 13 of 33

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


Figure 3. Snapshots of the energetic basins and barriers, as marked in Figure 2, for the unbinding transition of HI-6 (shown in red color) in polarizable water medium are shown. Only those residues of AChE binding gorge that stabilizes these states and are responsible for drug release microkinetics are depicted. Neither the water molecules forming WBs, nor any HBs and cation-π interactions between AChE and HI-6 are shown for the sake of clarity. See below for further details.



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

Figure 4. Same as Figure 3 but in nonpolarizable medium.


Page 14 of 33

Page 15 of 33

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


Figure 5. The well-tempered metadynamics free energy landscape. 2D contour plot of the unbinding transitions of the oxime drug from the binding gorge of the enzyme in (a) polarizable and (b) nonpolarizable water media are shown, as obtained from reweighting the configurations of the one-dimensional welltempered metadynamics run, employing protein-oxime COM separation as CV. Isoenergy lines are drawn every 5 kJ/mol as indicated by colour bar. While reweighting, the sum of cation-π interactions between CAS pyridinium N: with the aromatic rings of W86, Y337, and Y341 and PAS pyridinium N: with the aromatic ring of W286 are considered.

By a reweighting algorithm 29 for recovering equilibrium Boltzmann distribution of any unbiased variable from a well-tempered metadynamics simulation, we have obtained the twodimensional free energy landscape of the unbinding transition. The results are presented in Figure 5. The plots are prepared from a monodimensional metadynamics in COM separation CV, which was then reweighted to obtain the cation-π metacoordinate (s) distribution. The employed metacoordinate is the sum of cation-π interactions prevailing in the bound state of the complex (cf. Figures S2 and S3). As can be seen in Figure S2, the cation-π order parameter in the equlibrium bound state has a strong dependence on the angle of approach to the centre of the aromatic ring, which are further modified depending on the polarizability of water model used. In particular, the difference in the plots of W86-CASP and Y337-CASP interactions are noteworthy. The lower dipole moment of water in the protein interior causes a strong cation-π order parameter to prevail in the polarizable water model. This is also reflected in 2D PMF plot of unbinding transition (Figure 5). For example, in contrast to nonpolarizable water model, two metastable states along the gorge unbinding pathway could be 15 ACS Paragon Plus Environment


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 16 of 33

clearly located in polarizable water. This result clearly shows the dependence of protein-drug interactions and hence the oxime’s unbinding pathway on the treatment of water polarization.

Table 1. Exponential fitting parameters for the protein-oxime water bridge lifetime correlation functions at various basins and barriers of the PMF profile during the oxime unbinding transition in polarizable water medium. The unbinding transition is performed employing centre of mass separation between AChE and HI-6 as CV during the metadynamics run. The major protein residues participating in the water bridge formations are indicated.

a

Basin/ Barrier

Protein Residue

Time Constants (ps)

Amplitude (%)

Lifetime (ps)

Basin 1

GLU285 GLU292 SER298

1974.2 4.8 ; 137.6 2285.1

100 55 ; 45 100

-a 71.72 494.66

Barrier 1

ASP283 GLU285 GLU292 SER298

1.6 ; 33.6 190.7 3.2 ; 29.3 53.9 ; 1.6

44 ; 56 100 63 ; 37 34 ; 66

24.56 109.21 15.68 46.75

Basin 2

ASP74 ASP283 GLU285 SER298

3.4 ; 132.8 30.1 ; 3.7 25.7 ; 504.8 47.9 ; 1.8

30 53 48 23

70 47 52 77

12.21 18.39 187.95 20.66

Barrier 2

ASP74

1.7 ; 31.6

87 ; 13

13.23

Basin 3

ASP74

8.2

100

8.46

; ; ; ;

insufficient simulation data

3.2. Hydration Dynamics During Unbinding Transition In bulk, in the protein surface, and in the protein interior, such as in the active site gorge of AChE, water manifests different dipole moments. In effect, water molecules buried deep inside a hydrophobic protein binding pocket result into a decreased dielectric constant and subsequent stabilization of the electrostatic interaction, such as HB prevailing between protein and drug molecule. Consequently, rupture of this type of HBs involves an energetically penalized transition state (HBs can be up to 1.2 kcal/mol stronger in hydrophobic environments 34 ) and thus can dictate the lifetime and binding kinetics of a protein-drug complex. Water molecules facilitate the unbinding 35 by inserting into the broken protein-drug HB and forming WB using two of its hydrogen atoms to bond with the liberated protein residue 16 ACS Paragon Plus Environment

Page 17 of 33

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


and the drug. The presence of water thus decreases the energy barriers that need to be surmounted to rupture the HBs. Water has been peddled as the "lubricant" that acts as an expeditor of the protein conformational fluctuations involved in the structural relaxation. 36 Table 2. Same as Table 1 but in nonpolarizable water medium.

a

Basin/ Barrier

Protein Residue

Time Constants (ps)

Amplitude (%)

Lifetime (ps)

Basin 1

GLU285 SER293 SER298 TYR124

-a 1.9 ; 70.2 -a 298.3 ; 31.3

-a 70 ; 30 -a 54 ; 46

-a 28.59 -a 168.51

Barrier 1

GLU285 SER293 SER298 TYR124

74.1 ; 3524.5 537.6 2.3 ; 72.5 2.5 ; 65.6

41 ; 59 100 64 ; 36 56 ; 44

2126.7 -a 84.65 54.75

Basin 2

HIP287 SER293 TRP286 TYR124

2.6 ; 215.3 0.03 ; 11.3 133.3 ; 0.2 31.7 ; 2.1

54 50 38 30

; ; ; ;

46 50 62 70

112.62 15.59 61.59 57.53

Barrier 2

HIP287 SER293 TRP286 TYR72 TYR124

2.3 ; 79.6 62.5 ; 2.3 170.5 ; 0.3 312.6 ; 28.0 230 ; 10

78 14 26 37 25

; ; ; ; ;

22 86 74 63 75

26.76 13.1 55.65 83.5 19.93

Basin 3

HIP287 SER298 TYR72

2.3 187.6 ; 76.4 6.2

100 50 ; 50 100

3.54 134.85 12.85


In view of this, we have identified the residues participating in HB and WB formation along the oxime unbinding pathways by conducting equilibrium MD simulations of the structures at various basins and barriers along the gorge. Each of the NVT simulations are run for 10 ns, the last 8 ns of which are analyzed. The coordinates of the solvated system are saved every 1 ps. The structures at the barriers were subjected to a harmonic constraining potential to hold the protein-ligand system at the barriers’ top. Then we have investigated the microscopic aspects of hydration dynamics by evaluating the WB and HB lifetime correlation functions. 17 ACS Paragon Plus Environment


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

Figure 6. Time dependence of protein-oxime water bridge lifetime correlation function, C WB (t) in (a) polarizable and (b) nonpolarizable water medium. Major residues of AChE responsible for the water bridge formation at various basins and barriers along the oxime’s unbinding pathway are shown. Their lifetime are reported in Table 1 and Table 2, respectively.

To elucidate the hydration dynamics we use a geometric criterion for hydrogen bonding, according to which a hydrogen bond donor (D) and acceptor (A) are considered hydrogen bonded if the D-A distance is less than 3.5 Å, and the D-H-A angle is less than 35◦ . The 18 ACS Paragon Plus Environment

Page 18 of 33

Page 19 of 33

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


WB network relaxation time is defined in terms of the decay of the hydrogen bonding bridges

Table 3. Exponential fitting parameters for the correlation functions of the mean lifetime of proteinoxime hydrogen bond in (a) polarizable and (b) nonpolarizable water media at various basins and barriers along the unbinding pathways. The unbinding transitions are performed employing centre of mass separation between AChE and HI-6 as CV during the metadynamics run. The protein residues participating in the hydrogen bond formations are indicated. Basin/ Barrier

Protein Residues

Time Constants (ps)

Amplitude (%)

a) Polarizable Water Basin 1

ASP74, GLU285, GLU292, SER298, TYR77, TYR124, TYR337

5350.7

100

Barrier 1

ARG296, ASP283, GLU285, HIP287, PRO78, SER298, TRP286, TYR72, TYR77

153.8 ; 3262.2

25 ; 75

Basin 2

ASP74, ASP283, GLU285, TRP286, TYR72, TYR124

223.7 ; 864.2

21 ; 79

Barrier 2

ASP74, GLU292, HIP 287, PRO290, SER293, TRP286, TYR72, TYR124

9207. 4 ; 458.8

15 ; 85

Basin 3

ASP74, TYR72, TYR77

201.6 ; 702.2

45 ; 55

b) Non-polarizable Water

a

Basin 1

GLU285, SER298, TYR124, TYR337

1641

100

Barrier 1

ASP74, GLU285, SER293, SER298, THR83, TRP86, TRP286, TYR72, TYR124

-a

-a

Basin 2

ARG296, HIP287, TRP286, TYR124, TYR341

720.1

100

Barrier 2

ASP74, ASP283, HIP287, PRO78, SER293, SER298, TRP286, TYR72, TYR77, TYR124, TYR341

2271.9

100

Basin 3

ASP283, GLU292, HIP287, THR75

256.8 ; 9.2

19 ; 81


correlation function, C WB (t) = h(0)h(t) / h , where h(t) is a WB population operator, which is equal to one when the two associated HBs of the bridge are intact at time t , but it becomes zero when either of the HB breaks, and the angular brackets denote an average over all bridge forming water molecules. Similar to this, the protein-drug HB lifetime dynamics 19 ACS Paragon Plus Environment


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

can be measured by the correlation function, C HB (t) =

ij

Page 20 of 33

Sij (t0 )Sij (t + t0 )/

ij

Sij (t0 ),

where like h(t), Sij (t) is a HB population operator. The values of Sij specifies HB between all D-A pairs involving all the participating residues and the oxime polar atoms at time t. The functions C WB (t) and C HB (t) is the probability that a random pair that is water bridged or hydrogen bonded at time zero is still bonded at time t, regardless of whether or not the bond was broken at intermediate times. Thus, beyond an initial transient period, their decay reflects the rearrangement of the underlying network.

Figure 7. Time dependence of protein-oxime hydrogen bond lifetime correlation function, C HB (t) in (a) polarizable and (b) nonpolarizable water medium. Hydrogen bonds between all possible donor-acceptor pairs are considered. Major residues of AChE responsible for the hydrogen bond formation at various basins and barriers along the oxime’s unbinding pathway are reported in Table 3.

We found that both the lifetime correlation functions are in general multiexponential, with the average lifetime depending significantly on the specificity. For the two water models, list of residues participating in the hydration dynamics, the fitting parameters of the 20 ACS Paragon Plus Environment

Page 21 of 33

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


correlation functions, their amplitude and average lifetime are provided in Tables 1-3 and Figures 6 and 7. Subtlety of the inclusion of polarizable water model are clearly manifested in the hydration dynamics, which eventually dictates the unbinding pathways. In the bound state, the oxime can form strong HB with the residues lining the gorge wall. The hydrophobic environment in the gorge has been found to be correctly reflected in the polarizable water model; the HB lifetime correlation function in all the three binding basins decays much more slowly compared to its nonpolarizable version. The intrusion of water into protein-oxime HB region, breaks up the old HB and form WB in place of it. The nonpolarizable model treats the water molecules without any distinction in protein and in bulk. In effect, the WB formed by these molecules are seen to follow long lifetimes when compared with its polarizable counterpart in the hydrophobic binding pockets of the protein (for example, see Table 1 and 2 for the E285 and S298 WB data that can be used for direct comparison in basin1 and barrier 1). In summary, treating electronic polarizability in the unbinding simulation ensures that water trapped inside the protein binding pocket results into weaker WB, which in effect causes stronger protein-drug HB, both being facilitated by its lower dipole moment in the protein interior.

3.3. Cation-π Interaction During Unbinding Transition The existence of the noncovalent force due to cation-π interactions in shaping up the proteindrug unbinding kinetics is irrefutable. 1,37 The interactions of cationic ligands with the PAS and the CAS region of AChE, rich in aromatic residues, is a paradigmatic case to cite. 33,38 The acceptance of this interaction in structural biology, due to the works of Dougherty 37 has led to the currently accepted electrostatic model of the said interaction, in which a cation is steered electrostatically by the permanent quadruple moment offered by the π system. The bispyridinium oximes can simultaneously form two cation-π interactions through the cationic



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

Figure 8. Scatter plot of the cation-π interaction functional in the (d, h) plane during oxime’s unbinding. The residues interacting with the CASP:N and PASP:N of the oxime are as mentioned. Their values (s) are represented by colour code as indicated by spectral scale. The data accumulated over last 10 ns of 12 ns equilibrium production runs in each of the basins and barriers along the pathway, for polarizable (left panel) and nonpolarizable (right panel) water media are shown. In all the plots, (0,0) represents the six-membered aromatic ring centroid of the tryptophan and tyrosine residues. Star symbol on the line, h = 0 indicates equidistance to any carbon atom of six membered aromatic ring, while open circle represents distance to the remotest carbon atom of the five membered ring in the tryptophan residue.

nitrogen atoms of the CAS pyridinium (CASP:N) and the PAS pyridinium (PASP:N) rings. During the unbinding transition, we have explored all the possibilities of cation-π interactions between the two cationic nitrogens of HI-6 with W286, W86, Y72, Y77, Y124,


Page 22 of 33

Page 23 of 33

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


Y337, Y341, F295, F297, F338, and F346 and found that only tryptophan and tyrosine residues forms significant cation-π interaction. 39 Figure 8 shows the scatter plots of the most significant cation-π functionals (cf. Equation 4) during the oxime’s sojourn from the bound to the unbound state through the intervening barriers. In order to study the cation-π interaction that is electrostatic in nature, the inclusion of electronic screening effects through the effective polarizable model is certainly more realistic, which is shown here to steer the oxime’s unbinding differently when compared with its nonpolarizable counterpart. Namely, the decreased dipole moment of MFP water in the protein interior stabilizes the electrostatic interaction and hence the cation-π order parameters are in general seen to be in higher side along the gorge pathway. The cation-π interaction and hydration dynamics together show that the microenvironment through which ligand unbinding takes place is largely dictated by the influence of water on the protein dynamical transition.

3.4. Microkinetics of Unbinding Transitions With the demonstration of the effects of the inclusion of polarizability in place, we set out to determine the rate constants of the rare unbinding events by a combined Mtd and HD method. 15 In this combined approach, a single MtD run is used to discover a pre-optimized biasing potential (cf. Figure 2) that is subsequently "frozen" and is used as the nonnegative biasing potential, Vbias to accelerate the event in subsequent HD runs. 30 The dynamics in the auxiliary HD system becomes faster and allows one to map the unbiased state-to-state time evolution by a suitable renormalization of time through a boost factor, facc (cf. Equation 5). By conducting MtD and HD separately one can ensure that the bias potential vanish at the separatrix between one basin and the neighbouring states in the forward unbinding process, so that the transition state theory (TST) is applicable. Thus the rigorous posteriori assessment, as is required in the case of on-the-fly combination of both steps, 4 of fulfilling this requirement can be avoided. We have initiated several HD runs at different starting configurations from the bottom of each of the basins. The forward microkinetic rate is de23 ACS Paragon Plus Environment


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 24 of 33

fined as the fraction of reactants (say, the population in basin 1) that turns into products (say, the population in basin 2) per unit of time. The COM separation CV in MtD-HD run can be regarded as a reaction coordinate, ξ that separates the free energy wells (cf. Figure 2) through the barrier at ξ # . The transition state at ξ = ξ # divide the planes between two basins and are sparsely populated areas at the top of the barriers. The hydration dynamics 40 or for that matter, the cation-π interactions in the dividing plane can be thought to play a crucial role on the fate of this population and hence on the measureable rates of unbinding. Therefore the hyperdynamics rate that implicitly assumes each protein-drug configurations crossing the transition state will end up in the product basin (and hence contribute to the rate), needs to be corrected by the transmission coefficient in order to take into account the possiibility of recrossings at the barrier top. From the ensemble of 50 HD runs, we have identified the total number of escape trajectories in the forward direction (as reported in Table 4 below). These trajectories include configurations crossing, ξ # with a positive velocity rapidly recrosses the barrier before settling in the product well. The biased rate constants, kbias are calculated by constructing histograms of the first-passage-time (FPT) of the barrier crossing. 15 The decay profiles of the time dependent FPT follows a first order rate equation as can be seen in Figure 9. kbias is finally equated with the decay constant of the FPT profile. The HD rates are then found from kbias /facc . To compensate for the solvent influenced recrossing events at the barrier top, the HD rate is multiplied with transmission coefficient, κ. 40 From Onsager’s regression hypothesis it follows that 41

# ˙ # δ ξ(0) − ξ ξ(0)θ ξ(t) − ξ ,

κ(t) = ˙ ˙ ξ(0) δ [ξ(0) − ξ # ] ξ(0)θ


(6)

Page 25 of 33

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


Figure 9. Evaluation of HD simulated biased unbinding rate constants, kbias from basin 1 to basin 2 [ (a) and (c) ] and from basin 2 to basin 3 [ (b) and (d) ] in polarizable (upper panel) and nonpolarizable (lower panel) water medium. In each of the plots the histograms of the FPT, deduced from hyperdynamics runs and its single exponential decay patterns (in red colour) are shown. The time dependent transmission coefficients, κ(t) as obtained from long equilibrium simulations at the top of each barriers are shown in the insets. The plateau values of κ(t = ∞) as mentioned in the insets are used to find the rates.

where θ is the Heaviside step function. For time longer than the typical time of molecular motions but much shorter than the time constant of the reaction, κ will stabilize at a plateau value. At this point of time all the reactants that crossed the barrier at time 0 will be found either in the product well or in the reactant well. The numerator of Equation 6 will then contain contributions only from net reactive flux to the product basin and therefore will tend to reach a constant at long time:

κ = lim κ(t). t→∞


(7)


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 26 of 33

In order to calculate the transmission coefficient, the structures at the barriers’ top were used. For each of the configurations at the barrier top, additionally a positive forward velocity is assigned through uniform random number. The system is then followed forward and backward in time starting from each of these configurations. The trajectories so produced can be referred to as "Bennett-Chandler" trajectories. 41 Direct numerical implementation of Equation 6 suffers from efficiency problems. 42 For instance to record a transmission coefficient of 0.1 one would need to generate ∼ 104 trajectories in order to get a 10% relative error. In view of this, we began by generating 10 to 100 protein-drug configurations constrained at the ξ = ξ # dividing plane. Then we have started many fleeting trajectories for 100 ps to 2 ns from this constrained ensemble, initiated with randomized velocities drawn from a MaxwellBoltzmann distribution. Four different types trajectories are identified and their respective number of occurrences are counted: (1) NR→P , for trajectories that had originated from reactant well and end up in product well, (2) NR→R , for reactant to reactant trajectories, (3) NP →R , for product to reactant trajectories, (4) NP →P , for product to product trajectories. Clearly, NR→P makes a positive contribution to κ, NP →R a negative contribution, while NR→R and NP →P are unproductive events. An alternative to the numerical implementation of Equation 6 could thus be the straightforward evaluation of the ratio: κ(t) =

NR→P (t) − NP →R (t) . NR→P (t) + NP →R (t) + NR→R (t) + NP →P (t)

(8)

The insets in Figure 9 show the estimates of κ(t). Since in the limit t → 0+ all trajectories are confined to the reactant well, κ(0) = 1. Recrossings make κ smaller, and hence 0 ≤ κ ≤ 1. Finally the rate constant is given as,

k = κ × (kbias /facc ) .

(9)

In short, the MtD-HD simulations ensures that the bias potential is not added at the barriers’ top, separating any two basins, such that the rates (kbias /facc ) evaluated through this method can be regarded as TST rates. 15 However, the collective character of the solvent 26 ACS Paragon Plus Environment

Page 27 of 33

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


motion and the protein-drug cation-π interactions at barriers’ top can cause some of the forward-flux reactive configurations to recross the barrier into the reactant well. This can result into the actual forward rate constants markedly deviate from the TST rates. 40 Transmission coefficient (cf. Equation 8), evaluated from unbiased equilibrium simulations of an ensemble of configurations at the barriers’ top is a safety parameter that takes care of the recrossing events. Finally, the actual unbinding rate is given as the product of transmission coefficient and the TST rate (cf. Equation 9). Table 4. Evaluation of rate constants of unbinding transitions in a) polarizable and b) nonpolarizable water media Basin i → Basin j

Number of escape events

1→2 2→3

1289 150

1→2 2→3

5075 5476

k bias (s-1 )

Boost Factor, f acc

a) Polarizable water 47 2.83×108 1.85×109 1697 b) Nonpolarizable water 33 2.13×109 1.65×109 28

κ(t = ∞)

k (s-1 )

0.047 0.11

2.8×105 1.2×105

0.072 0.123

4.6×106 7.2×106

Table 4 lists the number of escape events and the biased rate constants as obtained through the HD runs, the boost factor and the long-time limit of the transmission coefficient for the unbinding transitions in the two water media. The robustness 15 of the MtD-HD combined approach is further checked by conducting similar sets of HD simulations taking truncated bias potentials’ locations (cf. Figure 2 dashed lines) at two other positions (Supporting Information, S4). We found that the rates are within ± 15 % limit of the values reported here in. In all these studies, the unbinding rate constants are found to be about ten times lower (and hence their associated mean FPT are larger by an equal magnitude) in polarizable water in comparison to nonpolarizable medium. This is entirely in line with our expectations that polarizable water molecules inside a protein binding pocket can provide additional stability towards protein-drug HB and cation-π interactions, which can act as a possible kinetic traps in the course of unbinding. To the best of our knowledge till date there is no report on the experimental value of the microkinetic rates of the studied system. It 27 ACS Paragon Plus Environment


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

would be interesting to compare our theoretical findings with experimental data whenever they are available .

4. CONCLUSIONS The major new aspects of the present theoretical study may be highlighted as follows. Classical molecular simulations with effective polarization in mean-field polarizable model of TIP3P water is employed to thoroughly investigate the oxime (HI-6) unbinding from the active site gorge of AChE, known to be strongly influenced by intermolecular cation-π, HB and WB interactions. We do observe major changes in the HB and WB network and also in the cation-π interactions with the introduction of polarizable mean-field model of water in comparison to the nonpolarizable water. Potential of mean force is calculated by employing well-tempered metadynamics for the unbinding of HI-6 from AChE both in polarizable and nonolarizable TIP3P water. It is clearly observed that the free energy barriers associated with the inter-basin transitions and the transition pathways are markedly different depending on the inclusion of the water polarization effects. 2D contour plots of the unbinding process are also generated by using a suitable reweighting method, where cation-π metacoordinate is used as one of the collective variable. Two metastable states along the gorge unbinding pathway could be clearly located in polarizable water, in contrast to the one in nonpolarizable model. WB and HB lifetime correlation functions are in general found to be multiexponential, with the average lifetime depending significantly on the specificity of the solvent and the WB/HB interacting partners. The WB are seen to follow long lifetimes in nonpolarizable water medium when compared with its polarizable counterpart in the hydrophobic binding pockets of the protein. In polarizable water, the HB lifetime correlation function in all the three binding basins decays much more slowly compared to its nonpolarizable version. It is also observed that the decreased dipole moment of polarizable


Page 28 of 33

Page 29 of 33

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


water in the protein interior stabilizes the electrostatic interaction and hence the cation-π order parameters are in general seen to be in higher side along the gorge pathway of AChE. Rate of the unbinding process are also calculated from a combination of MtD and HD runs. The unbinding rate constants are found to be about ten times lower in polarizable water in comparison to nonpolarizable medium. Thus, we have convincingly shown the importance of the inclusion of water polarization in the molecular dynamics of protein-ligand complex in an otherwise effective polarizing environment that only takes note of ionic charge rescaling of the substrates. We expect that the present MtD-HD-based combined MD approach with the effective polarization provides a physically justifiable and efficient route to arrive at a structure-kinetic relationships, which can serve as an indispensible tool for lead optimization towards drug discovery. We believe that the present work is the first comprehensive treatment of microkinetic rates of protein-ligand unbinding transitions, dictated by the three major biological recognition mechanisms at play in polarizable aqueous medium.

Supporting Information Available The supporting material contains the chemical struture, atomic charges of HI-6, scatter plot of the cation-π interaction functional in the bound state of the AChE•HI-6 complex model system, and convergence of well-tempered metadynamics runs. The material is available free of charge via the Internet at http://pubs.acs.org.

Acknowledgement The authors thank computer centre, BARC for providing the ANUPAM parallel computational facility. One of us (TB) sincerely appreciates Prof. Igor Leontyev for several fruitful discussion and generously giving access to the GROMACS code, implementing MFP-TIP3P water model. The work was supported by DAE under project XII-N-R&D-02.04/Theoretical & Computational Chemistry of Complex Systems. The authors declare no conflict of interest.



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

References (1) Meyer, E. A.; Castellano, R. K.; Diederich, F. Interactions with Aromatic Rings in Chemical and Biological Recognition. Angew. Chem. Int. Ed. 2003, 42, 1210-1250. (2) Pal, S. K.; Zewail, A. H. Dynamics of Water in Biological Recognition. Chem. Rev. 2004, 104, 20992123. (3) Desiraju, G. R.; Steiner, T. The Weak Hydrogen Bond in Structural Chemistry and Biology. 1999 , Oxford University Press, New York. (4) Tiwary, P.; Limongelli, V.; Salvalaglio, M.; Parrinello, M. Kinetics of Protein-Ligand Unbinding: Predicting Pathways, Rates, and Rate-Limiting Steps. Proc. Natl. Acad. Sc. (USA) 2015, 112, E386-E391. (5) Ponder, P. J.; Case, D. A. Force Fields for Protein Simulations. Adv. Prot. Chem. 2003, 66, 27-85. (6) Leontyev, I. V.; Stuchebrukhov, A. A. Polarizable Mean-Field Model of Water for Biological Simulations with AMBER and CHARMM Force Fields. J. Chem. Theory Comput. 2012, 8, 3207-3216. (7) Leontyev, I. V.; Stuchebrukhov, A. A. Polarizable Molecular Interactions in Condensed Phase and their Equivalent Nonpolarizable Models. J. Chem. Phys. 2014, 141, 014103 (12 pages). (8) Quinn, D. M. Acetylcholinesterase: Enzyme Structure, Reaction Dynamics, and Virtual Transition States. Chem. Rev. 1987, 87, 955-979. (9) Bourne, Y.; Taylor, P.; Radic, Z.; Marchot, P. Structural Insights into Ligand Interactions at the Acetylcholinesterase Peripheral Anionic Site. EMBO J. 2003, 22, 1-12. (10) Newmark, J. The Birth of Nerve Agent Warfare: Lessons from Syed Abbas Foroutan. Neurology 2004, 62, 1590-1596. (11) Ekstrom, F.; Pang, Y-P.; Boman, M.; Artursson, E.; Akfur, C.; Borjegren, S. Crystal Structures of Aacetylcholinesterase in Complex with HI-6, Ortho-7 and Obidoxime: Structural Basis for Differences in the Ability to Reactivate Tabun Conjugates. Biochem Pharm 2006, 72, 597-607. (12) Kesharwani, M. K.; Ganguly, B.; Das, A.; Bandyopadhyay T. Differential Binding of Bispyridinium Oxime Drugs with Acetylcholinesterase. Acta Pharma. Sinicia 2010, 31, 313-328. (13) Copeland, R. A.; Pompliano, D. L.; Meek, T. D. Drug Target Residence Time and its Implications for Lead Optimization. Nat Rev Drug Discov. 2006, 5, 730-739. (14) Kohagen, M.; Lepsik, M.; Jungwirth, P. Calcium Binding to Calmodulin by Molecular Dynamics with Effective Polarization. J. Phys. Chem. Lett. 2014, 5, 3964-3969. (15) Schneider, J.; Reuter, Karsten. Efficient Calculation of Microscopic Dissolution Rate Constants: The Aspirin-Water Interface. J. Phys. Chem. Lett. 2014, 5, 3859-3862. (16) Minoux, H.; Chipot, C. Cation-π Interactions in Proteins: Can Simple Models Provide an Accurate Description? J. Am. Chem. Soc. 1999, 121, 10366-10372. (17) Branduardi, D.; Gervasio, F.L.; Cavalli, A.; Recanatini, M.; Parrinello M. The Role of the Peripheral Anionic Site and Cation-π Interactions in the Ligand Penetration of the Human AChE Gorge. J. Am. Chem. Soc. 2005, 127, 9147-9155. (18) Jorgensen, W. L.; Jenson, C. Temperature Dependence of TIP3P, SPC, and TIP4P Water From NPT Monte Carlo Simulations: Seeking Temperatures of Maximum Density. J. Comput. Chem. 1998, 19, 1179-1186.


Page 30 of 33

Page 31 of 33

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


(19) 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. (20) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.;, Roitberg, A.; Simmerling, C. Comparison of Multiple AMBER Force Fields and Development of Improved Protein Backbone Parameters. Proteins 2006, 65, 712-725. (21) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollamn, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157-1174. (22) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S. et al. General Atomic and Molecular Electronic Structure System. J. Copmut. Chem. 1993, 14, 1347-1363. (23) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31, 1695-1697. (24) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E M. LINCS: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18, 1463-1472. (25) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577-8593. (26) Barducci, A.; Bonomi, M.; Parrinello, M. Metadynamics. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2011, 1,826-843. (27) Barducci, A.; Bussi, G.; Parrinello, M. Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method. Phys. Rev. Lett. 2008, 100, 020603-020606. (28) Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R. A.; Parrinello, M. PLUMED: A Portable Plugin for Free-Energy Calculations with Molecular Dynamics. Comp. Phys. Comm. 2009, 180, 1961-1972. (29) Bonomi, M.; Barducci, A.; Parrinello, M. Reconstructing the Equilibrium Boltzmann Distribution from Well-Tempered Metadynamics. J. Comp. Chem. 2009, 30, 1615-1621. (30) Voter, A. F. A Method for Accelerating the Molecular Dynamics Simulation of Infrequent Events. J. Chem. Phys. 1997, 106, 4665-4677. (31) Shen, T.; Tai, K.; Henchman, R. H.; McCammon, J. A. Molecular Dynamics of Acetylcholinesterase. Acc. Chem. Res. 2002, 35, 332-340. (32) Barducci, A.; Bonomi, M.; Parrinello, M. Linking Well-Tempered Metadynamics Simulations with Experiments. Biophys. J. 2010, 98, L44-L46. (33) Pathak, A. K.; Bandyopadhyay, T. Unbinding Free energy of Acetylcholinesterase Bound Oxime Drugs Along the Gorge Pathway from Metadynamics-Umbrella Sampling Investigation. Proteins 2014, 82, 1799-1818. (34) Gao, J.; Bosco, D. A.; Powers, E. T.; Kelly, J. W. Localized Thermodynamic Coupling Between Hydrogen Bonding and Microenvironment Polarity Substantially Stabilizes Proteins. Nat. Struct. Mol. Biol. 2009, 16, 684-690. (35) Schmidtke, P.; Luque, F. J.; Murray, J. B.; Barril, X. Shielded Hydrogen Bonds as Structural Determinants of Binding Kinetics: Application in Drug Design. J. Am. Chem. Soc. 2011, 133, 18903-18910.



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

(36) Sinha, V.; Ganguly, B.; Bandyopadhyay, T. Energetics of Ortho-7 (Oxime Drug) Translocation through the Active-Site Gorge of Tabun Conjugated Acetylcholinesterase. Plos One 2012, 7, e40188 (14 pages). (37) Dougherty, D. A. Cation-π Interactions in Chemistry and Biology: A New View of Benzene, Phe, Tyr, and Trp. Science 1996, 271, 163-168. (38) Branduardi, D., Gervasio, F. L.; Cavalli, A.; Recanatini, M.; Parrinello, M. The Role of the Peripheral Anionic site and Cation-π Interactions in the Ligand Penetration of the Human AChE Gorge. J. Am. Chem. Soc. 2005, 127, 9147-9155. (39) Pathak, A. K.; Bandyopadhyay, T. Ortho-7 Bound to the Active-Site Gorge of Free and OP-Conjugated Acetylcholinesterase: Cation-π Interactions. Biopolymers 2015, DOI: 10.1002/bip.22712 (in press). (40) Ciccotti, G.; Ferrario, M.; Hynes, J. T. Dynamics of Ion Pair Interconversion in a Polar Solvent. J. Chem. Phys. 1990, 93, 7137-7147. (41) Chandler, D. Statistical Mechanics of Isomerization Dynamics in Liquids and the Transition State Approximation. J. Chem. Phys. 1978, 68, 2959-2970. (42) Ruiz-Montero, M. J.; Frenkel, D.; Brey, J. J. Efficient Schemes to Compute Diffusive Barrier Crossing Rates. Mol. Phys. 1997, 90, 925-942.


Page 32 of 33

Page 33 of 33

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


Table of Contents (TOC) Graphic:


Oxime Unbinding from AChE Gorge - ACS Publications - American

Recommend Documents