Dynamic Mechanism of a Fluorinated Oxime Reactivator Unbinding

Subscriber access provided by - Access paid by the | UCSB Libraries

B: Biophysical Chemistry and Biomolecules

Dynamic Mechanism of a Fluorinated Oxime Reactivator Unbinding from AChE Gorge in Polarizable Water Arup Kumar Pathak, and Tusar Bandyopadhyay J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.8b01171 • Publication Date (Web): 19 Mar 2018 Downloaded from http://pubs.acs.org on March 19, 2018

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

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 44 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

Dynamic Mechanism of a Fluorinated Oxime Reactivator Unbinding from AChE Gorge in Polarizable Water Arup K. Pathak∗ and Tusar Bandyopadhyay∗ Theoretical Chemistry Section, Bhabha Atomic Research Centre, Mumbai 400 085, INDIA E-mail: [email protected]; [email protected]

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

(Preprint submitted to J. Phys. Chem. B Friday 16th March, 2018)

∗

To whom correspondence should be addressed

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 A well-tempered metadynamics simulation is performed to study the unbinding process of a fluorinated oxime (FHI-6) drug from the active site gorge of acetylcholinesterase enzyme in a polarizable water medium. Cation-π interactions, water bridge, and hydrogen bridge formations between the protein and the drug molecule are found to strongly influence the unbinding process, forming basins and barriers along the gorge pathway. Distinct unbinding pathways are found when compared with its recently reported nonfluorinated analogue, HI-6. For example, due to their permanent positive charges on both the pyridinium rings of HI-6, it exhibits the minimum in the potential of mean force of the unbinding process in the gorge mouth (where the peripheral anion site, PAS of the enzyme is located), largely caused by cation-π interaction. Whereas, the same interaction, both in the catalytic active-site (CAS) and PAS region is found to be greatly enhanced in its lipophilic fluorinated analogue, FHI-6, causing a deep potential energy minimum in the bound state. This may render FHI-6 to be held more firmly in the CAS region of the gorge, as is also evidenced from the microkinetics of unbinding transitions, measured through a combination of metadynamics and hyperdynamics simulations.

2


Page 2 of 44

Page 3 of 44 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 Acetylcholinesterase (AChE) is a serine protease enzyme that hydrolyzes the neurotransmitter acetylcholine in order to terminate the nerve impulses. It catalyze the breakdown of neurotransmitter acetylcholine into acetic acid and choline in neuromuscular junctions and cholinergic brain synapses. 1,2 The high turnover number and catalytic activity of the primary cholinesterase enzyme, AChE is limited by the diffusion of the target. A catalytic triad consists of a serine, a glutamate and a histidine residue is responsible for its high catalytic activity. 2–5 The monomeric form of AChE has a ellipsoidal shape with a deep and narrow gorge of 20 Å length. 2,6–8 The catalytic active site (CAS) is situated at the base of the narrow gorge. For musmusculus AChE (mAChE) the catalytic triad is comprised of S203, E334 and H447 residues and an anionic subsite is comprised of W86, E202 and Y337 residues. The peripheral anionic site (PAS) comprises of many aromatic residues (Y72, Y124, W286 and Y341) and D74 residues which sits 20 Å away from the CAS at the rim of the active gorge. The region between the PAS and CAS forms a hydrophobic pocket and helps to bind ligand and aryl substrates. The substrate is captured at the rim of the gorge and then is migrated through the hydrophobic pocket to the catalytic site. It is quite known that hydrogen bond (HB), water bridge (WB) formations and cation-π interactions are omnipresent in many biological recognition processes. 9–11 These interactions also play essential role in oxime drug capture and binding/unbinding in the gorge of the protein. 12,13 The drug binding/unbinding pathways at molecular level and their corresponding microscopic rate constants are important for rational drug design and lead optimization. 13,14 The catalytic triad of the AChE enzyme is a vulnerable target of many organophosphorus (OP) compounds. 15 OP compounds are known as notorious poison that leads to severe public health hazard with large number of fatalities annually worldwide, 16 and are considered as looming threats to be used as chemical warfare agents in terrorist attack and armed conflicts. The small lipophilic OP compounds not only inhibit the catalytic activity of AChE at the neuromuscular junction but also in the central nervous system (CNS) by penetrating the blood brain barrier (BBB). 17–19 The BBB, although a highly selective permeable barrier with very high electrical resistivity that allows the passage of water and the diffusion of small

3



lipid-soluble molecules (e.g. O2 , CO2 , hormones etc.), glucose and amino acids, which are crucial to normal neural function, is also permeable to small lipophilic OP agents. These OP compounds form a covalent conjugate by phosphorylation with the serine residue of the catalytic triad of AChE and thus hindered it’s normal catalytic activity in CNS and also in neuromuscular junction. The hindrance of the catalytic activity leads to the death due to asphyxiation if note treated promptly. Prior to aging (by dealkylation or deamination of the phosphorus conjugate) of the OP inhibited AChE, the inhibited enzyme can be reactivated by strong nucleophiles, such as oximes. 20,21 Large number of oximes including mono-pyridinium and bis-pyridinium oximes are used to reactivate the OP inhibited AChE. 21 Out of them, bis-pyridinium oximes are most potent as they bind with the AChE reversibly by using both the pyidinium rings. One (PAS pyridinium; PASP) of its two pyridinium rings is used to bind with peripheral anionic site (PAS) and other one (CAS pyridinium; CASP ring) is engaged with the active catalytic site (CAS). The oxime of the CASP helps the reactivation of the inhibited enzyme. Many electrostatic interactions between oxime drug and AChE, namely, AChE-oxime cation-π , HB interactions and the WBs connecting the protein and the oxime drug molecules are responsible for binding the drug in the active-site gorge. However, effectiveness of all these oxime drugs is limited and one of their major drawbacks, as summarized by Mercey et. al., 21 is that their permanent positive charge prevent them from crossing BBB in order to reactivate inhibited AChE in CNS. As a result, the charged oximes that are very effective to reactivate the OP inhibited AChE in the periphery, their efficacy is limited in the CNS. 21,22 Hence, the neurological effects of OP poisoning, like seizures, convulsions and psychological effects to the infected individuals cannot be prevented by administering the existing oxime drugs. In view of this, many painfull strategies are developed that even include methods like direct injections into the brain. Introduction of uncharged reactivators e,g. RS41A, RS69N and RS194B containing five, six and seven membered rings 23–29 have witnessed little success. Recently it is revealed, 21 based on membrane permeability measurements, that the BBB crossing ability is increased in proportion to the number of fluorine atoms attached to a given drug molecule. This observation has opened up a new avenue of drug development that are designed to cross the BBB. Thus another effective way to reactivate the inhibited 4


Page 4 of 44

Page 5 of 44 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


AChE in the CNS would be the fluorination of existing oxime drugs. The introduction of fluorine atom into the heterocyclic ring of pyridinium oximes makes the oxime drugs more lipophilic and thus they can invade the BBB. So the reactivation chances of inhibited AChE in the CNS increase rapidly with the administration of fluorinated oxime drugs with enhanced lipophilicity. 30,31 However, the designed fluorinated oximes, apart from enhanced BBB permeability, must remain tightly bound to the active site gorge of AChE in order to orchestrate the enzyme’s reactivation. Fortunately, HI-6, a commonly used bis-pyridinium drug 20 when fluorinated (FHI-6) is found to be more effective than its parent compound (HI-6) for the reactivation of OP inhibited AChE. 21 The unbinding process of an ionic ligand, like FHI-6, in water solvent from the active-stie gorge of AChE, rich in aromatic residues, will certainly be dictated by cation-π, HB and WB interactions. Encouraged by our recent report 13 on the unbinding of its non-fluorinated counterpart (HI-6), here we set out to explore the unbinding process of FHI-6 from the active-site gorge of AChE with a primary aim to investigate the effect of fluorination on the CASP ring of drug molecule. Electronic polarizability plays a deciding role during the the course of the drug’s unbinding due to the above mentioned electrostatic interactions. In view of this, we have conducted all-atom molecular dynamics (MD) simulations with effective polarization in mean-field polarizable (MFP) model of TIP3P water. 32,33 Both one dimensional (1D) and two dimensional (2D) potential of mean force (PMF) are calculated for the drug’s unbinding by applying enhance sampling technique, namely, well-tempered version of the meta-dynamics (WT-MtD) simulation tool. 34 A combination of metadynamics and hyperdynamics simulations are then used on the WT-MtD discovered potential wells to calculate the micro-kinetics of unbinding transitions from one binding pocket to other. We have shown that the strength of cation-π, HB, WB interactions, the unbinding transition pathways, the PMF profiles and the micro-kinetics of unbinding transitions can remarkably differ due to single-atom fluorination of the oxime drug.

5



2. METHODOLOGY 2.1. System Setup The initial coordinates of the AChE•HI-6 complex is taken from protein data bank (PDB : 2GYU). 20 The HI-6 drug is removed from the complex before the drug-free protein is docked with FHI-6. Autodock is used to produce initial co-ordinates of the AChE•FHI-6 complex. 35 Different protein-drug (AChE•FHI-6) conformations are generated by employing Lamarkian genetic algorithm (LGA). Energy based autodock scoring function is used to find the possible binding site of the fluorinated oxime drug while the scoring function is comprised of electrostatic interaction, solvation and loss of entropy due to drug binding, short range van der Waals and hydrogen bonding interactions. Initially, the whole protein is considered for docking. The grid size of 130 × 130 × 130 along the x, y, z axes, respectively, with a grid spacing of 0.750 Å is used to find possible binding sites of FHI-6. After that, a more refined search is conducted by setting the grid size to 45 × 45 × 45 along the x, y, z axes with a grid spacing of 0.320 Å. The maximum number of energy evaluation is set to 2.5×106 in the LGA run. Root mean square (with 0.2 nm tolerance) clustering is carried out on the resulting conformers of the AChE•FHI-6 complex. The resulting conformers are then ranked according to their energy values. A relaxed conformation of the docked complex is shown in figure 1. The lowest energy conformer of the AChE•FHI-6 complex is then used as docked structure for subsequent MD simulations. The docked AChE•FHI-6 complex is solvated by employing MFP version of the three site water model, namely, MFP/TIP3P in a cubical box having dimension of 10.3 nm × 10.3 nm × 10.3 nm. In the MFP model, the charges on the atoms of the water are adjusted self consistently ‘on-the-fly’ to mimic the heterogeneous solvent environment. A well-known damping scheme is developed 32,33 to achieve the mean-field nature of the local polarization. The damping scheme is parameterized by the average time, τ , which is dictated by the HB life-time (≈ 1ps). The accuracy of the averaging in the MFP model is highly dependent on the HB making and breaking during two consecutive simulation steps. Therefore, in order to achieve a very good accuracy, several HBs must break and form during two simulation steps. τ =5 ps is used in the present study, since diffusion of water molecules in protein 6


Page 6 of 44

Page 7 of 44 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: The fluorinated analogue of HI-6 oxime drug, FHI-6 in complex with the protein, AChE. The bound complex is shown in a relaxed conformation. The aromatic residues of the protein involved in cation-π interactions with the FHI-6 drug are shown by pink colour surface representation. The FHI-6 drug is shown by violet colour licorice representation. The structures of FHI-6 oxime is also shown in the lower panel. The cationic site pyridinium ring and peripheral anionic site pyridinium ring are marked as CASP and PASP, respectively.

7



interior are restricted by the HB and WB formation. All the MD simulations are performed by employing GROMCAS 4.0.7 suite of program. 36 The protein, AChE was simulated using the AMBER99 force field parameters. 37 The fluorinated oxime drug, FHI-6 is modeled using the general AMBER force field (GAFF). 38 The atomic charges of the FHI-6 are calculated by using GAMESS 39 suite of ab initio program system at B3LYP/6-311++G** level of theory. 40 B3LYP hybrid density functional is comprised of Becke’s three parameter exchange and correlation due to Lee, Yang and Parr. The force field parameters of the drug can be obtained from the author on request. It is important to mention that a charge neutral AChE is maintained during the MD study by adding adequate number of sodium ions. The oxime drug has cationic site (nitrogen atom) in each of its pyridinium rings (see figure 1). The bound drug molecule experience strong cation-π interactions with the π-electron clouds of aromatic residues (tyrosine, phenyl alanine and tryptophan) lining the gorge wall (figure 1). In order to capture the protein-drug cation-π interactions properly that would reflect the electronic screening effect implicit in effective polarizable model, the charges of all ionized groups or atoms are rescaled (reduced) by the inverse of the electronic part of water dielectric constant i.e. by a factor of 0.75 √ (1/ el , el =1.78). 32,33 Hence the positive charges on the nitrogen atoms, N:CASP (0.3286 a.u.) and N:PASP (0.3656 a.u.) of the oxime drugs are rescaled to 0.2465 a.u. and 0.2742 a.u., respectively, for N:CASP and N:PASP. The excess charge is delocalized over the five carbon atoms in each pyridinium rings to compensate for the reduction of the charge on the cationic nitrogen atom of pyridinium ring due to the rescaling. The present model is shown to capture the major changes in binding characteristic in a polarizable water medium. 12,13 The charge neutral solvated system is finally consisted of approximately 0.11 million atoms. Note that in the protein interior, the aromatic residues in the binding pocket would typically correspond to a slightly higher el value of about 2.0. This would cause a change in the scaling factor by 0.04 (0.75-0.71), which in turn would cause a nominal change in difference (≈0.01 a.u.) of the scaled charges. In view of this and that there is rigorous proof 33 on what conditions polarizable Hamiltonians are reduced to nonpolarizable ones and how the scaling factor of 1.78 appears naturally, we have rescaled the charges by a factor of 0.75. Both steepest descent and conjugant gradient methods are used to energy minimize the 8


Page 8 of 44

Page 9 of 44 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


solvated system. The system is then slowly heated to a temperature of 298 K before the equilibration run in NPT ensemble. 10 ns equilibrium run is carried out at 298 K and 1 atm pressure. Nose-Hover thermostat 41,42 with a coupling constant of 0.1 ps is employed to maintain the temperature at 298K and the pressure is maintained by applying isotropic pressure coupling. The system is then subjected to a 20 ns NVT ensemble equilibrium run. All bond lengths involving hydrogen are constrained by linear constraint solver algorithm (LINCS). 43 Minimum image convention with periodic boundary conditions in all three directions is used for all the simulation. The long-range electrostatic interactions are modeled by employing Particle-mesh Ewald sum algorithm. 44 The well equilibrated solvated system of the protein-drug complex is then used for subsequent enhanced non-equilibrium simulation in order to capture the rare events efficiently within reasonable computational time.

2.2. Well-Tempered Metadynamics Due to the very complex nature, the protein-drug system may have large number of local minima (basins) and maxima (barriers) along the unbinding pathways, dictated by a large number of influencing parameters, namely, structural fluctuations, electrostatic interactions, solvation, van der Walls, HB and WB interactions. The system may spend large amount of time in a local potential minima before it escapes through a barrier to other basin. Hence, the protein-drug unbinding is a notoriously slow process. Non-equilibrium enhanced sampling technique, namely, metadynamics ( MtD) 45 can effectively be used to accelerate the unbinding process and to study such rare event in a reasonably time frame . The free-energy surface as a function of suitable descriptor of the system, namely, collective variable (CV) can be constructed by employing MtD simulation technique. To discourage the system revisiting the same configurational space that has already been sampled during the simulation, history dependent Gaussian shaped biasing potential is added to the system along a particular CV. However, due to the inherent non-equilibrium nature of the MtD technique, it is difficult to attain true equilibrium information, unless the "slow" build-up limit of Gaussian hill is considered. This well-known convergence problem of the MtD method is overcome by rescaling the Gaussian hill height (W) in well-tempered version of the MtD method (WT-MtD) as 34,45 9



W = ω 0 τG e

−

VG (CV,t) kB ∆T

.

Page 10 of 44

(1)

where, ω 0 is initial deposition rate, τ G is deposition stride, ∆T is a temperature and VG (CV, t) is the bias potential accumulated in CV, while k B is the Boltmann constant. In WT-MtD method, the initial deposition rate ω 0 of the bias potential decreases with the bias accumulated over time (t) along a CV. The CVs are sampled at a fictitious higher temperature T+∆T . ∆T →0 corresponds to equilibrium MD and ∆T → ∞ corresponds to the standard MtD. One can control the extent of free energy surface exploration by suitably tuning the ∆T parameter. Gaussain hills of width=0.05 nm and initial height=0.50 kJ/mol are added at 5 ps intervals with a bias factor of 10 during the WT-MtD simulation in the NVT ensemble. All the WT-MtD simulations are carried out by utilizing PLUMED 46 plugin to the GROMACS 36 suite of MD program. The centre-of-mass (COM) separation between the protein and oxime drug is used as a suitable CV to study the un-binding process. The free energy profiles are calculated directly at the end of the WT-MtD runs for the drug unbinding from the protein, AChE. The PMF profile is also calculated along other unbiased degrees of freedom, namely, cation-π meta-coordinate by employing an appropriate reweighting scheme. 47

2.3. Model for Cation-π Interaction In order to explore the unbinding PMF as a function of cation-π interaction 48–50 (by employing cation-π meta-coordinate as an unbiased CV during WT-MtD run) that the drug molecule experiences with the aromatic residues lining the gorge wall, we have defined a cation-π meta-coordinate [S(Oxime − Aromatic)] based on distance and angle between the aromatic residues and the drug molecule to be given as

S(Oxime − Aromatic) =

n X

fθ (θi ) × fr (~roxime − ~rAromatic ) ,

i=1

10


(2)

Page 11 of 44 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


where, the sum runs over the aromatic residues and the Fermi functions are expressed as,

fθ (θi ) =

1 1+

exp |θ|−0.785 0.785

;

(3)

and fr (~roxime − ~rAromatic ) =

1 1+

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

.

(4)

Here, ~roxime −~rAromatic is the distance vector and is defined as the distance between either of the cationic nitrogen atoms (N:CASP or N:PASP) of FHI-6 and the COM of aromatic residues (Y72, Y77, W86, Y124, W286, F295, F297, Y337, F338, Y341, F346) of the AChE. The inclination of the distance vector with the normal to the plane of the aromatic ring is measure by angle, θ, which can be calculated as,   ~ ~a × b · (~roxime − ~rAromatic ) , θ = Cos−1  ~ ~a × b |~roxime − ~rAromatic |

(5)

Where ~a and ~b are the two vectors which join the centroid of the aromatic ring with two of its adjacent carbon atoms.

2.4. Quantum Mechanical Calculations The Trp, Tyr, Phe residues in the rim and gorge are involved in cation-π interaction with the drug. To explore further, full geometry optimizations are carried out for the complexes of drug and the Trp, Tyr, Phe residues by employing B3LYP functional with split valance 6-31++G(d,p) Gaussian type basis functions. 40 Updated hessian based Newton-Raphson algorithm is used to find the most stable structure for each complex in the multidimensional potential energy surface. Electrostatic interaction energy between the oxime drug and aromatic residue for each of the complexes is calculated by employing localized molecular orbital energy decomposition analysis (LMOEDA) method. 51 GAMESS suit of program is used to perform all the quantum mechanical calculations. 39

11



Page 12 of 44

2.5. HB and WB Time Correlation Functions Hydration dynamics is measured through the protein-drug HB and WB time correlation functions. Based on geometric criterion, a donor (D) and an acceptor (A) are said to be hydrogen bonded when the distance between them is less than 3.5 Å and the D-H-A angle is less than 35◦ . For WB both the HBs of a water moleules connecting the protein and drug molecule must satisfy the above condition. The HB time correlation function, C HB (t) is defined as,

CHB (t) =

X hPmn (0)Pmn (t)i hPmn i

mn

,

(6)

where, Pmn (t) is the HB bond population operator. HB bonds between all D-A pairs involving polar atoms of FHI-6 and all the participating amino acid residues of the AChE at time t is described by Pmn (t). The HB bond population operator is equal to 1 when the particular HB is survived at time t and it becomes zero when the same HB bond does not survive. The angular bracket denote the average over all D-A pairs between AChE and drug in a particular basin or barrier. Similarly, WB time correlation function, C WB (t) is defined as,

CW B (t) =

hA(0)A(t)i , hAi

(7)

where, A(0) and A(t) are the WB population operator at time 0 and t, respectively. The WB population operator is equal to 1 when both the concerned HB are intact at a certain time and is equal to 0 when either of the associated HB does not survive. The average over

12


Page 13 of 44 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


all WB forming water molecules in a particular basin or barrier is denoted by the angular bracket.

2.6. Hyperdynamics Kinetics of protein- drug unbinding process is very important to study in order to ascertain the efficacy and suitability of a particular drug. The rate constant of protein-drug unbinding process can be calculated by employing hyperdynamics (HD) simulation. 13,52,53 The HD simulations are performed from the bottom of the WT-MtD discovered basin along the unbinding pathway. In order to ensure that no bias is added to the barrier region, the accumulated hills along the COM CV up to a certain time is used for the hyperdynamics run (see Fig. 2). The accumulated MtD biasing potential is then added to the physical potential as offered by the force-field parameters in order to carry out equlibrium MD runs. The elevated potential helps to remove meta-stabilities along the unbinding pathway and ergodicity is also satisfied. Due to this, the sampling of the diffusive dynamics along the chosen COM CV can be done with ease in a reasonable computational time. PLUMED plugin is used to read the biasing potential from an external file during the MD simulation. 60 independent starting configurations are used for the HD run and each HD run is carried for 25 ns. During the HD simulation time, the total number of escape events are counted. It is worthwhile to mention that all the forward trajectories that have once re-crossed the barrier from reverse direction are also included in that total escape events. The average boost factor, fac , due to the elevated potential can be calculated from the nonzero bias potential, Vb (R) in the bound state correspond to a particular basin and is expressed as,

n

βVb (R)

fac = e

1 X βVb [R(ti )] = e , n i=1

(8)

where R is the CV (COM), n is the total number of HD steps in which the bound state is 13



confined to its respective basin, and β is inverse temperature.

3. RESULTS AND DISCUSSION Here we first present the result for converged 1D unbinding free energy profile along the biased COM separation CV (subsection 3.1). The cation-π interactions between the aromatic residues lining the AChE gorge wall and the drug molecule is well known to shape up the oxime’s unbinding pathway and its kinetics. 13,50 In view of this, in subsection 3.2 we have parameterized the cation-π interaction based on geometry based criterion, and in subsection 3.3 its quantum mechanical considerations are discussed. Next the unbinding PMF in the COM separation–cation-π interaction space is presented (subsection 3.4). The HB and WB dynamics that plays crucial role in the unbinding process of the drug molecule is presented in subsection 3.5. Finally, in subsection 3.6 we present the kinetics of transition from one basin to the other of the 1D unbinding free-energy profile.

3.1. 1D Unbinding Free Energy Profile Before we embark upon the free-energy profiles, it is important to note that convergence of metadynamics method is a major issue. The well-tempered version 34 of the method has although solved this problem theoretically by adding progressively smaller bias as the simulation continues, the problem of estimating local convergence remains. 54 In order to ensure convergence of the free energy profile, at least in the discovered basins as shown in Figure 2, we have conducted additional analysis (see figures in Supporting Information, SI). They are (i) the heights of the Gaussian hills as a function of time (see SI, Fig. S1), which is seen to approach a value close to zero in a 20 ns run and is a signature of complete diffusive behaviour

14


Page 14 of 44

Page 15 of 44 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 profile for the unbinding transition by employing centre-of-mass separation between AChE and oxime drug, FHI-6 as collective variable is shown. The PMF profile is generated from the converged WT-MtD simulation. The basins and barriers along the unbinding path are marked by black and red numerals,respectively. The bias potential for HD run is indicated by dashed line. The locations of the corresponding potential minima for basins 1 and 2 of HI-6 unbinding 13 are indicated by vertical blue dotted lines.

of the system in the explored phase space; (ii) PMF profiles at different time intervals, which are aligned by setting the global minimum to zero (see SI, Fig. S2), and the PMF profile in the discovered basins is found to be unaltered after 10 ns till the metadynamics run is completed at 20 ns; (iii) time evolution of cation-π and COM separation CV space (see SI, Fig. S3), where the system is found to diffuse efficiently with no newer regions in the CV space can be seen to be explored after 10 ns; (iv) The free energy profile as obtained through reweighting 47 on the biased CV (COM separation) itself and that from sum- hills facility of PLUMED 46 plugin are compared (see SI, Fig. S4), where an agreement between them, at least in the bound states, indicates that the calculated profile is converged. The aromatic residues in the CAS and PAS site of the AChE are conserved among various species 55 and their structural fluctuations shapes up the cation-π interactions with the oxime drugs. During the drug unbinding old HB, WB and cation-π interactions must break and newer ones must form in place leading to many barrier and local minima along the drug passage pathway. The PMF profile for the unbinding transition of the FHI-6 drug is

15



Figure 3: Snapshots of the basins and barriers (marked in Fig. 2) along the unbinding pathway of FHI-6 drug (purple colour) are shown. Only few important residues responsible for the stabilization of a particular state are shown. The HB and WB forming residues are shown by blue and red colour CPK representations, respectively (see Table 2 and 3). Orange colour CPK representation is used for residues which form both HB and WB.

16


Page 16 of 44

Page 17 of 44 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


calculated from the converged WT-MtD run and shown in Fig. 2 is a testimony to this fact. When compared with the unbding tranistion of its non-fluorinated counterpart (HI-6), 13 the FHI-6 drug molecule shows distinctly different behavior. The unbinding of the FHI-6 consists of two basins and one barrier. On the other hand three basins and two barriers are observed in the unbinding pathway for the HI-6. The position of basin 1, appears at 1.38 nm and 1.51 nm along the protein-drug COM separation for FHI-6 and HI-6 drug, respectively. The most stable basin (deepest minima along the unbinding transition pathway), basin 2 in both the cases are also situated at different COM separation. For FHI-6, it appears at 1.96 nm and for HI-6 it is located at 2.26 nm along the COM separation coordinate. From this observation one can infer that both the basins (1 and 2) are situated at lesser COM separation for FHI-6 drug and that FHI-6 will penetrate more inside the gorge. Thus distinctive role of one-atom fluorine substitution in the HI-6 drug can clearly be established from the 1D PMF profiles. The snapshots of the structure of the various basins and barriers along the unbinding pathway are shown for FHI-6 drug in Figure 3. For clarity, only a few important residues are shown in each of the snapshots. It is worthwhile to mention that the unbinding pathways from the AChE gorge varies distinctly depending on the drug molecules. For example, anti-Alzheimer’s disease drug: Huperzine A,, 56,57 oxime drugs: OBI and Ortho-7, 12,58 showed vastly different binding/unbinding pathways along the AChE gorge.

3.2. Cation-π Interaction along the Unbinding Pathway The non-covalent cation-π interactions of the aromatic residues in the gorge and in the rim with the oxime drug play essential role in shaping up the oxime drug binding and unbinding kinetics. 13,50 A few tryptophan (Trp), tyrosine (Tyr) and phenylalanine (Phe) residues in the gorge and rim are engaged in the cation-π interactions with the cationic nitrogen atom of CAS (CAS:N) and PAS (PAS:N) pyridinium rings. 59,60 In the currently accepted electrostatic model for the cation-π non-covalent interaction due to Dougherty, it is considered that a cation is steered electrostatically by the permanent quadruple moment caused by the π electron cloud, like π electron of aromatic resides. 59 It is interesting to mention that the bis-pyridinium drug like, FHI-6 can form two cation-π interactions simultaneously 17



Figure 4: Scatter plots of the cation-π interaction functional in basins (BAS1 and BAS2) and barrirer (BAR1) in the (d; h) plane during FHI-6 unbinding. d and h are the projections of the distance vector (~roxime − ~rAromatic ) onto the aromatic ring and along the normal to the aromatic plane, respectively. The residues interacting with the CASP nitrogen and PASP nitrogen of the oxime are mentioned. Their values (S) are represented by colour code as indicated by spectral scale. The six-membered aromatic ring centroid of the tryptophan and tyrosine residues are represented by (0,0) in all the plots. 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.

18


Page 18 of 44

Page 19 of 44 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


due to the presence of two cationic nitrogen atoms (CAS:N and PAS:N). Here we explore all possible cation-π interactions by considering Y72, Y77, W86, Y124, W286, F295, F297, Y337, F338, Y341, F346 residues. To study the cation-π interactions, fresh equilibrium MD simulations for 20 ns are carried out at each basins (BAS1 and BAS2) and barrier (BAR1). The scatter plots of the most significant cation-π interaction functional during FHI-6 unbinding are shown in Fig. 4. It is worthwhile to mention that, cation-π interactions due to tyrosine, tryptophan and phenylalanine are significant for the FHI-6 drug. On the other hand, cation-π interaction due to phenylalanine is not significant in case of HI-6 drug. 13

Table 1. Calculated electrostatic interaction energy in kcal/mol based on localized molecular orbital energy decomposition analysis for minimum energy complexes of fluorinated HI-6 (FHI-6) and HI-6 drugs with various aromatic residues, namely, tyrosine (Tyr), tryptophan (Trp) and phenyl alanine (Phe) at B3LYP/6-31++G(d,p) level of theory.

Tyr-FHI-6 Tyr-HI-6 Trp-FHI-6 Trp-HI-6 Phe-FHI-6 Phe-HI-6 Electrostatic energy

-10.95

-10.04

-8.56

-7.61

-9.23

-8.73

3.3. Quantum Mechanical Calculations Electrostatic interaction energy between the oxime drug and aromatic protein residues also serves as a measure of cation-π interaction. Quantum mechanical calculations are employed to find out the electrostatic interaction energy between the drugs (HI-6 and FHI-6) and isolated aromatic residues, namely, Trp, Phe and Tyr. This will help us to categorize the Trp, Phe and Tyr residues (like Y72, Y77, W86, Y124, W286, F295, F297, Y337, F338, Y341, F346) responsible for cation-π interaction in the during the drug unbinding. Here, a total of six complexes of FHI-6 and HI-6 with the Tyr, Trp and Phe residues are considered. Full geometry optimizations are carried out for all the complexes by employing B3LYP functional with split valnce 6-31++G(d,p) bais set. 40 The electrostatic interaction energy are calculated for each of the complexes by employing LMOEDA 51 method and are reported in Table 1. 19



It is clear from Table 1, that FHI-6 interact with the aromatic residues more strongly than HI-6. This may be a probable reason behind the vast difference in the PMF profiles for their unbinding.

Figure 5: The 2D contour plot for the unbinding transition of the oxime drug, FHI-6 from the AChE gorge based on WT-MtD simulation. The 2D free energy landscape is generated from reweighting the configurations of the one-dimensional WT-MtD simulation, employing AChE-oxime COM separation as CV. Isoenergy lines are drawn every 10 kJ/mol. The sum of all possible cation-π interactions between the oxime drug and protein residues (Y72, Y77, W86, Y124, W286, F295, F297, Y337, F338, Y341, F346) are considered for reweighting. The numerals 1 and 2 indicate the the locations of two basins (c.f. Figure 2). For comparison the green star symbols indicate the locations of the two basins for HI-6 unbinding 13 in the same contour space.

3.4. 2D Potential of Mean Force for Unbinding As the cation-π interactions are very important to shape up the drug unbinding process from the gorge of the protein, it will be also very interesting to consider cation–π interactions as a CV in constructing two dimensional (2D) PMF profile. The PMF profile along any unbiased variable (cation–π CV) can be constructed from the knowledge of the biased CV ( COM CV) profile with the help of reweighting algorithm. 47 Here, sum of cation–π interactions with Y72, Y77, W86, Y124, W286, F295, F297, Y337, F338, Y341, F346 aromatic residues as a CV are considered to construct 2D PMF surface. The 2D contour plot is shown in Figure 5. The functional from of the cation–π metacoordinate (S) is provided in Eq. 2. The distributions 20


Page 20 of 44

Page 21 of 44 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


of the cation–π metacoordinate (S) are generated from the one dimensional WT-MtD run, biased along the COM separation CV. It is also interesting to mention that 2D PMF profile agrees well with the 1D PMF profile (c.f. Figure 2) which has two basins and a barrier along the unbinding pathway. The contour plots for FHI-6 and HI-6 13 drugs’ unbinding are markedly different from each other (compare between numerals 1 and 2 and the green star symbols in Figure 5). Calculated electrostatic energies between aromatic residues and the drug molecules (see Table 1) clearly indicate that they are always higher for FHI-6 than its non-fluorinated counterpart. This is also reflected in the cation-π order parameters in that they are higher for FHI-6 than that prevailing in AChE•HI-6 complex. As a result, both the CASP and PASP rings of the FHI-6 molecule are found to be strongly held in the CAS and in the PAS region of enzyme, respectively. On the other hand, the CASP ring of HI-6 is also found to interact with the PAS region of the enzyme in the bound state. 13 This causes the minimum of 1D PMF profile for HI-6 to be located at a higher distance along the COM separation CV, close to the PAS region of the enzyme. Whereas, the 1D PMF profile for FHI-6 is found to be located at much lower distance along the COM separation CV, closer to the CAS region. Put in other words, FHI-6 is more tightly bound inside the binding gorge than HI-6. This is also reflected in the calculated microkinetics of unbinding transition (see below) in that the unbinding rate constant of FHI-6 is about 10 times lower than that of HI-6.

3.5. Hydrogen Bond and Water Bridge Dynamics In water, the major driving forces for protein folding are the maximizing the interaction of polar residues with water and minimizing the interaction of hydrophobic residues with the water. 61 This phenomenon leads to the dielectric heterogeneity i.e. interior of the protein is quite hydrophobic in contrast to the hydrophilic surface. 10–14,62 Recent studies show that the dielectric constant of bulk water is different from small water clusters. 63,64 The dipole moment of water molecule is different in different environment e.g. in bulk, in the surface of a protein and in the interior of a protein. Thus, the water molecule buried deep inside the hydrophobic pocket e.g., gorge of the AChE shows reduced dielectric constant that that of the bulk. This leads to the stabilization of electrostatic interaction i.e. HB between the protein and drug. 21



Figure 6: Time dependence of protein - FHI-6 HB lifetime correlation function, CHB (t) in polarizable water medium at the basins and barrier along the gorge pathway are shown. Hydrogen bonds between all possible D-A pairs are considered. Major amino acid residues of AChE protein responsible for the HB formation at various basins and barrier are reported in Table 2. Single exponential fitting parameters (time constant and life time) are also reported in Table 2.

It is known that the breaking of HB in the hydrophobic environment requires 1.2 kcal/mol more energy. 65 During the unbinding transition of FHI-6 drug, hydrophobic HBs must break to facilitate the transition. Therefore, the lifetime and binding kinetics of the AChE•FHI-6 complex are dictated by the hydrophobic HBs due to their higher thermodynamic stability. However, water molecules present in the hydrophobic gorge act as lubricant by inserting between broken protein-drug HBs and forming two HBs (in place of it) with the protein residue and the drug molecule simultaneously using two of its hydrogen atoms. This helps in facilitating the unbinding of the drug. 66,67 Hence, lubricating water molecules decrease the energy barrier during the liberation process. Thus, the unbinding transition of drug is controlled by the dynamical balance between making and breaking of HBs and WBs. In view of this, additional equilibrium MD simulations are performed in each of the basins (1 and 2) and barrier (1) along the unbinding pathways (see Figure 2) to identify the protein residues involved in HB and WB interactions. All the equilibrium MD simulations are carried out in NVT ensemble for 20 ns. It is worthwhile to mention that the last 10 ns of the trajectory are analyzed. The coordinates of the solvated AChE•FHI-6 system are 22


Page 22 of 44

Page 23 of 44 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


saved in every 1 ps during the simulation. The lifetime correlation functions for different HB and WB (cf. Eqs. 6 and 7) are calculated to investigate the microscopic aspects of hydration dynamics during the unbinding of FHI-6 drug. The probability that a random pair that is either hydrogen bonded or water bridged at time ’0’ is still bonded at time t, regardless of whether or not the specific HB or WB was broken at some intermediate times is described by the lifetime correlation function. The rearrangement of the underlying WB or HB network can be tracked by the decay of HB and WB lifetime correlations functions. Thus, the hydration dynamics during the unbinding of the drug can be effectively monitored at the molecular level by studying the lifetime correlations functions, C HB (t) and C WB (t) The variations of the HB and WB lifetime correlation functions, namely, C HB (t) and C WB (t) are shown in Figures 6 and 7, respectively. The nature of the decay of C HB (t)) in all the basins and barrier is single exponential. On the other hand, C WB (t) shows multi-exponential decay in majority of instances, except for residues that show single exponential decay; V282 in basin1, H447 in barrier 1 and E285 in basin 2. Exponential fitting parameters (time constants and amplitudes) for the time correlation functions of HB and WB in polarizable water medium at various basins and barriers along with the major amino acid residues of the protein responsible for the HB and WB formations are indicated in Table 2 and Table 3, respectively.

Table 2. Single exponential fitting parameters for the time correlation functions of the hydrogen bond between AChE and FHI-6 at various basins and barriers along the drug unbinding pathway. The protein residues involved in the hydrogen bond formations are also indicated. Single word abbreviations of amino acids are used.

Basin/ Barrier

Protein residue

Basin 1 Barrier 1 Basin 2

G120, E285, S298, Y337, H447 D74, T75, V282, E285, W286„S298, Y341 T75, E285, V282, H287, S298

23


Time constant (ps) Lifetime (ps) 250.1 5802.2 3664.7

82.1 539.7 486.5


Page 24 of 44

It is clear from Figures 6 and 7 and Tables 2 and 3, the nature of the time correlations functions not only depend on the basins and barrier but also on the participating residues. The time constant (250.1 ps) and lifetime (82.1 ps) of HB in basin 1 is short, but the same of WBs are higher. Specially, the lifetime of WB involving E285 is extraordinarily higher (2330.8 ps). So the basin 1 is predominantly stabilized by WB bonding. On the other hand, time constants and life times of HBs are much higher in barrier 1 and basin 2 in comparison to the same of WBs. The WB lifetime correlations function in barrier 1 and basin 2 decays much faster than the corresponding HBs and their lifetime is below 100 ps for all the residues involved in WB formations. Thus the barrier 1 and basin 2 are predominantly stabilized by HBs.

Table 3. Exponential fitting parameters (time constants, amplitudes and life time) for the time correlation functions of the water bridge between AChE and FHI-6 drug in polarizable water medium at various basins and barriers along the drug unbinding pathway. The major protein residues involved in the water bridge formations are indicated. Single word abbreviations of amino acids are used.

Basin/ Barrier Protein residue Time constants (ps) Amplitude (%) Lifetime (ps)

Basin 1

E202 V282 E285

73.6, 7956.7 2889.7 5664.3, 5664.2

29, 71 100 50, 50

179.8; 377.2 2330.8

Barrier 1

V282 E285 S298 Y337 H447

249.7, 28.7 654.6, 66.9 1.0, 59.1 124.3, 506.6 28.7

72, 28 81, 19 88, 12 27, 73 100

45.9 99.0 21.8 59.0 16.7

Barrier 1

V282 E285 S298

79.4, 1.0 65.4 189.8, 189.7

53, 47 100 50, 50

43.7 38.2 80.9

24


Page 25 of 44 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 7: Time dependence of protein - FHI-6 WB lifetime correlation function, CWB (t) in polarizable water medium for the major amino acid residues of AChE responsible for the WB formation at various basins and barriers along the oxime drug unbinding pathway are shown. Exponential fitting parameters (time constants, amplitudes and life time) are also reported in Table 3.

25



3.6. Microkinetics of FHI-6 Unbinding from the AChE Gorge A combined Mtd and HD method 13,53 is used to calculate the rate constant of basin1 → basin2 unbinding transition. In this combined approach, one first conducts a single MtD simulation to find a customized biasing potential, which is then "frozen" and used to exploit the strong MD acceleration in the susequent HD runs. While customizing the MtD biasing potential (c.f. dashed line in Figure 2) that is to be used in HD runs, we have ensured that it vanishes at the separatrix between the two basins. An on-the-fly combination of both techniques has recently been suggested by Tiwary and Parrinello 68 in their pioneering work to obtain dynamics from metadynamics. Their "on-the-fly" method is although computationally less expensive than our present "two-step" approach, it runs the risks of depositing bias in the transition state region, requiring further rigorous assessment of the reliability of the discovered dynamics. 69 On the other hand, in the two-state approach, we first obtain a converged free-energy profile along a chosen CV by MtD run. A non-negative biasing potential, V bias (c.f. dashed line in Figure 2) is then constructed by collecting the accumulated bias up to a suitably chosen time of the converged MtD run such that the constructed potential lies well below the transition state region. V bias is then used in subsequent auxiliary HD runs 52 in order to extract kinetics following the principles of transition state theory. Thus although the present two separate calculations are computationally more expensive, they together safeguard the trustworthy of calculated unbinding kinetics (free from the risk of adding bias in transition state region and hysteresis in the free-energy profile) 69 along the chosen CV over converged well-defined free energy barrier. To reconstruct unbiased state-to-state time evolution we have used a suitable renormalization of time through a boost factor, f ac (see Eq. 8). 13 Several HD simulations with different starting configurations from the bottom of basin1 are performed. The fraction of reactants that converted to products per unit time can bs used to calculate the forward microkinetic rate constant. However, this measurement can be flawed because of the implicit assumption that each protein-drug configuration which crossed the transition state will end up in the product side and contribute to the forward rate. Whereas in reality the population in the dividing planes can be influenced by the cation-π interactions or hydration dynamics and thus play important role on the fate of the

26


Page 26 of 44

Page 27 of 44 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


measurement of the rate constants. 70 Therefore, transmission coefficient, κ is used to correct the HD rate in order to account for the possibility of re-crossings at the barrier top.

Figure 8: Evaluation of bias rate constant (kbias ) from basin 1 to basin 2 by employing HD simulation. The FPT is shown in histogram and its single exponential decay (red colour line) is also shown. The plateau value to calculate the transmission co-efficient is also shown in the inset.

Total number of escape trajectories in the forward direction are identified from the 60 HD runs and are shown in Table 4. The configurations in the COM separation coordinate that crosses the barrier top, ξ # with positive velocity may rapidly re-crosses the barrier before settling in the product basin are also included in the escape trajectories. Then a histogram of the first passage time (FPT) of the barrier crossing can be constructed to calculate the bias rate constant, kbias . 53 It is clear from the Figure 8 that the decay profile of the FPT follow a first order rate equation and decay constant can be equated to kbias . In order to correct the HD rate from solvent induced re-crossing events 70 at the barrier top, we have calculated the transmission coefficient, κ following a simple ratio,

κ(t) =

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

27


(9)


Page 28 of 44

Here, NR→P , NP →R , NR→R and NP →P are the number of occurrences of trajectories that had started from reactant well and end up in product well, from product well to reactant well, from reactant to reactant well and from product to product well, respectively. It is worthwhile to mention that NR→P makes positive contribution to κ and NP →R makes negative contribution while NR→R and NP →P are unproductive. The calculation of κ(t) is shown in Figure 8 inset. κ(t) can be ssen to vary from 0 to 1 and in accordance with Onsager regression hypothesis, it attains a plateau value at a time which is longer than the typical time of molecular motions but much shorter than the reaction rate constant i.e. κ = limt→∞ κ(t). Finally, the rate constant is given as,

k = κ × (kbias /fac ) .

(10)

Table 4. Rate constant, number of escape events, 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 of FHI-6 from the active site gorge of the AChE.

Transition

Number of escape events

k bias (s-1 )

f ac

κ (t→ ∞)

k(s-1 )

Basin 1→2

792

2.18 × 108

61

0.011

3.94 × 104

Rate constant (k ), number of escape events, the biased rate constants (kbias ) as obtained through the HD runs, the boost factor (fac ) and the long-time limit [κ(t → ∞)] of the transmission coefficient (κ) for the unbinding transitions of FHI-6 from the active site gorge of the AChE are shown in Table 4. A comparison with calculated rate constants for HI-6 13 in MFP/TIP3P water medium reveals that the unbinding rate constant of FHI-6 is about 10 times lower. This indicates that the introduction of one fluorine atom to the drug can provide extra stability toward protein-drug ineractions in the bound state. This also corroborates well with the calculated PMF profile. This additional stability, caused by cation-π, HB and WB interactions can act as a possible kinetic trap in the course of unbinding of the 28


Page 29 of 44 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


fluorinated drug and thus could be more effective towards inhibited enzyme repair than its non-fluorinated counterpart.

4. Conclusions All-atom MD simulations are conducted to investigate the effect of fluorination of bispyridinium oxime drug, HI-6 on the dynamic mechanism of drug’s unbinding from the active site gorge of the nerve enzyme, AChE. Effective polarization in mean-field polarizable model of TIP3P water is used in all the MD simulations. The unbinding process consisted of basins and barrier is found to be strongly influenced by cation-π interactions, WB, and HB bridge formation/breaking between the protein and the FHI-6 molecule. Major enhancements and modifications in the cation-π interactions and in intermolecular HB and WB network are observed with the introduction of one fluorine atom in the oxime drug. Enhance sampling method, namely, WT-MtD simulation is used to calculate 1D and 2D PMF profiles for the unbinding of fluorinated drug from the AChE gorge. Cation-π meta-coordinate is used as another CV along with COM separation CV to generate the 2D PMF profile using a reweighting scheme. Both the PMF profiles are markedly different from that of the non-fluorinated drug. 13 Following quantum mechanical method, it is also observed that FHI-6 interacts more strongly with the aromatic residues of the enzyme in comparison to HI-6. This is also reflected on the enhanced cation-π order parameter of FHI-6 interaction, calculated by post-processing the unbinding trajectory based on a geometry based criterion. Cation-π interaction, HB and WB dynamics of the protein-drug system finally dictates the binding profile of the drug molecule. An introduction of fluorine atom to HI-6 causes the CASP and PASP rings of the fluorinated drug (FHI-6) to interact with the CAS and PAS region of the enzyme, respectively. On the other hand, the CASP ring of HI-6, in addition to making cation-π interaction with aromatic residues in the CAS region is also found to interact with those in the PAS region of the enzyme. This causes the bottom of the 1D unbinding free energy profile to be shifted towards PAS region of the enzyme. 13 Whereas, as we have shown here for FHI29



6, the bottom of the said free energy profile is located towards the CAS region of the enzyme because of the more harmonious cation-π interactions: CASP ring of FHI-6 with CAS region of the enzyme and PASP ring with the PAS region. This result into FHI-6 being more firmly bound in the gorge, as is also evidenced from reduced rate constant of unbinding transitions, measured through a combination of metadynamics and hyperdynamics simulations. This may suggest that FHI-6 (in addition to its inherent BBB penetration property) will be more effective in repair of the inhibited enzyme than its non-fluorinated counterpart.

Supporting Information Additional figures in support of convergence test of well-tempered metadynamics run. The material is available free of charge via the Internet at http://pubs.acs.org. Acknowledgement The authors would like to thank computer centre, BARC for providing the ANUPAM parallel computational facility. Dr(s). A. K. Samanta and P. D. Naik are gratefully acknowledged for the constant encouragement. The work was supported by DAE under project XII-N-R&D02.04/Theoretical & Computational Chemistry of Complex Systems. The authors declare no conflict of interest.

References (1) Sarter, M.; Parikh, V. Choline transporters, cholinergic transmission and cognition. Nat. Rev. Neurosci. 2005, 6, 48-56. (2) Quinn, D. M. Acetylcholinesterase: enzyme structure, reaction dynamics, and virtual transition states. Chem. Rev. 1987, 87, 955-79. (3) Hedstrom, L. Serine protease mechanism and specificity. Chem. Rev. 2002, 102, 4501-24. (4) Taylor, P.; Radic, Z. The cholinesterases: from genes to proteins. Annu. Rev. Pharmaco. Toxico. 1994, 34, 281-320. (5) Tripathi, A.; Srivastva, U.C. Acetylcholinesterase: A versatile enzyme of nervous system. Ann.. Neuro Sci. 2008, 15, 106-111. (6) Sussman, J. L.; Harel, M.; Frolow, F.; Oefner, C.; Goldman, A.; Toker, L.; Silman, I. Atomic structure of acetylcholinesterase from torpedo californica: A prototypic acetylcholine-binding protein. Science 1991, 253, 872-879.

30


Page 30 of 44

Page 31 of 44 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


(7) Bourne, Y.; Taylor, P.; Marchot, P. Acetylcholinesterase inhibition by fasciculin: crystal structure of the complex. Cell 1995, 83, 503-512. (8) Radic, Z.; Gibney, G.; Kawamoto, S.; MacPhee-Quigley, K.; Bongiorno, C.; Taylor, P. Expression of recombinant acetylcholinesterase in a baculovirus system: Kinetic properties of glutamate 199 mutants. Biochemistry 1992, 31, 9760-7. (9) Desiraju, G. R.; Steiner, T. The weak hydrogen bond in structural chemistry and biology. Oxford University Press: New York, 1999. (10) Pal, S. K.; Zewail, A. H. Dynamics of water in biological recognition. Chem. Rev. 2004, 104, 2099-2123. (11) 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. (12) Pathak, A. K.; Bandyopadhyay, T. Unbinding of fluorinated oxime drug from the AChE gorge in polarizable water: A well-tempered metadynamics study. Phys. Chem. Chem. Phys. 2017, 19, 55605569. (13) Pathak, A. K.; Bandyopadhyay, T. Protein-drug interactions with effective polarization in polarizable water: Oxime unbinding from AChE gorge J. Phys. Chem. B 2015, 119, 14460-14471. (14) Tiwary, P.; Limongelli, V.; Salvalaglio, M.; Parrinello, M. Kinetics of protein-ligand unbinding: Predicting pathways, rates, and rate limiting steps. Proc. Natl. Acad. Sci. U. S. A. 2015, 112, E386-E391. (15) Newmark, J. The birth of nerve agent warfare: Lessons from Syed Abbas foroutan. Neurology 2004, 62, 1590-1596. (16) Eddleston, M.; Buckley, N. A.; Eyer, P.; Dawson, A.-H. Management of acute organophosphorus pesticide poisoning. Lancet 2008, 371, 597-607. (17) Butt,A. M.; Jones, H. C.; Abbott, N. J. Electrical resistance across the blood-brain barrier in anaesthetized rats: A developmental study. J. Physiol. 1990, 429, 47-62. (18) Ballabh, P. ; Bruan, A.; Nedergaard, M. The blood-brain barrier: an overview: structure, regulation, and clinical implications. Neurobiol. Dis. 2004, 16, 1-13. (19) Abbott, N. J.; Ronnback, L. Hansson, E. Astrocyte-endothelial interactions at the blood-brain barrier. Nat. Rev. Neurosci. 2006, 7, 41-53. (20) Ekstrom, F., Pang, Y.P., Boman, M., Artursson, E., Akfur, C., Borjegren, S. Crystal structures of acetylcholinesterase 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. (21) Mercey, G.; Verdelet, T.; Renou, J.; Kliachyna, M.; Baati, R.; Nachon, F.; Jean, L.; Renard, P. Y. Reactivators of acetylcholinesterase inhibited by organophosphorus nerve agents. Acc Chem Res. 2012, 45, 756-66. (22) Sakurada, K.; Matsubara, K.; Shimizu, K.; Shiono, H.; Seto, Y.; Tsuge, K.; Yoshino, M.; Sakai, I.; Mukoyama, H.; Takatori, T. Pralidoxime iodide (2-PAM) penetrates across the blood-brain barrier. Neurochem. Res. 2003, 28, 1401-1407. (23) Cook, A. M.; Mieure, K. D.; Owen, R. D.; Pesaturo, A. B.; Hatton, J. Intracerebroventricular administration of drugs. J. Pharmacotherapy 2009, 29, 832-845. (24) Bradley, W. G., Jr. R-guided focused ultrasound: A potentially disruptive technology. J. Am. Coll. Radiol. 2009, 6, 510-513.

31



(25) Neuwelt, E. A.; Frenkel, E. P.; Diehl, J. T.;Maravilla, K. R.; Vu, L. H.; Clark, W. K.; Rapoport, S. I.; Barnett, P. A.; Hill, S. A.; Lewis, S. E.; Ehle, A. L.; Beyer, C. W., Jr.; Moore, R. J. Osmotic blood-brain barrier disruption: A new means of increasing chemotherapeutic agent delivery. Trans. Am. Neurol. Assoc. 1979, 104, 256-260. (26) Carman, A. J.; Mills, J. H.; Krenz, A.; Kim, D. G.; Bynoe, M. S. Adenosine receptor signaling modulates permeability of the blood-brain barrier. J. Neurosci. 2011, 31, 13272-13280. (27) Radic, Z.; Sit, R. K.; Kovarik, Z.; Berend, S.; Garcia, E.; Zhang, L.; Amitai, G.; Green, C.; Radic, B.; Fokin, V. V.; Sharpless, K. B.; Taylor, P. Refinement of structural leads for centrally acting oxime reactivators of phosphylated cholinesterases. J. Biol. Chem. 2012, 287, 11798-11809. (28) Gorecki, L.; Korabecny, J.; Musilek, K.; Malinak, D.; Nepovimova, E.; Dolezal, R.; Jun, D.; Soukup, O.; Kuca, K. SAR study to find optimal cholinesterase reactivator against organophosphorous nerve agents and pesticides. Arch Toxicol. 2016, 90, 2831-2859. (29) Kovarik, Z.;, Macek, N.; Sit, R. K.; Radic, Z.; Fokin, V. V.; Sharpless, K. B.; Taylor, P. Centrally acting oximes in reactivation of tabun-phosphoramidated AChE. Chem. Biol. Interact. 2013, 203,77-80. (30) Jeong, H. C.; Kang, N. S.; Park, N. J.; Yum, E. K.; Jung, Y. S. Reactivation potency of fluorinated pyridinium oximes for acetylcholinesterases inhibited by paraoxon organophosphorus agent. Bioorg. Med. Chem. Lett. 2009, 19, 1214-1217. (31) Jeong, H. C.; Park, N. J.; Chae, C. H.; Musilek, K.; Kassa, J.; Kuca, K.; Jung, Y. S. Fluorinated pyridinium oximes as potential reactivators for acetylcholinesterases inhibited by paraoxon organophosphorus agent. Bioorg. Med. Chem. 2009, 17, 6213-6217. (32) 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. (33) Leontyev, I. V.; Stuchebrukhov, A. A. Polarizable molecular interactions in condensed phase and their equivalent nonpolarizable models. J. Chem. Phys. 2014, 141, 014103. (34) Barducci, A.; Bussi, G.; Parrinello, M. Well-tempered metadynamics: A smoothly converging and tunable free-energy method. Phys. Rev. Lett. 2008, 100, 020603-020606. (35) 34. Morris, G.M., Huey, R., Lindstrom, W., Sanner, M.F., Belew, R.K.,. Goodsell, D.S. and Olson, A.J. AutoDock4 and autoDockTools4: Automated docking with selective receptor flexiblity. J. Comput. Chem. 2009, 16, 2785-91. (36) 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. (37) 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: Struct., Funct., Genet. 2006, 65, 712-725. (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) 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. General atomic and molecular electronic structure system. J. Comput. Chem. 1993, 14, 1347-1363. (40) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648-5652.

32


Page 32 of 44

Page 33 of 44 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


(41) Nose, S. A molecular dynamics method for simulations in the canonical ensemble. J. Chem. Phys. 1984, 81, 511. (42) Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A. 1985, 31, 1695. (43) 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. (44) 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. (45) Barducci, A.; Bonomi, M.; Parrinello, M. Metadynamics. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2011, 1, 826-843. (46) 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. Comput. Phys. Commun. 2009, 180, 1961-1972. (47) Bonomi, M.; Barducci, A.; Parrinello, M. Reconstructing the equilibrium Boltzmann distribution from well-tempered metadynamics. J. Comput. Chem. 2009, 30, 1615-1621. (48) Minoux, H.; Chipot, C. Cation-π interactions in proteins: Can simple models provide an accurate description? J. Am. Chem. Soc. 1999, 121, 10366-10372. (49) 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. (50) Pathak, A. K.; Bandyopadhyay, T. Ortho-7 bound to the active-site gorge of free and OP-conjugated acetylcholinesterase: Cation-π interactions. Biopolymers, 2016, 105, 10-20. (51) Su, P.; Li, H. Energy decomposition analysis of covalent bonds and intermolecular interactions. J. Chem.Phys. 2009, 131, 014102. (52) Voter, A. F. A method for accelerating the molecular dynamics simulation of infrequent events. J. Chem. Phys. 1997, 106, 4665-4677. (53) Schneider, J.; Reuter, K. Efficient calculation of microscopic dissolution rate constants: The aspirinwater interface. J. Phys. Chem. Lett. 2014, 5, 3859-3862. (54) Tiwary, P.; Parrinello, M. A time-independent free Energy estimator for metadynamics. J. Phys. Chem. B 2015, 119, 736-742. (55) Shen, T.; Tai, K.; Henchman, R. H.; McCammon, J. A. Molecular dynamics of acetylcholinesterase. Acc. Chem. Res. 2002, 35, 332-340. (56) Bai, F.; Xu, Y.; Chen, J.; Liu, Q.; Gu, J.; Wang, X.; Ma, J.; Li, H. ; Onuchic, J. N. ; Jiang. H. Free energy landscape for the binding process of Huperzine A to acetylcholinesterase. Proc. Natl. Acad. Sci. 2013, 110, 4273-8. (57) Xu, Y; Shen, J.; Luo, X. ; Silman, I.; Sussman, J. L.; Chen, K. ; Jiang. H. How does huperzine A enter and leave the binding gorge of acetylcholinesterase? Steered molecular dynamics simulations. J. Am. Chem. Soc. 2003, 125, 11340-9. (58) Kesharwani, M. K.; Ganguly, B.; Das, A.; Bandyopadhyay, T. Differential binding of bispyridinium oxime drugs with acetylcholinesterase. Acta Pharmacol Sin. 2010, 31, 313-28.

33



(59) Dougherty, D. A. Cation-π interactions in chemistry and biology: A new view of benzene, Phe, Tyr, and Trp. Science 1996, 271, 163-168. (60) 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. (61) Rose, G.; Fleming, P.; Banavar, J.; Maritan, A. A backbone-based theory of protein folding. Proc.Natl. Acad. Sci. (USA.) 2006, 103, 16623-33. (62) Isom, D.G.; Cannon, B.R.; Castaneda, C.A.; Robinson;A.; Garcia-Moreno, E.B. High tolerance for ionizable residues in the hydrophobic interior of proteins . Proc. Natl. Acad. Sci. (USA) 2008, 105, 17784-88. (63) Pathak, A. K.; Samanta, A. K. Dielectric constant of solvents at nano to bulk regimes: Linear response theory and statistical mechanics-based approaches. Mole. Phys. 2017, 115, 2508-2514. (64) Pathak, A. K. A theoretical study on microhydration of chromate dianion: On the number of water molecules necessary to mimic the bulk. Mole. Phys. 2015, 113, 3345-3352. (65) 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. (66) 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. (67) 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. (68) Tiwary, P.; Parrinello, M. From metadynamics to dynamics. Phys. Rev. Lett. 2013, 111, 230602-230605. (69) Salvalaglio, M.; Tiwary, P.; Parrinello, M. Assessing the reliability of the dynamics reconstructed from metadynamics. J. Chem. Theory Comput. 2014, 10, 1420-1425. (70) Ciccotti, G.; Ferrario, M.; Hynes, J. T. Dynamics of ion pair interconversion in a polar solvent. J. Chem. Phys. 1990, 93, 7137-7147.

34


Page 34 of 44

Page 35 of 44 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


TO OC Graph hic



Table of Content


Page 36 of 44

Page 37 of 44 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 150x176mm (600 x 600 DPI)



Figure 2 254x190mm (96 x 96 DPI)


Page 38 of 44

Page 39 of 44 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 153x131mm (600 x 600 DPI)



Figure 4 189x184mm (600 x 600 DPI)


Page 40 of 44

Page 41 of 44 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 296x209mm (300 x 300 DPI)



Figure 6 79x63mm (300 x 300 DPI)


Page 42 of 44

Page 43 of 44 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 7 209x243mm (600 x 600 DPI)



Figure 8 271x192mm (150 x 150 DPI)


Page 44 of 44

Dynamic Mechanism of a Fluorinated Oxime Reactivator Unbinding

Recommend Documents